Thesis RGdeJong

3DP
A control system model for combined DP

station keeping and active roll reduction
Technische Universiteit Delft
R.G. de Jong
3DP
A control system model for combined DP
station keeping and active roll reduction
by
R.G. de Jong
to obtain the degree of Master of Science
at the Delft University of Technology,
to be defended publicly on Thursday July 5, 2018 at 14.00
Student number: 4195647

Project duration: October 2, 2017 – 5 July 2018
Thesis committee: Prof. dr. ir. A. P. van ’t Veer, TU Delft, chair
Dr. ir. P. R. Wellens, TU Delft, daily supervisor
Prof. dr. ir. J. W. van Wingerden, TU Delft
Dr. ir. M. Godjevac, TU Delft
Ir. T. G. Vos, Boskalis
Ir. R. Beindorff, Boskalis
An electronic version of this thesis is available at http://repository.tudelft.nl/

Preface
This thesis is the final deliverable in my pursuit of a master’s degree in Ship Hydromechanics at Delft
University of Technology. During the thesis research, I tried to develop a control system that enables
a ship to reduce its roll motion during Dynamic Positioning operations. The research combined my
interest in ship motions, thrusters and control systems. I enjoyed working on it for the past 8 months
and I am happy with the end result. I could not have done it without the help and support of several
people and I would like to thank them for their contributions.
First of all, I would like to thank Thijs Vos and Ruud Beindorff for their unrelenting enthusiasm, support
and for giving me the freedom to conduct the research according my own ideas.
Furthermore, Peter Wellens and Helio Bailly Guimaraes, thanks for teaching and inspiring me during
the final years of my education. You somehow always succeeded in feeding my curiosity and enthusi-
asm.
Last but not least, I would like to thank my family, Barbara, my friends and in particular Roel and
Casper for the good times both offshore and onshore.
R.G. de Jong
Papendrecht, June 2018
iii
Abstract
During DP operations in beam seas, the cable-lay vessel Ndurance experiences significant roll motions
that decrease the operability. This research aims to evaluate the possibility of reducing the ship roll
motion during DP operations by active control of the thrusters.
A vessel time domain model is constructed based on Cummins’ equation. The convolution term is
evaluated by using a state-space representation. Since no viscous effects are included when diffrac-
tion analysis is used, the viscous roll damping due to eddy making and skin friction are calculated in the
frequency domain by using empirical models. The viscous roll damping term in the time domain model
is subsequently determined by tuning a quadratic roll velocity coefficient in the time domain model until
the roll RAO matches the frequency domain roll RAO with viscous damping included.
A dynamic thruster model has been applied to the thrusters installed at the Ndurance. An empiri-
cal estimate of the added rotational inertia of the entrained water in the thruster has been included
for conservatism. To include the effect of both negative and positive thruster inflow velocities, a four-
quadrant model has been used to calculate the thrust and torque characteristics of the thrusters.
A DP system model has been constructed consisting out of a Kalman filter, PID controller and a thrust
allocation algorithm. The DP system model is merged with the vessel time domain model and is sub-
sequently subjected to first- and second-order wave forces, current forces and wind forces. The model
results are compared to time domain DP simulations conducted by MARIN.
In order to achieve combined DP station keeping and active roll reduction a new control model is
developed. The model is referred to as 3DP (3-dimensional DP). The control system model is based
on a hierarchical controller structure, consisting of high-level motion controllers and low-level azimuth
angle and shaft speed controllers. The high-level motion controller consists out of a combination of a
roll controller and a conventional DP controller. Both controller commands are merged by application
of proposed allocation matrices. The merged controller commands are subsequently used as setpoint
for the low-level shaft speed controller. The azimuth angle of the thrusters used to counteract the en-
vironmental loads in the surge direction is controlled by a proportional controller. The azimuth angle
controller is implemented to achieve both maximum roll reduction and sufficient station keeping ability.
Multiple simulations have been carried out with both the DP and the 3DP control model. From the
results of these simulations can be concluded that the proposed 3DP control model enables the pos-
sibility of thruster induced roll reduction. The maximum roll reduction is achieved during sea states
with low to moderate significant wave heights and wave periods in the natural frequency range. Dur-
ing these conditions, the roll reduction is around 0.5∘ roll RMS. From this result is concluded that the
3DP model is able to decrease the roll motion to a significant extent. The effect of thruster induced
roll reduction onto the station keeping performance of the ship is also investigated. From these results
can be concluded that the DP footprint of the vessel increases with a maximum of 1 meter when the
3DP model is applied. The DP station keeping capability results showed that the vessel is limited to a
beam current velocity of 0.9 m/s. The power consumption of the thrusters is calculated according the
simulation results. It is concluded that the 3DP control system consumes roughly 10 times more power
than the conventional DP system. Also, the workability increase has been calculated when the 3DP
model is used. For beam waves this resulted in a yearly workability increase of 5.2%. The workability
increase averaged over every environmental direction is 2.0%.
It is recommended that the 3DP control model is used as a back-up instrument to reduce the roll motion
of the vessel during critical offshore operations, like tool over-boardings or cable pull-ins. In this way,
the operational use of the system is limited to short periods of time and the wear and tear and increased
fuel consumption of the thrusters is of minimal nature.
v
Nomenclature
COG - Centre of Gravity
DOF - Degrees of Freedom
DP - Dynamic Positioning
IRF - Impulse Response Function
QTF - Quadratic Transfer Function
RAO - Response Amplitude Operator
RMS - Root Mean Square
RPM - Rounds Per Minute
RRR - Relative Roll Reduction
SS - State-Space
TF - Transfer Function
VFD - Variable Frequency Drive
WD - Water Depth
vii
List of Symbols
Greek Symbols
𝜂 vessel motion vector [m] or [rad]

𝜂 root mean square values of motion vector [m] or [rad]
𝜖 phase angle [rad]
𝜙 roll amplitude [rad] or [deg]
𝜓 yaw amplitude [rad] or [deg]
𝜃 pitch amplitude [rad] or [deg]
𝜔 wave frequency or propeller angular velocity [rad/s]
𝜔 cut-off wave frequency [rad/s]
𝜅 state-space convolution representation [-]
𝜉 wave amplitude [m]
𝜏 motion controller command [RPM]
𝜏 roll controller command [RPM]
𝜏 DP controller command [RPM]
𝜏 azimuth controller command [deg]
𝜏 shaft speed controller command [Nm]
𝛽 hydrodynamic pitch angle [deg]
𝜇 environmental direction [deg]
𝜌 sea water density [kg/m ]
𝜌 air density [kg/m ]
Roman Symbols
A infinite frequency added mass matrix [kg] or [kgm ]

A state-space system matrix [-]
A(𝜔) hydrodynamic added mass matrix [kg] or [kgm ]
𝐴 lateral wind area [m ]
𝐴 frontal wind area [m ]
𝑏 required forces and moments in the horizontal plane [kN] or [kNm]
B state-space input matrix [-]
B(𝜔) hydrodynamic damping matrix [kg/s] or [kgm /s]
𝐵 bilge keel roll damping coefficient [kgm/s rad]
𝐵 eddy making roll damping coefficient [kgm/s rad]
𝐵 equivalent roll damping coefficient [kgm/s rad]
𝐵 skin friction roll damping coefficient [kgm/s rad]
𝐵 lift roll damping coefficient [kgm/s rad]
𝐵 hydrodynamic damping at the cut-off frequency [kgm/s rad]
𝐵 quadratic viscous roll damping term [kgm/s rad]
𝐵 wave making roll damping [kgm/s rad]
C ship hydrostatic stiffness matrix [N/m] or [Nm/rad]
C state-space output matrix [-]
𝑐(𝑥) non-linear constraints vector [kN]
𝐶(𝜇) wind/current load coefficient [-]
𝑐 viscous quadratic damping tuning parameter [-]
𝐶 entrained water fitting parameter [-]
𝐶 four-quadrant torque coefficient [-]
ix
x List of Symbols
𝐶 four-quadrant thrust coefficient [-]

𝐷 propeller diameter [m]
𝑒 setpoint error [m] or [deg]
𝐸𝐴𝑅 propeller expanded area ratio [-]
F(𝜔) first-order wave excitation force amplitude vector [N] or [Nm]
𝐹 (x) thrust allocation objective function [-]
𝐻 significant wave height [m]
𝐼 added rotational inertia of entrained water in propeller [kg m ]
𝐼 propeller rotational inertia [kg m ]
𝐼 propeller shaft rotational inertia [kg m ]
𝐼 rotor rotational inertia [kg m ]
𝐼 total thruster system rotational inertia [kg m ]
𝐾 empirical coefficient entrained water added rotational inertia [-]
𝐾 proportional controller coefficient [-]
𝐾 integral controller coefficient [-]
𝐾 derivative controller coefficient [-]
𝑘 gearbox ratio [-]
𝐿 vessel length between perpendiculars [m]
M ship inertia matrix [kg] or [kgm ]
𝑁 total number of frequencies [-]
𝑁 total number of thrusters [-]
𝑃/𝐷 propeller pitch-diameter ratio [-]
𝑄 commanded motor torque [Nm]
𝑄 motor torque [Nm]
𝑄 propeller torque [Nm]
𝑄 commanded motor torque [Nm]
𝑅 maximum DP excursion in the horizontal plane [m] or [deg]
𝑆 (𝜔) motion response spectrum [m s] or [rad s]
𝑆 (𝜔) wave spectrum [m s]
𝑡 time vector [s]
𝑇 vessel draft [m]
𝑇 cut-off scaling time constant [-]
𝑇 thruster time constant [s]
𝑇 maximum thrust magnitude [kN]
𝑇 wave peak period [s]
𝑣 propeller advance velocity [m/s]
𝑊 workability when 3DP system is engaged [%]
𝑊 workability when DP system is engaged [%]
𝑊 workability increase [%]
𝑥 vessel surge displacement [m]
𝑦 vessel sway displacement [m]
𝑦 maximum value in model result range [-]
𝑦 minimum value in model result range [-]
𝑧 vessel heave displacement [m]
𝑍 number of propeller blades [-]
Contents
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Research objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 Research approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.5 Report structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Ship Model 5
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Ship motions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Frequency domain model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Time domain model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3.1 The convolution term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3.2 IRF in the time domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.3 IRF in the frequency domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.4 Transfer function fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.5 Time domain model verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Viscous roll damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.1 Ikeda viscous damping model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.2 Frequency domain model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.3 Time domain model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3 Thruster Model 23
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2 Dynamic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.1 Entrained water added rotational moment of inertia . . . . . . . . . . . . . . . . . 24
3.2.2 Four-quadrant model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Model verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Thruster curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Transient response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.4 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4 DP System Model 31
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2 DP system configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3 Environmental forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3.1 Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.3.2 Wind. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.3 Current . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.4 Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.5 DP controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.6 Thrust allocation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.7 Model verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.8 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
xi
xii Contents
5 3DP System Model 45

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2 Conventional DP control system vs. 3DP control system . . . . . . . . . . . . . . . . . . 45
5.3 Allocation model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.3.1 Roll reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.3.2 Station keeping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.4 Controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.4.1 Roll reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.4.2 Station keeping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4.3 Azimuth angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4.4 Shaft speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.5 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6 Results 57
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.2 Roll reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.3 DP footprint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.4 DP capability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.5 Thruster power consumption. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.6 Workability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
7 Conclusions and Recommendations 67
7.1 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7.2 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
A Appendix 71
A.1 Vessel particulars. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
A.2 Ramp function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
B Appendix 73
B.1 Thruster and electrical drive data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
B.2 Four-quadrant Fourier coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
B.3 Thrust deduction factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
C Appendix 77
C.1 Wind and current coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
C.2 DP model tuning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
C.3 Forbidden zones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
D Appendix 79
D.1 Tuning coefficients of 3DP model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
D.2 Motion time traces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
E Appendix 81
E.1 Yearly wave scatter Borssele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
E.2 Operability tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
F Appendix 85
F.1 Empirical bilge keel damping model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
F.2 Bilge keel model implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
F.3 Bilge keel model results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
List of Figures 89
List of Tables 91
Bibliography 93
1
Introduction
1.1. Background
The results of the UN Paris climate agreement in 2015 and the globally increasing environmental aware-
ness among people, made the need for clean and sustainable energy bigger then ever before. As a
consequence, the offshore wind industry has been rapidly expanding. The annual European offshore
wind turbine installations rose from 318.4 MW to 3018 MW between the years 2007 and 2015 [26].
Offshore contractors are hired to install these wind turbines. They are involved in the offshore con-
struction of, for example, foundations, monopiles, subsea power cables and transformer platforms. To
carry out these construction activities, a range of ships is used. These ships are all subjected to the
wind, waves and current present in the offshore environment. To enable safe and economical ship
operation, it is important to be able to predict the ship motion response to the environmental forces
present. Even more important would be the ability to actively control the ship motion response during
offshore operations. This enables the ship to operate more safely and economically.
Active ship motion control started in 1961 with the world’s first automatically controlled dynamically
positioned (DP) drillship Eureka, operated by Shell [20]. The Eureka was equipped with two thrusters
to control its position in the horizontal plane. It was a very successful ship. The Eureka drilled 9 cores
a day up to 150 meters in the sea floor at water depths up to 1200 meter, where comparably sized
anchored ships were able to do this in 2 to 4 days at a limited water depth of 100 meter.
Since the introduction of the Eureka, dynamically positioned ships became the standard for offshore
construction ships. Not only in the oil and gas sector, but also in the offshore renewable industry. Nowa-
days, numerous FPSOs, cable-laying, pipe-laying, heavy-lift and offshore supply ships are equipped
with a DP system to actively control the horizontal motions.
DP systems are currently used to control motions in the horizontal plane only. However, the safety
and operability of several offshore operations, like helicopter landings and lifting operations, can be in-
creased when also motions in the vertical plane can be controlled. A recent study conducted by Rudaa
[28] at NTNU showed that it is also possible to use the thrusters of the DP system to reduce the roll
motion of a ship.
Figure 1.1: Schematic visualization of the Ndurance
1
2 1. Introduction
1.2. Problem statement

The Ndurance is a DP offshore construction vessel operated by Boskalis. The vessel is currently used
for laying, connecting and repairing subsea power cables, rock bag installation, survey and mass flow
excavation operations. A visualization of the Ndurance is given in Figure 1.1. The vessel experiences
significant roll motions during beam waves in the natural frequency range. Several offshore DP op-
erations require a fixed heading and therefore the vessel is not always able to avoid beam waves by
weathervaning. The roll response during these circumstances can be so significant that it prevents
the crew from safely operating the equipment on deck, forcing the ship to wait for more advantageous
weather conditions. This phenomena is reducing the operability of the vessel.
To increase the operability of the vessel, the roll motion needs to be reduced. This can possibly be
achieved by using the thrusters of the installed DP system. Since the COG of the vessel is located
above the thrusters, the produced thrust results in a roll moment which can be used to counteract the
wave induced roll moment. Thruster controlled roll reduction is possible as concluded from mathemat-
ical and experimental model tests conducted at NTNU. However, the effect of active roll reduction on
the station keeping performance of the ship has not been investigated, since the available thrust is fully
used for roll reduction purposes. Also, the application of thruster induced roll reduction has only been
evaluated for an offshore supply vessel, which hull shape and DP system significantly differs from the
cable-lay vessel Ndurance.
1.3. Research objective

The main objective is to evaluate whether it is possible to reduce the significant roll amplitude of the
Ndurance in beam seas during DP operations by active control of the thrusters. To be able to achieve
the main objective, a merged mathematical ship, thruster and control system model with 3-dimensional
motion control capability has to be constructed (3DP). With this model it is also possible to evaluate
the effect of active roll reduction onto the station keeping performance of the ship. Recommendations
regarding possible thruster alternatives, ship lay-out and DP configuration that influence the possibility
and magnitude of actively controlled roll reduction are given. The sub-research objectives are formu-
lated as:
1. Develop a 6 DOF ship model in the time domain that captures the ship response to environmental
forces and validate this model with frequency domain data.
2. Develop a thruster model that captures the thruster dynamics in the time domain, this will include
modeling of the thruster drive train system and thrust variations due to oscillating inflow velocities.
3. Develop a simulation model that is able to simulate DP system behavior, this includes modelling
of the DP system and environmental forces.
4. Develop a control system model that enables combined DP station keeping and active roll reduc-
tion.
5. Identify and quantify the possibility and magnitude of thruster induced roll reduction and its effect
on the station keeping performance of the ship.
6. Give recommendations regarding possible thruster alternatives, the ship lay-out and DP system
configuration that influence the possibility and magnitude of actively controlled roll reduction.
The report is structured in such a way that the results of every sub-research objective will be presented
in the conclusions of the dedicated chapter.
1.4. Research approach 3
1.4. Research approach

First, the ship is mathematically modeled in the time domain in 6 DOF, according to the damping and
added mass as calculated by diffraction analysis. The time domain model results are compared to fre-
quency domain model results to verify that the implementation is correct. This is done by comparison
of the motion amplitude results when both models are subjected to regular harmonic waves. Subse-
quently, a viscous roll damping model is implemented in both the frequency and time domain model to
include viscous effects not captured by diffraction analysis.
Secondly, the thruster behavior in the time domain is modeled. Since the thrusters of the Ndurance are
fixed pitch propellers, the produced thrust is directly linked to the RPM of the electrical motor driving
the thrusters. To model the thruster response to control inputs and to capture the ramp-up and ramp-
down of the thrusters, a dynamic thruster model is constructed which includes the total thruster system
inertia. Next, the influence of oscillating inflow velocities, due to the roll velocity of the ship, on the
produced thrust and torque are captured by using the four-quadrant method.
Thirdly, the mathematical models of the thrusters and the ship are merged and the environmental
forces induced by wind, waves and current are generated in the time domain. Also, a control model is
implemented for station keeping (DP). This model enables the simulation of the 6 DOF ship motions
during DP operations. The model is verified with MARIN data.
Fourthly, a new control model system is proposed that combines anti-roll control and station keep-
ing control. This model is used to evaluate the possible roll reduction and the impact of the roll control
mode onto the DP footprint, station keeping capability, power consumption and workability.
Finally, the results of the simulations are analyzed by comparing the root mean square values of the
roll amplitudes with and without active roll control. Also, the DP footprint, DP capability, thruster power
consumption and vessel workability is evaluated for different environmental conditions. According the
results, recommendations will be given regarding the influence of different ship parameters on the
possibility and magnitude of thruster induced roll reduction.
1.5. Report structure

The report structure is shortly summarized in the enumeration below.
• In Chapter 2 (Ship Model), models to calculate a ship’s motion response are presented and
reviewed. A time domain ship model is constructed and verified.
• In Chapter 3 (Thruster Model), a dynamic thruster model is presented and implemented. Also,
the thruster behavior in oscillating flow is included in the dynamic model. The results are verified
with manufacturer data and expert indication
• In Chapter 4 (DP System Model), the DP control system simulation model is developed. The
results are verified with time domain test data.
• In Chapter 5 (3DP System Model), the anti-roll controller is introduced and merged with the DP
system to obtain a new simulation model for roll reduction and station keeping combined. This
system is called 3DP.
• In Chapter 6 (Results), the results of the 3DP system model regarding the achieved roll reduction
are evaluated. Also, the DP footprint, DP capability, thruster power consumption and workability
of the 3DP model is evaluated and compared to the performance of the conventional DP system
model.
• In Chapter 7 (Conclusions and Recommendations), conclusions and recommendations are
given regarding the model results and the conducted research.
2
Ship Model
In order to evaluate the possibility and magnitude of thruster induced roll reduction and its effect on the
DP footprint of the Ndurance, it is important to construct a model that is able to calculate the 6 DOF
ship motion response, resulting from environmental and thruster loads. The method and models to
estimate and simulate ship motions are presented in this chapter.
2.1. Introduction
In literature there exists two different models to capture the physical behavior of a ship in a mathematical
hydromechanic model: the frequency domain model and the time domain model. Both models and their
underlaying assumptions are discussed. Subsequently, since no viscous effects are included in either
models, models to capture these viscous effects are presented and discussed.
2.1.1. Ship motions

A ship in a seaway has the ability to move in 6 directions, so called degrees of freedom. These degrees
of freedom, or shortly DOFs, are visualized in Figure 2.1 and described in Table 2.1. The origin of the
coordinate system is located at the ship’s center of gravity (COG).
Motion Name Symbol Unit

1 Surge x [m]
2 Sway y [m]
3 Heave z [m]
4 Roll 𝜙 [rad]
5 Pitch 𝜃 [rad]
6 Yaw 𝜓 [rad]
Figure 2.1: Ship motions in 6 degrees of freedom [20] Table 2.1: Description of ship motions and symbols
2.2. Frequency domain model

The frequency domain model is the most widely used method to calculate the response of a ship in
a seaway. The approach is based on the assumption that the ship is excited by purely sinusoidal
wave forces at a particular well defined frequency and that the resulting motions are small relative to
the cross-sectional dimensions of the ship. Only the steady-state response can be be obtained using
the frequency domain method. The model to calculate the ship response at zero-speed is equivalent
to a frequency-dependent mass-damper-spring system. The model is mathematically defined by an
’apparent’ differential equation according Cummins [7]:
[M + A(𝜔)]𝜂̈ + B(𝜔)𝜂̇ + C𝜂 = F(𝜔) cos(𝜔𝑡 + 𝜖) (2.1)
5
6 2. Ship Model
in which:
M = ship inertia matrix

A(𝜔) = ship hydrodynamic added mass matrix
B(𝜔) = ship hydrodynamic damping matrix
C = ship hydrostatic stiffness matrix
F(𝜔) = first-order wave excitation force amplitude vector
𝜔 = wave frequency
𝜂 = ship motion amplitude vector, see Table 2.1
𝜖 = phase angle vector
𝑡 = time vector
The ship is considered laterally symmetric. Therefore, no coupling terms other than the surge-heave-
pitch and sway-roll-yaw coupling terms are included in the model. The frequency independent inertia
matrix and hydrostatic stiffness matrix become:
𝑚 0 0 0 0 0 0 0 0 0 0 0
⎡ ⎤ ⎡ ⎤
0 𝑚 0 0 0 0 0 0 0 0 0 0
⎢ ⎥ ⎢ ⎥
0 0 𝑚 0 0 0 0 0 𝑐 0 𝑐 0⎥
M=⎢ ⎥ C=⎢
⎢ 0 0 0 𝑚 0 0 ⎥ ⎢0 0 0 𝑐 0 0⎥
⎢ 0 0 0 0 𝑚 0 ⎥ ⎢0 0 𝑐 0 𝑐 0⎥
⎣ 0 0 0 0 0 𝑚 ⎦ ⎣0 0 0 0 0 0⎦
Where the first three diagonal terms in the inertia matrix represent masses and the last three terms
represent rotational moments of inertia.
The energy that is carried away from the system due to the generation of surface waves as a result of
the ship motions is captured in the radiation force, existing out of the added mass A(𝜔) and potential
damping B(𝜔) according Fossen [11]:
𝑎 0 𝑎 0 𝑎 0 𝑏 0 𝑏 0 𝑏 0
⎡ ⎤ ⎡ ⎤
0 𝑎 0 𝑎 0 𝑎 0 𝑏 0 𝑏 0 𝑏
⎢ ⎥ ⎢ ⎥
𝑎 0 𝑎 0 𝑎 0 𝑏 0 𝑏 0 𝑏 0
A(𝜔) = ⎢ ⎥ B(𝜔) = ⎢ ⎥
⎢ 0 𝑎 0 𝑎 0 𝑎 ⎥ ⎢ 0 𝑏 0 𝑏 0 𝑏 ⎥
⎢𝑎 0 𝑎 0 𝑎 0 ⎥ ⎢𝑏 0 𝑏 0 𝑏 0 ⎥
⎣ 0 𝑎 0 𝑎 0 𝑎 ⎦ ⎣ 0 𝑏 0 𝑏 0 𝑏 ⎦
Due to lateral symmetry, the coupling terms in the above matrices are equal, so 𝑐 = 𝑐 , this also
applies to 𝑎 and 𝑏.
The first order wave excitation force amplitude is calculated by integration of the wave pressure acting
on the submerged ship hull in a undisturbed wave. This is commonly referred to as the Froude-Krilov
force according Journée [20]. The Froude-Krilov force is corrected since a part of the waves will be
diffracted due to the presence of the hull in the fluid. The corrected first-order wave excitation force
amplitude is denoted by F(𝜔):
𝐹
⎡ ⎤
𝐹
⎢ ⎥
𝐹
F(𝜔) = ⎢ ⎥
⎢𝐹 ⎥
⎢𝐹 ⎥
⎣𝐹 ⎦
Note that the last three terms in the first order wave excitation force amplitude vector represent mo-
ments, since the corresponding motions are rotations.
2.2. Frequency domain model 7
The added mass coefficients, damping coefficients, Froude-Krilov and diffraction forces can be ob-
tained by using diffraction analysis software, such as for example, AQWA or WAMIT. In diffraction
analysis, the ship geometry is represented by a discretized panel model. The theory underlying these
solvers is very clearly and deliberately described by Journée [20] and Newman [21]. The discretized
panel model of the Ndurance as implemented in AQWA is visualized in Figure 2.2. The vessel partic-
ulars are given in Appendix A.
Figure 2.2: Panel model of the Ndurance in AQWA
Diffraction analysis is based on the assumption of potential flow, implying that the flow is assumed to
be and to stay irrotational and non-viscous. As a consequence, it is not possible to calculate viscous
effects by using diffraction analysis software.
Since a linear model is used to calculate the ship motions, the response of the ship in irregular waves
can be calculated by using the superposition principle. The superposition principle allows the calcu-
lation of a motion response spectrum by using linear Response Amplitude Operators (RAO) and a
defined wave spectrum according Journée [20]:
𝜂(𝜔)
𝑆 (𝜔) = | | ⋅ 𝑆 (𝜔) (2.2)
𝜉
in which:
𝑆 (𝜔) = motion response spectrum

𝜂(𝜔)
= motion RAO (motion amplitude / wave amplitude)
𝜉
𝑆 (𝜔) = wave spectrum
The RAO depends on the water depth, wave direction, draft and loading condition of the ship. A visu-
alization of the ship motion RAOs are given in Figure 2.3.
Several statistical parameters of the ship motion response can be calculated according to the mo-
tion response spectrum. For example, the root mean square (RMS) value of the motion amplitudes, is
calculated according Journée [20] by:
𝜂 = √𝑚 = √∫ 𝑆 (𝜔)𝑑𝜔 (2.3)
The root mean square value of the motion amplitude can be used to assess the operability of a specific
offshore operation for each specific environmental condition characterized by a wave spectrum.
8 2. Ship Model
Surge RAO Sway RAO

1 1
Surge amplitude [m/m]
Sway amplitude [m/m]

0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Frequency [rad/s] Frequency [rad/s]
Heave RAO Roll RAO

1.2 10
Heave amplitude [m/m]
Roll amplitude [deg/m]

1 8
0.8
6
0.6
4
0.4
0.2 2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Frequency [rad/s]
Pitch RAO Yaw RAO
2 1
Pitch amplitude [deg/m]
Yaw amplitude [deg/m]
0.8
1.5
0.6
1
0.4
0.5
0.2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Figure 2.3: Ndurance Motion RAOs ( ∘, T=4.7m, WD=1000m)
2.3. Time domain model

A limitation of the frequency domain model is its disability to calculate the ship response when it is sub-
jected to non-linear or transient loads. To overcome this limitation, the time domain model as proposed
by Cummins [7] can be used to calculate the transient ship response at zero speed by an integro-
differential equation:
[M + A ]𝜂̈ + ∫ K(𝑡 − 𝜏)𝜂𝑑𝜏

̇ + C𝜂 = 𝐹(𝑡) (2.4)
in which:
A = infinite frequency added mass matrix

K = retardation function matrix
𝐹(𝑡) = external excitation force
In the time domain model, the hydrodynamic damping term is represented by the convolution of the re-
tardation function and the ship velocities. The infinite frequency added mass matrix and the retardation
2.3. Time domain model 9
function matrix are defined by:
𝑎 0 𝑎 0 𝑎 0 𝐾 0 𝐾 0 𝐾 0
⎡ ⎤ ⎡ ⎤
0 𝑎 0 𝑎 0 𝑎 0 𝐾 0 𝐾 0 𝐾
⎢ ⎥ ⎢ ⎥
𝑎 0 𝑎 0 𝑎 0 𝐾 0 𝐾 0 𝐾 0
A =⎢ ⎥ K=⎢ ⎥
⎢ 0 𝑎 0 𝑎 0 𝑎 ⎥ ⎢ 0 𝐾 0 𝐾 0 𝐾 ⎥
⎢𝑎 0 𝑎 0 𝑎 0 ⎥ ⎢𝐾 0 𝐾 0 𝐾 0 ⎥
⎣ 0 𝑎 0 𝑎 0 𝑎 ⎦ ⎣ 0 𝐾 0 𝐾 0 𝐾 ⎦
The retardation function is the impulse response function of the ship. It takes into account the so-
called fluid memory. The retardation function can be calculated for every ship motion and coupling
term according Fossen [11] by:
2
K(𝑡) = ∫ B(𝜔) cos(𝜔𝑡)𝑑𝜔 (2.5)
𝜋
The infinite frequency added mass is the hydrodynamic added mass evaluated at infinite frequency. It
is defined according Fossen [11] by:
1
A = A = A(𝜔) + ∫ K(𝑡) sin(𝜔𝑡)𝑑𝑡 (2.6)
𝜔
As can be derived from equation 2.6, the infinite frequency added mass is calculated by averaging over
the complete frequency range. It can be derived from the equations of the retardation function and the
infinite frequency added mass, that the time domain approach is constructed from hydrodynamic coef-
ficients in the frequency domain. Therefore, the added mass and damping coefficients, calculated by
diffraction analysis can be used to develop a time domain model.
Because of the ability to calculate the vessel transient response when subjected to non-linear forces
and for simulation and controller design purposes, see Fossen [11], it is decided to simulate the motion
response of the Ndurance in the time domain.
2.3.1. The convolution term

The convolution integral in equation 2.4, can make the simulation of the ship motion time and memory
consuming according Armesto [1]. This is the result of the fact that the convolution integral equation
requires to integrate all the way back to the start of the simulation at every time step. This causes
simulations to run increasingly slower when the simulation time increases. In literature, three different
approaches to compute the radiation term in Cummins’ equation can be found:
• Direct integration of the convolution term by the Impulse Response Function (IRF), see Ricci [27].
• Computation of the convolution term by a State-Space (SS) model, see Perez [25].
• Computation of the convolution term by using Prony’s coefficients to estimate the IRF, see De
Backer [8].
In literature, the state-space method is commonly used to model the dynamic ship response in the
time domain. This is due to the computational performance of the method and its suitability for control
purposes, see Hatecke [13]. Due to these advantageous properties it is decided to use the state-space
method to evaluate the convolution term in Cummins’ equation.
State-space representation requires that a dynamic system is linear and time-invariant. Since the radi-
ation term in Cummins’ equation satisfies these conditions, it is possible to compute the fluid-memory
effects by an approximated state-space model:
𝑥̇ = A 𝑥 + B 𝜂̇
𝜅 = ∫ K(𝑡 − 𝜏)𝜂𝑑𝜏
̇ ≃| (2.7)
𝜅= C 𝑥
10 2. Ship Model
in which:
𝑥̇ = state vector
A = state-space system matrix
B = state-space input matrix
C = state-space output matrix
The strategy to construct the state-space models is based on the work of Perez and Fossen [24]. First,
the retardation functions, or IRFs, are calculated in the time domain in order to calculate the values
of the infinite frequency added mass matrix. According these results, the IRF is transformed to the
frequency domain. This feels like a natural approach, since the added mass and damping coefficients
as calculated by diffraction analysis are already represented in the frequency domain. Next, a transfer
function (TF) is fitted through the IRF and the fitted transfer function is subsequently transformed back
to the time domain by the state-space representation. A diagram of the strategy to construct the state-
space models is given in Figure 2.4.
Time domain IRF(t) SS(t)
Frequency domain IRF(𝜔) TF(𝜔)
Figure 2.4: Strategy to obtain state-space representation
It may seem a more logical approach to directly transform the IRF to the state-space model in the
time domain. However, this is significantly more complex, due to the fact that the order of the system
is hard to guess by looking at the IRF alone and it is not possible to enforce the model to dynamic
system constraints. The calculation steps visualized in Figure 2.4 will be presented in more detail in
the following sections.
2.3.2. IRF in the time domain

The IRF needs to be calculated in the time domain to obtain an expression for the infinite added mass
terms, see Equation 2.6. The formula for the IRF is evaluated at infinite frequency. This yields prob-
lems, since the IRF needs to be calculated by a numerical model. Therefore, the function is only
evaluated until a cut-off frequency, in this case 4.5 rad/s. The cut-off frequency is determined by the
diffraction software, since computational-wise, the software is not able to calculate higher frequency
data for the panel model used.
Two different numerical models to calculate the IRF of the ship motion and coupling terms are evalu-
ated in this study. One is used in the simulation program OrcaFlex [23] and the other is proposed by
Journée [16].
The model used in OrcaFlex [23] is formulated by:
−3𝑡
K(𝑡) = exp [( ) ] ⋅ ∑ 2Δ𝐵 (𝜔)Δ𝜔 cos(𝜔𝑡) (2.8)
𝑇
in which:
𝑁 = total number of frequencies
𝑡 = time vector
Δ𝐵 = hydrodynamic damping, 𝐵 − 𝐵
Δ𝜔 = wave frequency, 𝜔 − 𝜔
𝑇 = cut-off scaling time constant
The variable 𝑇 can be adjusted to smoothly scale down the IRF function. A high value results in low
smoothing of the IRF.
According to Journée [16], the discretized model of the IRF is:
2 Δ𝐵 2
K(𝑡) = ⋅∑( [ cos(𝜔 𝑡) − cos(𝜔 𝑡)]) + ⋅ 𝐵 sin(𝜔 𝑡) (2.9)
𝜋𝑡 Δ𝜔 𝜋𝑡
in which:
𝐵 = hydrodynamic damping at the cut-off frequency
𝜔 = cut-off frequency
Both models have been implemented to check whether they are in agreement and to confirm the correct
implementation of the numerical IRF model. The calculated IRFs for the six motion components using
both methods are given in Figure 2.5.
105 IRF Surge 106 IRF Sway

6 6
OrcaFlex OrcaFlex
4 Journee 4 Journee
2 2
K11
K22
0 0
-2 -2
-4 -4
0 5 10 15 20 25 0 5 10 15 20 25
Time [s] Time [s]
106 IRF Heave 107 IRF Roll

10 8
OrcaFlex OrcaFlex
6
Journee Journee
5 4
K33
K44
0 0
-2
-5 -4
0 5 10 15 20 25 0 5 10 15 20 25
Time [s] Time [s]
109 IRF Pitch 109 IRF Yaw

6 2
OrcaFlex OrcaFlex
4 Journee Journee
1
2
K55
K66
0
0
-1
-2
-4 -2
0 5 10 15 20 25 0 5 10 15 20 25
Time [s] Time [s]
Figure 2.5: Calculated motion IRFs using the OrcaFlex and Journée models, =50s
From Figure 2.5 can be concluded that both models result in the same estimation of the motion IRFs.
The effect of the smoothing function defined by the variable 𝑇 is clearly visible in the results. The
OrcaFlex model estimates smoothen out significantly quicker compared to the Journée model. To
increase the accuracy of the infinite frequency added mass calculation, it is decided to use the IRF
model as proposed by Journée [16].
12 2. Ship Model
2.3.3. IRF in the frequency domain

The calculated IRF in the time domain needs to be subsequently transferred to the frequency domain.
This is achieved according Fossen [11] by:
K(𝜔) = B(𝜔) + 𝑖𝜔[A(𝜔) − A ] (2.10)
From Equation 2.10 can be concluded that first the infinite frequency added mass matrix needs to be
calculated. This is done according the calculated IRFs for the ship motion and coupling terms, see
Equation 2.6. Again the integration needs to be cut-off, since integration at infinity is not possible using
numerical techniques. The results of the infinite frequency added mass terms for the six motion com-
ponents are calculated and plotted together with the frequency-dependent added mass estimates as
calculated by the AQWA diffraction software in Figure 2.6.
105 Surge added mass 106 Sway added mass
10 AQWA data AQWA data

A 6 A
Added mass [kg]
Added mass [kg]

8
6 4
4
2
2
0 0
0 1 2 3 4 5 0 1 2 3 4 5
107 Heave added mass 108 Roll added mass

7 10
AQWA data AQWA data
6
Added mass [kgm 2 ]
A 9.5 A
Added mass [kg]
5
9
4
8.5
3
8
2
7.5
0 1 2 3 4 5 0 1 2 3 4 5
1010 Pitch added mass 109 Yaw added mass

3.5
AQWA data 3 AQWA data
2.5
Added mass [kgm 2 ]
Added mass [kgm 2 ]
A A
2.5
2 2
1.5
1.5 1
0.5
1 0
0 1 2 3 4 5 0 1 2 3 4 5
Figure 2.6: Frequency-dependent added mass for six ship motions and calculated infinity frequency added mass estimates
From Figure 2.6 can be observed that the added mass estimates as calculated by the diffraction soft-
ware AQWA approach the calculated infinite frequency added mass estimate in the high frequency
limit. From this can be concluded that the numerical model to calculate the infinity frequency added
mass is implemented correctly and produces reliable results.
2.3.4. Transfer function fitting

The next step is to fit transfer functions through the calculated IRFs in the frequency domain. An
extensive paper regarding the fitting method and procedure is written by Duarte [10]. To obtain the
coefficients matrices A , B and C of the state-space models, a transfer function fitting algorithm
is used. The steps executed in the algorithm are described in the following enumeration:
1. A parametric transfer function is fitted through the retardation function in the frequency domain
defined by:
𝑃(𝑠) 𝑝 𝑠 +𝑝 𝑠 + ... + 𝑝
̃
𝐾(𝑠) = = (2.11)
𝑄(𝑠) 𝑠 +𝑞 𝑠 + ... + 𝑞
in which:
𝑠 = Laplace complex variable

𝑝 = nominator coefficients
𝑞 = denominator coefficients
𝑚 = order of nominator
𝑛 = order of denominator
2. A least squares method is solved to minimize the error between the fitted transfer function and
the retardation function:
min ∑ 𝑤 (𝑘)|𝐾(𝑘)𝑄(𝑤(𝑘)) − 𝑃(𝑤(𝑘))| (2.12)
In which 𝑤 is a vector with weight factors and 𝑛 is the number of frequency points. This method
is implemented in the Matlab function ’invfreqs’.
3. The weight factors are calculated by:
1
𝑤 (𝑘) = (2.13)
|𝑝(𝑘) + 𝑖𝜔(𝑘)|
4. The quality of the fit is calculated by:

̃
∑(𝐾(𝑠) − 𝐾(𝑠))
𝑅 =1− (2.14)
∑(𝐾(𝑠) − 𝐾(𝑠))
According to Perez [25], some specific properties of the convolution terms have implications on the
parametric transfer function models. The convolution term properties and their implications are listed
in Table 2.2.
Property Implication on parametric model

lim K(i𝜔)=-B(∞) K(s) is zero at s=0
→
lim K(i𝜔)=0 TF is strictly proper
→
K(t=0 )=∫ [𝐵(𝜔) − 𝐵(∞)]𝑑𝜔 ≠ 0 TF relative degree is 1
lim K(t)=0 TF is BIBO stable
→
Re(K(i𝜔))≥0 K(s) is passive, diagonal terms are positive real
Table 2.2: Convolution term properties and corresponding implications on the fitted transfer function (TF) according Perez [25]
The properties as given in Table 2.2 are incorporated in the transfer function fitting algorithm. The
minimum order of the transfer function is set to 2, since it is not possible to fit a first-order transfer
function to the retardation function. In the fitting procedure this order is increased and the weight
vector is re-evaluated until the transfer function estimate has a fit quality of at least 0.9. The results of
the fitting procedure are given in Figure 2.7.
14 2. Ship Model
105 IRF Surge 106 IRF Sway

6 5
K( ) K( )
5 4
TF fit TF fit
4
3
K11
K22
3
2
2
1 1
0 0
10-2 10-1 100 101 10-2 10-1 100 101
107 IRF Heave 107 IRF Roll

2 8
K( ) K( )
TF fit TF fit
1.5 6
K33
K44
1 4
0.5 2
0 0
10-2 10-1 100 101 10-2 10-1 100 101
109 IRF Pitch 109 IRF Yaw
8 2.5
K( ) K( )
TF fit 2 TF fit
6
1.5
K55
K66
4
1
2
0.5
0 0
10-2 10-1 100 101 10-2 10-1 100 101
Figure 2.7: Frequency domain motion IRFs and fitted transfer functions
A transfer function is fitted through every motion and coupling term component. As can be derived from
Figure 2.7, the fitted transfer functions are approximating the IRF functions quite well. The roll and yaw
motions show some slight underestimation of the peaks at the higher frequencies, but it is expected
that these effects are negligible. The transfer functions are subsequently transformed to a state-space
model using the ’tf2ss’ function in Matlab.
2.3.5. Time domain model verification

The time domain model in which the convolution terms are replaced by equivalent state-space repre-
sentations is defined by:
[M + A ]𝜂̈ + K 𝜂̇ + C𝜂 = 𝐹(𝑡) (2.15)
where,
𝜅 0 𝜅 0 𝜅 0
⎡ ⎤
0 𝜅 0 𝜅 0 𝜅
⎢ ⎥
𝜅 0 𝜅 0 𝜅 0
K =⎢ ⎥
⎢ 0 𝜅 0 𝜅 0 𝜅 ⎥
⎢𝜅 0 𝜅 0 𝜅 0 ⎥
⎣ 0 𝜅 0 𝜅 0 𝜅 ⎦
The 𝜅-terms represent the state-space representations of the motion and coupling terms defined by
the state-space coefficients matrices A , B and C .
The time domain model has been constructed in the simulation software Simulink. Both the frequency
and time domain model are subjected to equal harmonic quartering wave forces and moments over a
frequency interval of 0.1 to 2.5 rad/s. In this way the RAOs as visualized in Figure 2.3 are reconstructed
and the time domain model is checked for the complete range of frequencies. The motion amplitude
results as calculated by both models are compared to check whether the results are in agreement. The
results are visualized in Figure 2.8.
Surge RAO Sway RAO

1 1
Frequency domain model Frequency domain model
Surge amplitude [m/m]
Sway amplitude [m/m]
0.8 Time domain model 0.8 Time domain model
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Heave RAO Roll RAO

1.2 10
Heave amplitude [m/m]
Roll amplitude [deg/m]
1 8
Time domain model Time domain model
0.8
6
0.6
4
0.4
0.2 2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Frequency [rad/s]
Pitch RAO Yaw RAO
2 1
Pitch amplitude [deg/m]
Yaw amplitude [deg/m]
Time domain model 0.8 Time domain model

1.5
0.6
1
0.4
0.5
0.2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Figure 2.8: Comparison of calculated RAOs by both frequency and time domain model ( ∘, T=4.7m, WD=1000m)
16 2. Ship Model
From Figure 2.8 can be observed that the calculated motion amplitudes are in good agreement for
all 6 motion components. From these results is concluded that the constructed time domain model is
implemented correctly, since it produces motion amplitude estimates similar to the frequency domain
model when subjected to harmonic forces and moments.
To quantify the agreement between the motion estimates obtained by both models, the Relative Mean
Absolute Error (RMAE) of each motion is calculated. The RMAE is defined by:
∑ |𝑦 − 𝑦̂ |
RMAE = ⋅ 100% (2.16)
|𝑦 −𝑦 |
in which
𝑛 = total number of data points

𝑦 = frequency domain result
𝑦̂ = time domain result
𝑦 = maximum value in result range
𝑦 = minimum value in result range
The RMAE value for each motion component visualized in Figure 2.8 is calculated and are given in
Table 2.3.
Motion RMAE [%]

Surge 0.36
Sway 0.25
Heave 0.59
Roll 0.05
Pitch 0.59
Yaw 0.33
Table 2.3: RMAE values of indicated ship motions
From Table 2.3 can be derived that that the maximum RMAE value is 0.59% for both the heave and
pitch motion. A magnitude of 0.59% is considered negligible and therefore it is concluded that the
time domain model with state-space representation is capable to correctly calculate the vessel motion
response.
2.4. Viscous roll damping 17
2.4. Viscous roll damping

Viscous effects are neglected when potential theory is used to calculate ship motions. This means that
the damping only represents the energy dissipated by the radiated waves. Since the motions of a ship
are largely inertia driven, the results obtained by potential theory are in general sufficiently accurate to
correctly calculate the motion response.
For the roll motion however, viscous effects are of more importance. A significant amount of energy is
dissipated by the formation of eddies at the bilges. This energy contributes to the roll damping of the
ship. The contribution of the viscous damping to the total roll damping of the ship is found to be difficult
to determine. Two approaches are commonly seen in literature:
1. Experimental model tests
2. Empirical models
Since no model test results of the Ndurance are available, an empirical model approach is used to take
into account the effect of viscous roll damping.
2.4.1. Ikeda viscous damping model

As stated by Chakrabarti [5], the viscous roll damping term is a function of the roll velocity and can be
expressed by:
̇ = 𝐵 𝜙̇ + 𝐵 |𝜙|̇ 𝜙̇ + 𝐵 𝜙
𝐵(𝜙) (2.17)
where 𝐵 , 𝐵 and 𝐵 represent linear, quadratic and cubic viscous damping coefficients. Since no higher
order terms can be incorporated in the frequency domain model, the total damping is subsequently
approximated by a linear equivalent damping term:
̇ = 𝐵 𝜙̇
𝐵(𝜙) (2.18)
According to Himeno [14], the equivalent viscous damping coefficient can be expressed in terms of its
various contributions:
𝐵 =𝐵 +𝐵 +𝐵 +𝐵 (2.19)
in which:
𝐵 = skin friction damping

𝐵 = eddy making damping
𝐵 = lift damping
𝐵 = bilge keel damping
For each of this contributions a specific empirical model has been constructed according to experimen-
tal results. Since the Ndurance has no bilge keels and only the zero-speed condition is of interest, this
equation can be reduced to:
𝐵 =𝐵 +𝐵 (2.20)
A schematic visualization of the eddy making and skin friction damping is given in Figure 2.9.
z
Figure 2.9: Schematic visualization of eddy making damping and skin friction damping
18 2. Ship Model
Eddy making damping

Eddy making damping arises due to the separation of vortices at the bilges of the ship. For ships with
a barge-like rectangular midsection like the Ndurance, the contribution of the eddy making damping is
most significant amidships. According to Ikeda [15], the eddy making damping coefficient of the bare
ship hull with a rectangular cross-section is formulated by:
2
𝐵 = 𝜌𝐿𝐷 (𝐻 + 1 − 𝑂𝐺/𝐷)[𝐻 + (1 − 𝑂𝐺/𝐷) ]𝑅 𝜔 (2.21)
𝜋
in which:
𝜌 = sea water density

𝐿 = ship length
𝐷 = draft
𝐻 = half of the beam-draft ratio (B/2D)
𝑅 = roll motion amplitude
𝑂𝐺 = vertical distance from still water line (SWL) to COG
The eddy making damping model assumes that the hull has a rectangular (90∘ ) bilge and that the rect-
angular cross sections is present over the complete length of the ship (section coefficient 𝐶 = 1).
The Ndurance has more ship-like sections in the aft and front of the ship. Next to this, the Ndurance
has a slightly rounded bilge. Therefore, the length of the part of the ship at which the sections have a
𝐶 ≥ 0.99 has been selected as an equivalent length to be used in Equation 2.21. The length for which
𝐶 ≥ 0.99 is 54 meters.
Skin friction damping

Friction damping is the result of skin friction stress that acts on the ship hull. The skin friction stress
is induced by the sea water viscosity and the roll motion of the ship. According Himeno [14], the skin
friction damping coefficient can be expressed by:
.
𝑟 𝑅 𝜔
𝐵 = 0.787𝜌𝑆𝑟 √𝜔𝜈{1 + 0.00814( ) } (2.22)
𝜈
in which:

𝑆 = wetted hull surface
𝑟 = effective bilge radius
𝜈 = kinematic viscosity of water
The effective bilge radius 𝑟 is expressed by:
1 𝑆
𝑟 = [(0.887 + 0.145𝐶 ) − 2𝑂𝐺] (2.23)
𝜋 𝐿
in which:
𝐶 = ship block coefficient

2.4. Viscous roll damping 19
2.4.2. Frequency domain model

As can be derived from Equation 2.21 and Equation 2.22, both the eddy making damping coefficient
and the skin friction damping coefficient depend on the wave frequency and the roll motion amplitude.
However, the roll motion amplitude also depends on the damping coefficients. To include this interac-
tion, the damping coefficients are calculated by using an iterative approach. The damping coefficients
are calculated by using the roll motion amplitude estimates of the frequency domain model without
the implementation of a viscous damping term. Subsequently, the viscous damping coefficients are
calculated and implemented in the frequency domain model:
[M + A(𝜔)]𝜂̈ + B(𝜔)𝜂̇ + 𝐵 𝜙̇ + C𝜂 = F(𝜔) cos(𝜔𝑡 + 𝜖) (2.24)
The calculated roll motion amplitudes are then used to calculate a new viscous damping coefficient
estimate for every frequency. This process is repeated until the viscous damping term is converged.
The results of the iterative procedure are visualized in Figure 2.10.
108 Eddy making damping 106 Skin friction damping

14 5
Viscous roll damping [kgm/s rad]

Converged damping RAO Converged damping RAO
12 Iterations Iterations
4
10
8 3
6 2
4
1
2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
108 Converged total viscous damping Converged roll RAO

2.5 10
2 8
Roll RAO [deg/m]
1.5 6
1 4
0.5 2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Figure 2.10: Results of iterative calculation of viscous damping RAOs and final roll motion RAO in the frequency domain model
( ∘ , T=4.7m, WD=1000m)
From Figure 2.10 can be observed that both the damping RAOs and the roll motion RAO converge
to a final value. Also, the total converged viscous damping term RAO and the roll motion RAO are
visualized.
20 2. Ship Model
2.4.3. Time domain model

In the time domain model it is possible to incorporate a quadratic damping term. This is considered a
more physical approach, since it is known that drag forces have a quadratic relation with velocity, see
for example Chakrabarti [5]. Therefore a quadratic viscous damping term is included, defined by:
𝐵 = 𝑐 |𝜙|̇ 𝜙̇ (2.25)
where 𝑐 is a fitting parameter. The value of this fitting parameter is obtained by fitting the time domain
roll motion RAO to the frequency domain roll motion RAO with a linearized viscous roll damping term
included. The time domain model with a viscous roll damping term included is now defined by:
[M + A ]𝜂̈ + K 𝜂̇ + 𝑐 |𝜙|̇ 𝜙̇ + C𝜂 = 𝐹(𝑡) (2.26)
By adjusting the value of 𝑐 , the time domain roll motion RAO is fitted to the frequency domain roll
motion RAO. The result is visualized in Figure 2.11
9
Frequency domain model
8 Time domain model
6
Roll RAO [deg/m]
0
0 0.5 1 1.5 2 2.5
Frequency [rad/s]
Figure 2.11: Results of time domain model RAO fitting to frequency domain RAO ( ∘, T=4.7m, WD=1000m)
From Figure 2.11 can be observed that the roll motion RAO as calculated by both models are in good
agreement. The value of the tuning parameter is determined as 𝑐 = 2.2⋅10 . To quantify the agreement
between the frequency domain model and the time domain model, the RMAE is calculated according
Equation 2.16. From this calculation is found that the RMAE of the results is 0.11%. This is considered
negligibly small and therefore it is concluded that the match between the time domain model and the
frequency domain model is excellent.
One should be aware of the fact that the frequency and time domain roll motion RAO only match for
harmonic waves with unit amplitude. For higher wave amplitudes, the roll motion velocity will increase
and therefore the time domain model will overestimate the viscous roll damping, due to the quadratic
term, in comparison to the frequency domain model. For lower wave amplitudes, the opposite applies.
This behavior is visualized in Figure 2.12.
2.5
Viscous damping trend [-]

2 Time domain model
1.5
0.5
0
0 0.5 1 1.5
Wave amplitude a
[m]
Figure 2.12: Viscous roll damping trend visualization in both frequency and time domain model
Since it is known that the viscous damping term has a quadratic relation regarding the roll velocity, it is
assumed that the time domain model yields more physically correct results.
2.5. Conclusions 21
2.5. Conclusions
The first sub-research objective is formulated as:
”Develop a 6 DOF ship model in the time domain that captures the ship response to environmental
forces and validate this model with frequency domain data.”
Two models to calculate the ship motion response have been presented and evaluated in this chap-
ter. The time domain model is constructed based on state-space representation of the convolution
term. Empirical models are used to calculate the viscous damping terms in the frequency domain. The
viscous roll damping term is incorporated in the time domain by inclusion of a quadratic fitting coefficient.
From this research is concluded that:

• It is possible to evaluate the convolution term in the time domain model by an equivalent state-
space model.
• The maximum relative mean absolute error between the frequency domain model motion results
and the time domain model motion results is 0.59%. This is considered small enough to conclude
that the time domain model is verified by frequency domain results.
• It is possible to use a quadratic parameter in order to match the frequency domain and the time
domain model roll motion RAO with a viscous damping term included.
• The maximum relative mean absolute error between the frequency domain roll RAO and the tuned
time domain roll RAO is 0.11%. This error is considered negligible and therefore it is concluded
that time domain model tuning procedure is verified for unit wave amplitudes.
3
Thruster Model
To achieve thruster induced roll reduction, the thrusters need to counteract the first-order roll moment.
The dynamic response of the thrusters is therefore important to model. Methods and models to estimate
and simulate the thruster performance in transient condition are presented in this chapter.
3.1. Introduction
The Ndurance is equipped with a diesel-electric propulsion system. The thruster system lay-out struc-
ture is schematically visualized in Figure 3.1
Dynamic model
Generator VFD Motor Shaft Propeller
Figure 3.1: Schematic structure of thruster system lay-out
The generator is driven by a diesel engine that ensures the system is fed with electrical energy. The
variable frequency drive (VFD) is a controller that regulates the motor torque and speed by controlling
the stator magnetic motor flux. The required flux is calculated according the so-called vector control
method. See for example Bargmeyer [2]. The electric motor is a three-phase asynchronous motor.
The motor is connected to a drive shaft which transfers the rotational movement and torque of the
motor to the gearbox. The gearbox reduces the rotational speed and increases the torque, which are
subsequently transferred to the propeller. The rotation of the propeller results in thrust. Since the pro-
peller is a fixed pitch propeller, the system dynamics are characterized by the components indicated
in Figure 3.1. A model is presented which is able to capture the dynamic behavior of the thruster system.
The propeller will operate in fluctuating inflow conditions due to the roll motion of the vessel, see Figure
3.2.
z
+
- - +
Figure 3.2: Schematic visualization of oscillating thruster inflow velocities due to the rolling motion of the vessel
As a result of this, the thrust and torque production will also change. A method to capture the thruster
behavior when operating in these fluctuating inflow conditions will be presented and implemented.
23
24 3. Thruster Model
3.2. Dynamic model

During roll reduction mode, the thrusters are used to counteract the first-order roll moment. This en-
tails that the thrusters will be ramping up and down during one wave period. The response time of
the thrusters is considered an important parameter regarding the possibility and magnitude of thruster
induced roll reduction. To analyze the response time of the thrusters, a dynamic thruster model is con-
structed.
In literature there exist several modelling approaches. Cooke [6] proposes an energy-based dynamic
model based on lumped parameters:
𝜔̇ = 𝛽𝑄 − 𝛼𝜔|𝜔| (3.1)
in which:
𝜔 = propeller angular velocity

𝛼, 𝛽 = constant model parameters
𝑄 = input torque
As can be derived from above equation, the propeller angular acceleration is depended on the differ-
ence between the torque delivered by the drive system and the torque required by the propeller. It
is assumed that the parameters 𝛼 and 𝛽 remain constant during operation. Smogeli [32], proposes a
more elaborate model that also includes a shaft friction term:
𝐼 𝜔̇ = 𝑘 𝑄 − 𝑄 (𝜔, 𝛽) − 𝑄 (𝜔) (3.2)
1
𝑄̇ = (𝑄 −𝑄 ) (3.3)
𝑇
in which:
𝑘 = gearbox ratio
𝑄 (𝜔, 𝛽) = propeller torque
𝑄 = commanded motor torque
𝑇 = time constant
𝑄 = motor torque
𝐼 = rotational inertia of motor, shaft and propeller
𝑄 (𝜔) = shaft friction term
𝛽 = hydrodynamic pitch angle
The shaft friction term is defined by:
𝑄 (𝜔) = sign(𝜔)𝑄 + 𝐾 𝜔 (3.4)
where 𝑄 and 𝐾 represent respectively the static shaft friction term and the linear friction term.
3.2.1. Entrained water added rotational moment of inertia

Smogeli [32] acknowledges the fact that there exists an added mass term due to the hydrodynamic
forces in phase with the propeller rotational acceleration. However, in his work this term is neglected.
This will give a reduced total system inertia, and hence faster dynamics. For conservatism, an empiri-
cal estimate of the inertia term is included in the total inertia of the thruster system 𝐼 .
According to Saunders [29], the entrained water moment of inertia is defined by:
𝐼 =𝐾⋅𝐼 (3.5)
where 𝐼 is the propeller moment of inertia and 𝐾 is a factor in the range of 0.25-0.50, where a factor
of 0.25 is often used as a representative value. This method has a big uncertainty. Therefore, the
3.2. Dynamic model 25
empirical estimates of the added inertia term are calculated based on the work of MacPherson [19] and
Schwaneke [30]. In their work, the added rotational inertia term of the entrained water is defined by:
𝐼 = 𝐶 𝜌𝐷 (3.6)
in which:
𝐶 = fitting parameter
𝜌 = water density
𝐷 = propeller diameter
The fitting parameter 𝐶 is defined according Schwanecke [30] by:
0.0703(𝑃/𝐷) 𝐸𝐴𝑅
𝐶 = (3.7)
𝜋𝑍
in which:
𝑃/𝐷 = propeller pitch-diameter ratio

𝐸𝐴𝑅 = propeller expanded area ratio
𝑍 = total number of propeller blades
MacPherson [19] proposes the following expression for the fitting parameter 𝐶 :
𝐶 = 𝐶 𝐸𝐴𝑅(𝑃/𝐷) − 𝐶 (3.8)
where 𝐶 and 𝐶 are fitting parameters given in Table 3.1.
Parameter Z=3 Z=4 Z=5 Z=6

𝐶 0.00477 0.00394 0.00359 0.00344
𝐶 0.00093 0.00087 0.00080 0.00076
Table 3.1: Fitting coefficients according MacPherson [19]
Since both empirical models are based on model test results carried out using different propeller types,
the average of the model results is used as a representative value for the inertia term.
According Smogeli [32], the total rotational moment of inertia can be calculated from a motor or pro-
peller point of view. It is decided to use the motor point of view, since this is the component that is
controlled. Using this approach, the total system rotational moment of inertia is defined by:
𝐼 +𝐼 +𝐼
𝐼 = +𝐼 (3.9)
𝑘
in which:
𝐼 = propeller rotational moment

𝐼 = propeller shaft rotational moment
𝐼 = rotor rotational moment
3.2.2. Four-quadrant model

As can be observed from Figure 3.2, the thrusters will experience oscillating inflow velocities due to the
rolling motion of the vessel. These oscillating inflow velocities will affect the thrust and torque produc-
tion of the thrusters.
The thrust and torque characteristics of a wide range of propellers and ducted propellers is tested
experimentally. The results are captured by so-called 𝐶 and 𝐶 diagrams. These diagrams show the
thrust (𝐶 ) and torque (𝐶 ) coefficients as a function of the ship’s advance velocity, see for example
Stapersma [34]. However, these diagrams do not capture the thrust and torque behavior when the
propeller is subjected to a reversed inflow, or when the propeller is in braking mode. MARIN [34] de-
veloped the so called four-quadrant model to extend the 𝐶 and 𝐶 diagrams to the full operating range
of a propeller. The four quadrants are divided as shown in Table 3.2.
Quadrant 𝛽 [deg] Ship speed 𝑉 Propeller speed 𝑛 Description

1 0-90 + + normal forward operation
2 90-180 - + forward thrust at negative speed
3 180-270 + - reversing propeller during stopping
4 270-360 - - normal astern operation
Table 3.2: Four quadrant operation of a propeller [34]
The hydroydynamic pitch angle 𝛽 is defined by:
𝑣
𝛽 = arctan( ) (3.10)
0.7𝜋 ⋅ 𝑛 ⋅ 𝐷
where 𝑛 and 𝑣 represent the rotational speed of the propeller and the propeller advance velocity. The
propeller advance velocity as a result of the vessel roll motion is calculated by:
𝑣 = 𝜙̇ √𝑧 + 𝑦 + 𝑣 (3.11)
where 𝑧 , 𝑦 , 𝜙̇ and 𝑣 represent the vertical distance from the thruster centre to the CoG of the vessel,
horizontal distance from the thruster centre to the CoG of the vessel, the roll velocity amplitude and the
current velocity. An impression of a four-quadrant model is given in Figure 3.3.
Figure 3.3: Four-quadrant diagram of Ka 4-70 thruster with different P/D ratios as found in Oosterveld [22]
3.3. Model verification 27
The torque and thrust produced by the propeller are defined by:
𝜋
𝑄 (𝛽, 𝑛) = 𝐶 (𝛽)𝜌(𝑣 + (0.7𝜋𝑛𝐷) ) 𝐷 (3.12)
8
𝜋
𝑇 (𝛽, 𝑛) = 𝐶 (𝛽)𝜌(𝑣 + (0.7𝜋𝑛𝐷) ) 𝐷 (3.13)
8
The torque and thrust coefficients are approximated by using a 20-th order Fourier series:
𝐶 (𝛽) ≈ ∑ [𝐴 (𝑘) cos(𝛽𝑘) + 𝐵 (𝑘) sin(𝛽𝑘)] (3.14)
𝐶 (𝛽) ≈ ∑ [𝐴 (𝑘) cos(𝛽𝑘) + 𝐵 (𝑘) sin(𝛽𝑘)] (3.15)
where 𝐴 , 𝐵 , 𝐴 and 𝐵 are Fourier coefficients based on experimental data. The values can be found
in Appendix B.
The roll moment induced by the thruster is calculated by multiplying the produced thrust by the roll
moment arm. The roll moment arm is defined according Rudaa [28] as the vertical distance between
the thruster propeller centre and the centre of gravity of the vessel, see Appendix B.
To take into account the effect of thruster-hull interaction, the calculated open water thrust is reduced
by applying a thrust deduction factor 𝑡 = 0.04. This is the thrust deduction factor found by Wichers,
[35] for an azimuth angle of 90∘ , see Appendix B.
3.3. Model verification

The dynamic thruster model has been applied to the azimuthing thrusters installed at the Ndurance.
The stern thrusters are of the thruster type TH1500 and the forward thrusters are type TH1250. The
thruster parameters of both types can be found in Appendix B. The results regarding thrust production
are compared to manufacturer thruster curve data to ensure the validity of the modelling approach.
Subsequently, the transient response of the thrusters is verified by expert indication.
3.3.1. Thruster curves

The implementation of the four-quadrant model is verified by comparing the thruster curves during dif-
ferent advance velocities as provided by the manufacturer with model output. The results are given in
Figure 3.4.
TH1250 TH1500
110 110
Manufacturer data Manufacturer data
100 100
V = 0 knts V = 0 knts
a a
90 V = 2 knts 90 V = 3 knts
a a
a a
80 80
a a
70 70
Thrust [%]
Thrust [%]
60 60
50 50
40 40
30 30
20 20
10 10
0 0
0 20 40 60 80 100 0 20 40 60 80 100
Thruster RPM [%] Thruster RPM [%]
Figure 3.4: Model results for different advance velocities compared to manufacturer data for thruster type 1250 and 1500
As can be observed in Figure 3.4, the model output is in good agreement with the manufacturer data.
To quantify the agreement, the relative mean absolute error between the model estimates and manu-
facturer data is calculated according Equation 2.16. The results are given in Table 3.3.
TH1250 RMAE [%] TH1500 RMAE [%]

𝑉 =0 knts 1.06 𝑉 =0 knts 0.48
𝑉 =2 knts 0.81 𝑉 =3 knts 1.42
𝑉 =4 knts 0.91 𝑉 =6 knts 0.76
𝑉 =5 knts 0.91 𝑉 =10 knts 0.95
Table 3.3: RMAE values for thruster type TH1250 and TH1500 during different advance velocities
From Table 3.3 can be observed that the maximum RMAE value is 1.42%. This is considered suffi-
ciently small to consider the mathematical model verified.
3.3.2. Transient response

The transient response of the thrusters is analyzed to asses whether the model is able to simulate
realistic thruster behavior. In order to do so, a step impulse from 30% to 100% thrust is applied to the
dynamic thruster model. The transient response of both thruster types is visualized in Figure 3.5.
TH1250 TH1500
110 110
100 100
90 90
80 80
70 70
Thrust [%]
Thrust [%]
60 60
50 50
40 40
30 30
20 Step load 20 Step load

Thruster response Thruster response
10 10
9 10 11 12 13 14 15 9 10 11 12 13 14 15
Time [s] Time [s]
Figure 3.5: Transient response for thruster type 1250 and 1500 models
Figure 3.5 shows that both thrusters take around 3 seconds to ramp from 30% to 100% thrust. The
behavior was hard to verify by literature, since no data was available for this specific thruster set-up.
The Chief Engineer of the Ndurance estimated that it would take about 2 seconds for the thrusters to
ramp from 30% to 100% thrust. This estimate is in agreement with the model results. Therefore, it is
assumed that the model is able to simulate realistic dynamic thruster behavior.
3.4. Conclusions 29
3.4. Conclusions
The second sub-research objective is defined as:
”Develop a thruster model that captures the thruster dynamics in the time domain, this will include
modeling of the thruster drive train system and thrust variations due to oscillating inflow velocities.”
A dynamic thruster model is presented in this chapter. An empirical estimate of the entrained wa-
ter added rotational inertia has been incorporated in the dynamic model to ensure conservatism. The
thrust and torque produced by the thruster during oscillating inflow velocities is modelled by using a
four-quadrant model.

• The four-quadrant thruster model thrust results have a maximum relative mean absolute error
value of 1.42% when compared to manufacturer data. This error is considered sufficiently small
to consider the mathematical model verified.
• The step load response of the dynamic thruster model is in line with the expectations of the
Ndurance Chief Engineer. From this result is concluded that the thruster model is able to simu-
late realistic behavior.
4
DP System Model
In this chapter a DP system model is presented. The DP system model is developed in order to simulate
the behavior of the DP system installed on board the Ndurance. By using this model it is possible to
investigate the DP footprint, DP capability and thruster power consumption of the ship.
4.1. Introduction
The horizontal motions surge, sway and yaw as calculated by the time domain model are based on
a body-fixed reference frame and earth-fixed reference frame. The vessel motions, velocities and
accelerations are calculated in the body-fixed frame and are subsequently translated to the earth-fixed
reference frame by using the rotation matrix:
cos(𝜓) − sin(𝜓) 0
R(𝜓) = [ sin(𝜓) cos(𝜓) 0 ] (4.1)
0 0 1
The reference frames used in this report are illustrated in Figure 4.1. Also the orientation of the envi-
ronmental directions is indicated.
μ=90°
u

μ=0°
x
r
v
y
Figure 4.1: Definition of reference frames and environmental force direction
The Ndurance is equipped with a DP system. The DP system is a control system that enables the ship
to maintain a specified position and heading when conducting operations in the offshore environment.
The installed DP system on board the Ndurance uses 5 thrusters, of which 4 azimuthing thrusters and
1 tunnel thruster, to maintain its position and heading. The thruster lay-out is given in Figure 4.2.
31
32 4. DP System Model
T1
T2 T3
CoG
T4 T5
Figure 4.2: Schematic visualization of the thruster system lay-out of the Ndurance
In this chapter, the components out of which a DP system consists are explained. Next to this, an
approach to model a DP system is presented and implemented. This modeling approach is largely
based on the work of Serraris [31]. The chapter concludes with a verification of the constructed DP
system model. This is done by comparing the model results to a DP time domain study conducted by
MARIN.
4.2. DP system configuration

A DP system consist out of several components. In general, the main components are:
• Position and heading sensors
• Kalman filter, to filter out the first-order ship motions in the sensor signal
• DP Controller, to estimate the required thrust to ensure the ship maintains heading and position
• Allocation algorithm, to divide the required thrust over the available thrusters
• Thrusters, to produce the required thrust for station keeping
Since a mathematical model is constructed, the ship motions as a result of the environmental forces
can not be observed by a sensor. Therefore, the sensor is left out and the calculated ship motions as
explained in Chapter 2 are directly fed to the Kalman filter. A schematic overview of the DP system
model is given in Figure 4.3.
Vessel
time domain model
Allocation Kalman
DP
filter
Figure 4.3: Schematic overview of DP system model
As can be observed in Figure 4.3, the symbol 𝜏 is used to indicate the DP controller. In the up-
coming sections, each system component of the DP system model and its modelling approach will be
discussed.
4.3. Environmental forces 33
4.3. Environmental forces

A ship in the offshore environment is subjected to forces and moments induced by waves, wind and
current. Since a mathematical model is constructed, these forces and moment need to be generated
numerically. For each component a model is presented.
4.3.1. Waves
The forces and moments induced by sea waves have two components: a first and a second-order
component. The first-order wave forces and moments are oscillatory with a zero mean. The first-order
wave forces result in first-order ship motions. These motions are those one experiences when sailing
in a seaway. The first-order forces and moments are calculated in the time domain according Journée
[20] by using:
𝐹 (𝑡) = ∑ √2 ⋅ 𝑆 , (𝜔)Δ𝜔 cos(𝜔 𝑡 + 𝜖 + 𝜖 , ) (4.2)
in which:
𝑁 = total number of wave frequencies

𝜖 = random phase angle in the domain [0 , 2𝜋]
𝜖 , = phase shift of the wave load
𝑆 , (𝜔) = wave force response spectrum
𝑡 = time vector
The wave force response spectrum 𝑆 , (𝜔) is defined by the wave force RAO squared multiplied by a
wave spectrum:
𝐹
𝑆 , (𝜔) = | (𝜔)| 𝑆 (𝜔) (4.3)
𝜉
The wave spectrum is defined by a JONSWAP spectrum as found in Journée [20]:
320𝐻 −1950
𝑆 (𝜔) = 𝜔 exp{ 𝜔 }𝛾 (4.4)
𝑇 𝑇
in which:
𝐻 = significant wave height

𝑇 = wave peak period
𝛾 = 3.3 (peak factor)
𝜔/𝜔 − 1
𝐴 = exp{( ) }
𝜎√2
𝜎 = a step function of 𝜔 ∶ if 𝜔 < 𝜔 then 𝜎 = 0.07

if 𝜔 > 𝜔 then 𝜎 = 0.09
The second-order wave forces and moments arise from the square of the oscillatory components of the
first-order wave components. This results in mean, high- and low-frequency wave forces and moments.
The magnitude of these forces is smaller compared to the first-order wave forces. However, they can
result in significant motions due to the low amount of damping at these frequencies and the fact that
the frequency of these forces and moments can be near the low natural frequency of a moored ship.
The mean component of the second-order wave force is commonly referred to as the mean drift force.
The low-frequency and mean wave drift forces in irregular waves are calculated according Newman
[21] in the time domain by using:
𝐹 (𝑡) = ∑ ∑ 𝜉 𝜉 𝑃 cos((𝜔 − 𝜔 )𝑡 + (𝜖 − 𝜖 )) + ∑ ∑ 𝜉 𝜉 𝑄 sin((𝜖 − 𝜖 )𝑡 + (𝜖 − 𝜖 )) (4.5)
in which:
𝜉 = 𝑖 wave amplitude
𝜉 = 𝑗 wave amplitude
𝜔 = 𝑖 wave frequency
𝜔 = 𝑗 wave frequency
𝑃 = in-phase quadratic transfer function (QTF)
𝑄 = out-of-phase quadratic transfer function (QTF)
𝑡 = time vector
The mean components are calculated when 𝑖 = 𝑗. Note that the high-frequency, or sum, components
are not incorporated, since these forces and moments are considered negligible for a softly moored
structure like a ship with a DP control system.
The QTFs 𝑄 and 𝑃 are obtained from ANSYS AQWA. The panel method uses the far-field method to
calculate the mean components of QTF and the near-field method to calculate the sum and difference
components of the QTF. For a more elaborate review of both methods, see Journée [20].
4.3.2. Wind
The part of the ship that is not submerged is subjected to forces and moments induced by the wind.
The wind forces and moments for surge, sway and yaw are calculated according to Serraris [31] by:
1
𝐹 = 𝐶 (𝜇) 𝜌 𝑣 𝐴 (4.6)
2
1
𝐹 = 𝐶 (𝜇) 𝜌 𝑣 𝐴 (4.7)
2
1
𝑀 = 𝐶 (𝜇) 𝜌 𝑣 𝐴 𝐿 (4.8)
2
in which:
𝐶(𝜇) = wind load coefficient
𝜇 = environmental direction
𝜌 = air density, typically 0.00125 ton/m
𝑣 = wind speed
𝐴 = frontal wind area
𝐴 = lateral wind area
𝐿 = length between perpendiculars
In this environmental force model it is assumed that the wind speed is constant. The wind load coef-
ficients can be obtained by conducting wind tunnel test or by the use of empirical models. The wind
coefficients of the Ndurance are given in Appendix C.
4.3.3. Current
The submerged part of the ship is subjected to forces and moments induced by the water current. The
environmental current force model for surge, sway and yaw are calculated in a similar approach as the
4.4. Kalman filter 35
wind force as stated in Serraris [31]:
1
𝐹 = 𝐶 (𝜇) 𝜌𝑣 𝐿 𝑇 (4.9)
2
1
𝐹 = 𝐶 (𝜇) 𝜌𝑣 𝐿 𝑇 (4.10)
2
1
𝑀 = 𝐶 (𝜇) 𝜌𝑣 𝐿 𝑇 (4.11)
2
in which:
𝐶(𝜇) = current load coefficient

𝜌 = sea water density, typically 1.025 ton/m
𝑣 = current speed
𝑇 = draft
Similar to the wind speed, it is assumed that the current speed is constant. The current load coeffi-
cients can be obtained by conducting wind tunnel test or by the use of empirical models. The current
coefficients of the Ndurance are given in Appendix C.
4.4. Kalman filter

A Kalman filter is used in a DP system to filter out the first-order motions of the ship. The filtering is
done in order to make sure that the thrusters do not respond to the high frequency motions. This would
result in increased fuel consumption and wear and tear of the thrusters.
The benefits of a Kalman filter are minimum phase lag and the ability to continue estimation of a signal
when the sensor signal is lost. The last property can prevent the DP system from a run-off in case of
a sensor failure.
The Kalman filter estimates the vessel position by calculating the response of the ship according to
its mass and the low-frequent excitation force. This force is derived from the measured thruster forces.
The estimated response is compared to the measured response and a weighted portion of the differ-
ence between the position estimate and the measured position is added to the signal. The Kalman
filter is based on a prediction-correction algorithm. First the prediction step is carried out as defined by
Cadet [4]:
𝑥̂ = A𝑥̂ + B𝑢
{ (4.12)
P = AP A +Q
in which:
𝑥̂ = state that is estimated at timestep k

𝑢 = control input
A = state transition matrix
B = control input matrix
P = error covariance matrix at timestep k
Q = process noise covariance matrix
For every motion that is controlled by the DP system, a seperate Kalman filter has been implemented.
In the report the yaw motion is described, but exactly the same approach has been used for the surge
and sway motion. For the yaw motion, the estimated state vector, state transition matrix, control input
vector and process noise matrix become:

𝜓
𝑥̂ = [ ̇] (4.13)
𝜓
0
𝑢 = [ ] (4.14)
𝑀
1 Δ𝑡
A = [ ] (4.15)
0 1
1 0
B = [0 ] (4.16)
( )
𝑄 0
Q = [ ] (4.17)
0 𝑄
In Equation 4.14, 𝑀 represent the yaw moment induced by the thrusters. The 𝑄 term in Equation 4.17
represents the process noise covariance. This is a tuning parameter of the Kalman filter. A high Q
value will put more emphasis on the measurement signal.
Next is the correction step:

K = P H (HPH + R)
{ 𝑥̂ = 𝑥̂ + K (𝑧 − H𝑥̂ ) (4.18)
P = (I − K H)P
in which:
𝑧 = measurement signal
K = Kalman gain
H = mapping matrix from measurement to state
I = 2x2 identity matrix
R = measurement noise covariance matrix
For the yaw motion, the measurement signal vector, mapping matrix and measurement noise co-
variance matrix become:
𝜓
𝑧 = [ ̇] (4.19)
𝜓
1 0
H = [ ] (4.20)
0 1
𝑅 0
R = [ ] (4.21)
0 𝑅
In Equation 4.21, 𝑅 represents the measurement noise covariance, which is also a tuning parameter.
When the noise covariance is set at a high value, the Kalman filter puts less emphasis on the measure-
ment signal.
The Kalman filter uses a recursive method to estimate the filtered signal. Therefore, it needs to be
initialized with first estimates of the state and error covariance matrix. The measurement and process
noise covariances are tuned manually to obtain the desired filter performance. Tuning of these param-
eters affect the over-all performance of the DP system. When the Kalman filter is tuned to be very
responsive, the thruster response will also be responsive, resulting in a decreased DP footprint. More
relaxed Kalman tuning parameters will cause the thrusters to react more slowly, but will also increase
the DP footprint. DP operations like float-overs require a small DP footprint, where drilling ships or
FPSOs allow a bigger DP footprint. Therefore, it is important to keep in mind the desired performance
of the DP system when tuning the Kalman filter.
4.5. DP controller 37
4.5. DP controller
After the Kalman filter has been applied, the position and heading signals are send to the DP controller.
The DP controller consists out of 3 stand-alone closed loop PID controllers, responsible for controlling
the surge, sway and yaw motions of the vessel. An schematic example of a standard closed-loop PID
controller is given in Figure 4.4.
Figure 4.4: Block diagram of closed-loop PID controller
The PID controllers determine the total thrust in the horizontal plane required by the propulsion system
to maintain position. The required thrust is calculated according to Serraris [31] by:
𝜏 = 𝐾 𝑒 + 𝐾 ∫ 𝑒(𝑡)𝑑𝑡 + 𝐾 𝑒̇ (4.22)
in which:
𝐾 = proportional coefficient vector

𝐾 = integral coefficient vector
𝐾 = derivative coefficient vector
𝑒 = position error vector
The position error is defined by:

𝑥 −𝑥
𝑒 =[ 𝑦 −𝑦 ] (4.23)
𝜓 −𝜓
where 𝑥 , 𝑦 and 𝜓 represent the surge, sway and yaw setpoint during the DP operation. The variables
𝑥, 𝑦 and 𝜓 represent the actual vessel position and heading. Note that the PID coefficients are vectors,
since the coefficient value is not equal for each motion component. The proportional coefficient 𝐾
determines the stiffness of the DP system. This is comparable to the spring stiffness of a mass-spring-
damper system. The derivative 𝐾 coefficient represents the damping of the DP system, likewise the
damping term in a mass-spring-spring system. The integral term 𝐾 is the integration coefficient. This
coefficient enables compensation of the mean error offset of the system by integrating the error over
time.
The P, I and D coefficients determine the behavior of the DP system. The coefficient values are tunable.
Similarly to the tuning of the Kalman filter, it is important to keep in mind the desired performance of
the DP system when tuning the PID coefficients. The initial values are often obtained by rule of thumb
estimates and manually tuned to obtain the desired system behavior. Serraris [31] recommends the
following rule of thumb estimates for the P and D coefficients:
𝑁 ⋅𝑇
𝑃= (4.24)
(50 − 70%)𝑅
𝐷 = (50 − 70%) ⋅ 2√(𝑀 + 𝑎)𝑃 (4.25)
in which:
𝑁 = total number of thrusters

𝑇 = maximum thrust per thruster
𝑅 = maximum allowable horizontal excursion
𝑀 = ship mass
𝑎 = ship added mass
These rules of thumb provide a first estimate of the P and D coefficients. The I coefficient is often set to
zero, since model tests showed that instability of the system might occur according Serraris [31]. Next
to this, the mean offset, which is reduced by the I coefficient, can be neglected when one is interested
in the station keeping accuracy of the DP system.
4.6. Thrust allocation

The total thrust in the horizontal plane as calculated by the DP controller is fed to a thrust allocation
algorithm. Thrust allocation serves multiple purposes. The main purpose is to divide the required thrust
over the available thrusters in such a way that the power consumption is minimized. Next to this, the
thrust allocation algorithm also prevents thrusters from operating in so called ’forbidden zones’. For-
bidden zones are the azimuth angles at which a thruster wake blows in the direction of a neighbouring
thruster. When this happens, the efficiency of the downstream thruster is significantly reduced, which
can result in a run-off of the DP system. Finally, the thrust allocation algorithm makes sure that the
thrusters of the DP system are not overloaded, by making sure that no additional thrust is allocated to
thrusters that have reached their saturation limit.
The thrust allocation algorithm tries to find the minimum of a objective function while satisfying both
linear and non-linear constraints. The thrust allocation algorithm can be mathematically formulated by:
minimize 𝐹 (𝑥)
subject to 𝑐(𝑥) ≤ 0 (4.26)

A⋅𝑥 =𝑏
in which:
𝐹 (𝑥) = objective function

𝑥 = unknown thrust vector in x- and y-direction
𝑐(𝑥) = vector with non-linear constraints
A = matrix with linear constraints
𝑏 = required forces and moment in the horizontal plane
The objective function is defined as described in Serraris [31]:
𝑇 𝑇
𝐹 (𝑥) = ∑ | | + ∑ 𝑤 (( ) − 1) (4.27)
𝑇 𝑇
The last summation term is a penalty function that is included in the objective function when a thruster
saturates. In that case 𝑤 = 10000, else the value of 𝑤 equals 0.
The total number of unknowns is calculated by:
𝑁 = 2𝑁 +𝑁 (4.28)
In this equation 𝑁 represents the number of azimuth thrusters and 𝑁 the number of bow thrusters
that are part of the DP system. Since an azimuth thruster can deliver thrust in both the x- and y-direction,
4.6. Thrust allocation 39
it contributes to 2 unknown thrust components. Since the Ndurance has 4 azimuth thrusters and 1 bow
thruster, the total number of unknowns equals 9. The unknown thrust vector becomes:
𝑇
⎡ ⎤
⎢𝑇 ⎥
⎢𝑇 ⎥
⎢𝑇 ⎥
𝑥 = ⎢𝑇 ⎥ (4.29)
⎢𝑇 ⎥
⎢ ⎥
⎢𝑇 ⎥
⎢𝑇 ⎥
⎣𝑇 ⎦
In this vector, the term in the subscript corresponds to the thruster number and the thrust direction. The
vector with non-linear constraints makes sure that the function is minimized without thruster overload-
ing. The constraints are defined by:
⎡√𝑇 +𝑇 −𝑇 ⎤
⎢ ⎥
⎢√𝑇 +𝑇 −𝑇 ⎥
𝑐(𝑥) = ⎢ ⎥ (4.30)
⎢√𝑇 +𝑇 −𝑇 ⎥
⎢ ⎥
⎢ +𝑇 −𝑇 ⎥
√𝑇
⎣ ⎦
Note that these constraints are non-linear due to the root term. The linear constraints matrix A is
defined by Serraris [31]:
0 1 0 1 0 1 0 1 0
A = [ 1 0 1 0 1 0 1 0 1 ] (4.31)
𝑥 −𝑦 𝑥 𝑦 𝑥 −𝑦 −𝑥 𝑦 −𝑥
The terms in the third row of A represent the thruster x- and y-position with regards to the centre of
gravity of the ship. These terms have been added to calculate the moments created by each thruster.
The required forces and moments in the horizontal plane that the thruster system has to realize to
maintain position are captured in vector 𝑏:
𝐹
𝑏 = [𝐹 ] (4.32)
𝑀
Vector 𝑏 is directly fed from the DP controller and is used as the input variables of the allocation al-
gorithm. The allocation algorithm uses the ’fmincon’ Matlab function as solver. The mathematical
background of the non-linear programming method, which forms the basis for the fmincon function, is
elaborately discussed in Byrd [3].
The forbidden zones are taken into account by using a similar approach as proposed by Serraris [31].
The algorithm checks whether the solution contains azimuth angles located in a forbidden zone. When
this is the case, the algorithm fixes the forbidden azimuth angle to the nearest boundary of the forbidden
zone and reruns the allocation. This procedure is repeated until the solution satisfies the constraints
and no azimuth angles are located in forbidden zones. The performance of the algorithm is checked
by subjecting the vessel to environmental forces from 0∘ to 360∘ . The results are given in Appendix C.
The numerical values of the forbidden zones are given in Table 4.1.
Thruster Forbidden zone [deg]

2 256.4 - 283.6
3 76.4 - 103.6
4 255.6 - 284.4
5 75.6 - 104.4
Table 4.1: Thruster azimuth angle forbidden zones
The allocation algorithm is executed during every time step. As a result of this, it is possible that the
allocated azimuth angle differs significantly between two time steps. In reality a thruster is not able to
switch its azimuth angle instantly. According to Serraris [31], it takes an azimuth thruster around 30
seconds to turn 360∘ . Therefore, the thruster rotation rate is limited to 12∘ per second.
4.7. Model verification

Now the building blocks of a DP system model have been constructed, the blocks are connected and
merged into a complete DP system model. The constructed DP system model needs to be verified.
This is done by comparing the model output to DP time domain simulations conducted by MARIN.
MARIN uses its in-house developed simulator AnySim to conduct time domain DP simulations. The
underlaying principles and assumptions of AnySim are presented in Serraris [31].
Three simulation cases have been selected for comparison. The selected cases are given in Table
4.2.
Case Draft Vw 𝜇 Vc 𝜇 Hs Tp 𝜇 Spectrum

[#] [m] [m/s] [deg] [m/s] [deg] [m] [s] [deg] [-]
1 4.7 12.86 45 1.286 45 2.0 5.7 45 JONSWAP
2 4.7 12.86 180 1.286 180 2.0 5.7 180 JONSWAP
3 4.7 12.86 90 1.286 165 2.0 5.7 165 JONSWAP
Table 4.2: Overview of simulation cases and corresponding environmental conditions for model verification
The ship, thruster and system tuning parameters used are given in Appendix C.
One should note that there exist many differences between the constructed DP model and the MARIN
model. The differences considered the most significant are:
• Kalman filter, the MARIN model uses a coupled non-linear state transition model. Therefore, an
Extended Kalman filter is used. The DP system model uses a linear Kalman filter. Next to this, it
is unknown how the filter is tuned in the MARIN model.
• Allocation algorithm, the MARIN model uses a Lagrange multiplier method to solve the alloca-
tion problem, whereas the DP system uses a non-linear trust-region method to solve the allocation
problem.
• Ship motions, it is unknown which motion coupling terms have been included in the MARIN
model. Also, it is not known whether or not the MARIN model includes a viscous roll damping
term.
• Wave forces, the calculation method of wave forces in the time domain that is used by the MARIN
model is unknown.
The three selected cases have been simulated by the constructed DP system model. The mean, stan-
dard deviation, maximum value and minimum value of the calculated time signal of the ship motions,
thrust and azimuth angles are given in Table 4.3, Table 4.4 and Table 4.5.
4.7. Model verification 41
Case 1 MODEL MARIN

Motion Mean Std Max Min Mean Std Max Min
Surge [m] 0.46 0.12 0.69 -0.34 0.55 0.64 5.56 -1.13
Sway [m] 1.38 0.15 1.69 0.00 1.73 1.11 12.03 0.00
Heave [m] 0.00 0.05 0.19 -0.19 0.00 0.04 0.18 -0.15
Roll [deg] 0.00 0.18 0.65 -0.68 -0.11 0.26 1.03 -1.21
Pitch [deg] 0.00 0.27 1.11 -1.08 0.00 0.27 1.15 -1.08
Yaw [deg] -0.51 0.15 0.09 -1.08 -0.66 0.79 3.28 -5.06
Thrust
T1 [kN] -25.01 4.94 0.00 -42.23 -42.79 9.37 0.69 -59.93
T2 [kN] 64.27 11.25 104.95 0.01 59.58 20.32 128.05 0.65
T3 [kN] 61.25 12.44 99.81 0.01 57.04 20.51 130.38 0.65
T4 [kN] 96.97 13.04 181.38 0.00 93.37 22.03 183.30 1.12
T5 [kN] 95.49 13.75 183.69 0.01 96.42 23.74 198.68 1.12
Azimuth
T1 [deg] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
T2 [deg] -104.68 10.21 0.00 -141.98 -116.68 5.19 0.00 -152.22
T3 [deg] -107.79 9.87 0.00 -136.44 -109.97 3.90 0.00 -147.53
T4 [deg] -106.79 13.80 0.00 -148.16 -106.38 6.29 0.00 -167.86
T5 [deg] -110.24 13.94 0.00 -165.32 -103.69 8.42 63.15 -155.78
Table 4.3: Comparison of simulation results for case 1
Case 2 MODEL MARIN

Surge [m] -0.488 0.064 0.001 -0.659 -0.697 0.439 0.210 -2.997
Sway [m] 0.000 0.000 0.002 0.000 0.000 0.021 0.481 -0.216
Heave [m] 0.000 0.038 0.121 -0.132 0.000 0.057 0.211 -0.212
Roll [deg] 0.000 0.000 0.000 0.000 0.000 0.002 0.097 -0.037
Pitch [deg] 0.000 0.161 0.636 -0.681 -0.002 0.205 0.839 -0.895
Yaw [deg] 0.072 0.007 0.077 0.000 -0.005 0.035 0.462 -0.630
Thrust
T1 [kN] -0.827 0.062 -1.056 0.000 0.007 0.613 11.220 -17.300
T2 [kN] 19.256 5.670 60.718 0.004 27.201 5.601 57.101 -2.200
T3 [kN] 18.733 5.663 60.344 0.004 27.201 5.598 56.571 -2.200
T4 [kN] 34.142 10.106 107.947 0.004 23.596 5.337 51.772 -4.006
T5 [kN] 35.193 10.095 107.292 0.004 23.596 5.331 51.409 -4.006
Azimuth
T1 [deg] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
T2 [deg] -0.109 0.155 -0.019 -3.142 0.008 4.567 179.892 -178.233
T3 [deg] -0.112 0.156 -0.019 -3.142 0.008 4.567 179.892 -178.223
T4 [deg] 0.069 0.129 3.096 0.000 -0.007 4.594 179.861 -178.264
T5 [deg] 0.072 0.130 3.097 0.000 -0.007 4.593 179.861 -178.264
Several conclusions can be drawn from the comparison of the results given in Table 4.3, Table 4.4 and
Table 4.5. First the ship motion results are considered. Next, the thrust results and finally the resulting
azimuth angles.
Ship motion
The mean values of the ship motions are in good agreement for all cases. The maximum deviation is
0.35 meter, which is considered negligible for a ship with a length of 99 meters. From this can be con-
cluded that the mean environmental forces, caused by the wave drift force, wind force and current force
are generally in agreement in both models. However, the standard deviation, maximum and minimum
Case 3 MODEL MARIN

Surge [m] -0.272 0.085 0.111 -0.561 -0.299 0.501 2.045 -1.941
Sway [m] 1.340 0.143 1.555 0.000 1.680 0.320 3.294 0.000
Heave [m] 0.000 0.038 0.118 -0.112 0.000 0.057 0.218 -0.211
Roll [deg] 0.000 0.068 0.279 -0.291 -0.093 0.110 0.333 -0.567
Pitch [deg] 0.001 0.169 0.568 -0.574 -0.001 0.202 0.839 -0.841
Yaw [deg] 0.252 0.062 0.425 -0.097 0.240 0.292 1.777 -0.671
Thrust
T1 [kN] -30.291 3.734 0.404 -49.928 -41.585 1.988 1.125 -44.542
T2 [kN] 73.299 8.520 117.535 0.000 62.892 7.312 109.486 0.000
T3 [kN] 62.497 8.015 104.139 0.000 62.988 7.474 108.111 0.000
T4 [kN] 65.557 9.890 101.525 0.000 61.683 7.676 106.247 0.000
T5 [kN] 60.118 9.224 95.373 0.000 62.223 7.710 107.740 0.000
Azimuth
T1 [deg] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
T2 [deg] -79.771 15.960 0.000 -179.332 -80.487 6.275 0.000 -106.470
T3 [deg] -82.344 18.119 0.000 -179.259 -75.329 2.418 0.000 -103.600
T4 [deg] -80.045 18.675 0.000 -153.712 -81.033 6.065 0.000 -102.928
T5 [deg] -80.059 19.265 0.000 -153.722 -75.001 2.081 0.000 -104.400
value of the calculated ship motions differs significantly. The tuning of the Kalman filter, simulation time
and initial conditions can be reasons for these deviations. In case 1, a drift-off of 12 meters in the sway
direction is observed in the MARIN results. This would be unacceptable during a real DP operation.
The DP model shows much lower maximum and minimum values, resulting subsequently in a lower
standard deviation.
Thrust
In general the mean thrust values are in good agreement. However, the MARIN model allocates sig-
nificantly more thrust to the bow thruster. The reason is unknown, but probably this is the result of a
different allocation algorithm. During case 2, the MARIN model allocates more thrust to the forward
thrusters compared to the aft thrusters. This is unexpected, since the aft thrusters have a bigger thrust
capacity and it is therefore powerwise more efficient to use the aft thrusters to compensate the head
waves, wind and current. The DP model behavior is more as expected.
Azimuth angles
The mean azimuth angles are generally in agreement. However, the standard deviations tend to de-
viate. During case 2, the MARIN model has a significant higher standard deviation, whereas the DP
model shows higher standard deviations during the other cases. Reasons for this different behavior
could be the Kalman filter tuning and the different allocation algorithm.
To quantify the agreement between both models, the RMAE value of the horizontal motions and the
thrust and azimuth angle mean estimates are calculated. The results are given in Table 4.6.
Parameter RMAE [%]

Surge 8.7
Sway 13.3
Yaw 8.7
Thrust 6.7
Azimuth 2.1
Table 4.6: RMAE values for model results compared to MARIN results
4.8. Conclusions 43
4.8. Conclusions
The third sub-research objective is defined as:
”Develop a simulation model that is able to simulate DP system behavior, this includes modelling of the
DP system and environmental forces.”
A DP system model consisting out of a Kalman filter, DP controller and a thrust allocation algorithm is
developed. The DP model is merged with the vessel time domain model. Different models to calculate
the environmental forces due to wind, waves and current are presented.

• The developed DP model is able to simulate stable DP station keeping.
• The maximum relative mean absolute motion error when the simulation results are compared to
MARIN simulations is 13.3% for the mean sway motion. This error is equal to 0.35 m in the sway
direction, which is considered acceptable.
• The relative mean absolute error of the thrust and azimuth angle mean simulation results, when
compared to MARIN data, is respectively 6.7% and 2.1%. This is considered acceptable.
• The significant differences in standard deviation, maxima and minima between the constructed
DP model and the MARIN model are most probably the result of different Kalman filter tuning
coefficients and a different thrust allocation algorithm.
• The developed DP system model merged with the environmental force models and the vessel
time domain model is able to simulate DP system behavior.
5
3DP System Model
In this chapter, a novel control model is proposed for combined thruster induced roll reduction and DP
station keeping. Also, the tuning procedure of the controllers is discussed.
5.1. Introduction
The main objective of this research is to investigate whether it is possible to reduce the significant roll
amplitude of the Ndurance in beam seas during DP operations by active control of the thrusters.
Koschorrek et al. [17] proposed a method to combine DP and roll control for a vessel equipped with
Voith Schneider propellers. Sørensen and Strand [33] and Xu et al. [36] proved that it is possible to
actively damp the unintentional roll-pitch motion induced by the DP system of a semi-submersible. The
proposed control algorithm also increased the DP position keeping performance. A combined control
model for first-order roll reduction combined with DP station keeping using conventional thrusters has
not been investigated to the best of the author’s knowledge.
The development of such a control model is presented in this chapter. In the previous chapters, mod-
els to calculate the vessel motions, thruster response and DP system behavior have been presented.
These models combined provide a mathematical playground, which is used to develop and implement
a novel control model system that is able to combine DP station keeping and active roll reduction. A
general control structure is proposed and a novel allocation approach is presented. Also, the tuning of
the control systems in the model is discussed.
5.2. Conventional DP control system vs. 3DP control system

To achieve DP station keeping and active roll reduction combined, a new control model needs to be
implemented. The differences between a conventional DP control system and the anti-roll control
system are summarized in Table 5.1.
DP control system 3DP control system

System objective station keeping roll reduction and station keeping
System application all conditions beam seas
Controlled motions x,y,𝜓 x,y,𝜓, 𝜙
Controller bandwith low frequency high frequency
Allocation objective power minimization roll reduction maximization
Table 5.1: Differences between a conventional DP model and the 3DP model
From Table 5.1 can be observed that in order to develop an anti-roll control system, the allocation and
controller part of a conventional DP control system needs to be adapted. The newly proposed allocation
and controller models are presented in the following sections.
5.3. Allocation model

In a conventional DP system, the allocation algorithm allocates the required thrust to the available
thrusters in such a way that the thruster power consumption is minimized. During active roll reduction,
allocation based on minimization of power consumption is not applicable anymore, since a lot of thruster
power will be used to achieve roll reduction. Therefore, a new allocation method is proposed to enable
both roll reduction and station keeping.
45
46 5. 3DP System Model
5.3.1. Roll reduction

Rudaa [28] proposes to use thruster pairs to counteract the wave induced roll moment. When the
thruster pairs are fixed pointing in opposing direction, this yields three main advantages:
• Yaw stability, when the thruster pair consist out of a forward and aft thruster, the yaw moments
induced by both thrusters are balanced.
• Sway stability, when both a port side and starboard thruster pair is used, the sway motion induced
by the port side thruster pair is balanced by the starboard thruster pair.
• No thruster interaction, since the thruster wakes are pointed away from each other, there is no
possibility of thruster-thruster interaction
A schematic visualization of the thruster pairs and wave direction is given in Figure 5.1.
T1
T2 T3
CoG
T4 T5
Figure 5.1: Thruster pair configuration 3DP control model for roll reduction
The first thruster pair consist out of T2 and T4. The second thruster pair consist out of T3 and T5. As
can be observed from Figure 5.1, no thruster-thruster interaction can occur, since all thruster wakes
are pointed away from each other. To obtain maximum roll reduction, the thruster pairs need to work
in anti-phase. This means that when thruster pair 1 is ramping up to counteract the wave induced
roll moment, thruster pair 2 is ramping down. This method of operation also ensures yaw and sway
stability, since the averaged yaw moment and sway force produced by both thruster pairs over time
equals zero when operating perfectly in anti-phase. The allocation vector H is used to indicate which
thrusters are used for roll reduction:
0
⎡ ⎤
⎢ 1 ⎥
H =⎢ 1 ⎥ (5.1)
⎢ 1 ⎥
⎣ 1 ⎦
From H can be observed that all azimuthing thruster are used for roll reduction purposes.
5.3.2. Station keeping

Next to roll reduction, the thrusters are also used for station keeping. In beam waves, the wave forces
that need to be counteracted by the control model are mainly in the sway direction. Since thruster T2
and T4 are already pointed in the direction of the incoming wave force, it is a logical approach to use
these thrusters for station keeping in the sway direction.
The yaw motion is controlled by the bow thruster, this makes sense since this thruster has a long
yaw moment arm and is not used for roll reduction purposes. Therefore, the complete thruster capacity
can be used to control the yaw motion.
Since all thrusters are aligned in the sway direction, compensation of environmental forces in the surge
direction is not possible. To solve this problem, the azimuth angles of the thrusters used for surge
control are controlled by an azimuth controller. Thruster T3 is used for compensation of positive surge
5.3. Allocation model 47
forces and thruster T5 is used to compensate negative surge forces. A schematic visualization of the
thrusters used for station keeping is given in Figure 5.2.
T1 T1 T1
T2 T3 T2 T3 T2 T3
L2
CoG CoG CoG
L4 T4 T5 T4 T5 T4 T5
+
Figure 5.2: Thruster configuration 3DP control model for station keeping (sway, surge and yaw motion)
The allocation vector H is used to indicate which thrusters are used to control the surge, sway and
yaw motions:
0 0 1
⎡ ⎤
⎢ 0 0 ⎥
H =⎢ 𝑋 0 0 ⎥ (5.2)
⎢ ⎥
0 0
⎢ ⎥
⎣ 𝑋 0 0 ⎦
where 𝐿 and 𝐿 represent the thruster moment arms of respectively thruster T2 and thruster T4, see
Figure 5.2. 𝑋 and 𝑋 represent two Boolean variables defined by:
1 if 𝑇 ≤ 0
𝑋 ={
0 otherwise
1 if 𝑇 > 0
𝑋 ={
0 otherwise
where 𝑇 represent the required thrust in the surge direction to maintain position as calculated by
the surge DP controller. Allocation matrix H indicates that thruster T1 is used for yaw control. The
required thrust in the sway direction is divided over thruster T2 and T4. The division is done according
the thruster yaw moment arms. By doing so, the yaw moment of both thruster is balanced and no
resulting yaw moment is induced. Negative thrust in the surge direction is delivered by T2, whereas T5
delivers positive thrust in the surge direction when required.
5.4. Controllers
The proposed control model consist out of a high-level motion controller and low-level shaft speed and
azimuth angle controllers. The control structure is schematically visualized in Figure 5.3.
Vessel
time domain model
azi Kalman
DP
filter
Thruster
s Allocation
model
Figure 5.3: 3DP control model structure
In this section the high-level motion controllers and the low-level shaft speed controller and azimuth
angle controller are discussed.
5.4.1. Roll reduction

The purpose of the anti-roll controller is to actively reduce the roll motion of the vessel. The high-level
motion controller therefore needs to counteract the first-order roll moment induced by the waves. The
controller output is a thruster RPM command, since this is the input set point variable for the low-
level shaft speed controller. Rudaa [28] proposes to use a non-linear proportional derivative (PD) type
controller:
|𝜙|̇
̇ √
𝜏 = −𝐾 sign(𝜙) − 𝐾 𝜙̈ (5.3)
𝑐
where 𝐾 , 𝑐 and 𝐾 represent the proportional coefficient, a tuning parameter and the derivative co-
efficient. Rudaa [28] includes the non-linear square root term to make sure the contribution from the
thrusters is equally large for all roll rates. Also the thrusters are limited to a specific limit, to increase
their reaction time. Systematic numerical simulations are carried out to analyze the influence of the
tuning coefficients and the thruster limits.
Tuning coefficients
In order to analyze the influence of the coefficients 𝐾 , 𝑐 and 𝐾 onto the controller performance, sys-
tematic numerical simulations have been carried out with different tuning coefficients at one particular
sea state. The selected simulation cases and corresponding tuning coefficients are visualized in Figure
5.4.
105
Kd
100
105
105
0
10
100
Kp c
Figure 5.4: Selected search space for controller tuning coefficients in range [0.01 10000], . m and . s
5.4. Controllers 49
The roll reduction for every simulation case is calculated by subtracting the RMS value of the roll angle
time trace with active roll reduction by the RMS value of the roll angles without active roll reduction.
The RMS value of the roll angle time trace is calculated by:
1
𝜙 = √ (𝜙 + 𝜙 + ... + 𝜙 ) (5.4)
𝑛
where 𝑛 is the total number of samples and 𝜙 is the roll angle amplitude.
To determine the sensitivity of the tuning parameters to the roll reduction, the maximum partial deriva-
tive of the RMS reduction is calculated with respect to every tuning coefficient in the calculated range:
𝜕𝑅𝑀𝑆 𝜕𝑅𝑀𝑆 𝜕𝑅𝑀𝑆

𝜕𝑐 𝜕𝐾 𝜕𝐾
The other tuning parameters are considered equal 1 when calculating the partial derivatives. The result
is given in Figure 5.5.
RMS roll reduction Maximum sensitivity

0.6 0.45
c
K 0.4
p
0.5 K
d
0.35
Max. partial derivative [-]
RMS roll reduction [deg]
0.4 0.3
0.25
0.3
0.2
0.2 0.15
0.1
0.1
0.05
0 0
10-2 100 102 104 c Kp Kd
Tuning coefficient [-]
Figure 5.5: Roll reduction for different tuning coefficients and maximum sensitivity of every tuning coefficient ( . m and
. s)
As can be observed in Figure 5.5, the effect of the tuning variables 𝑐 and 𝐾 is of a negligible magnitude
in the current controller model. The roll reduction depends on 𝐾 solely, therefore it is decided to use
a much more simplified proportional anti-roll controller 𝜏 defined by:
𝜏 = 𝐾 𝜙̇ (5.5)
𝜏 ≤𝜏 ≤𝜏
where 𝜏 and 𝜏 represent the thruster RPM limits. In order to analyze the optimum value of the
controller proportional coefficient 𝐾 , systematic numerical simulations have been carried out for three
different sea states and increase values of 𝐾 . The magnitude of the roll reduction is quantified by
calculating the relative roll reduction (RRR):
𝜙 −𝜙
𝑅𝑅𝑅 = (5.6)
𝜙
where 𝜙 is the RMS value of the roll angle time trace without active roll reduction and 𝜙 is the
RMS value of the roll angle time trace with roll control activated. The results are given in Figure 5.6.
100
H =0.25m
s
Hs=0.50m
80
H =0.75m
s
RMS roll reduction [%]

60
40
20
0
0 0.5 1 1.5 2 2.5
Kp 106
Figure 5.6: Roll reduction percentage for different values ( . s)
From Figure 5.6 can be observed that the roll reduction percentages converge to a constant value
when the coefficient 𝐾 increases. As can be derived from Figure 5.6, the roll reduction percentage
converges quicker for a high sea state. This can be explained by the fact that controller reaches the
saturation limit 𝜏 earlier as a result of the higher roll velocities during increased environmental loads.
On the basis of these results, it is decided to use a controller coefficient of a value of 𝐾 = 2.5 ⋅ 10 ,
since at this value, the roll reduction has been converged to a maximum value for all sea states.
Thruster limits
The thruster RPM is limited to a maximum and a minimum value. The maximum RPM value is incorpo-
rated to avoid overloading the thruster above its maximum rated RPM. The minimum RPM is included
to increase the response time of the thruster. Ramping up a thruster from 0% RPM to 90% RPM takes
more time then, for example, ramping from 15% RPM to 90% RPM. See Figure 5.7.
Figure 5.7: Response time of thruster TH1250 and TH1500 for different minimum limit RPMs
On the other hand, increasing the minimum RPM of the thruster will decrease the magnitude of the
counteracting roll moment that a thruster pair can realize. This is the result of the thruster pair that
is working in anti-phase at a particular minimum RPM, decreasing the net counteracting roll moment.
Therefore, there exist a trade-off between thruster response time and the minimum thruster RPM vari-
able. To investigate the optimum minimum thruster RPM limit, systematic numerical simulations have
been carried out for different minimum RPM limits and sea states. The results are given in Figure 5.8.
5.4. Controllers 51
100
Hs=0.25m
90
H =0.50m
s
80 H =0.75m
s
70
60
50
40
30
20
10
0
0 5 10 15 20 25 30
min
[%]
Figure 5.8: Roll reduction percentage for different as percentage of ( . s)
Figure 5.8 indicates that the optimum lower limit of the thruster RPM is 5% for all three sea states.
A clear trend is visible that when the lower RPM limit is increased, the roll reduction percentage de-
creases. This is as expected, since the net counteracting roll moment produced by the thruster pairs
is decreased. In the simulations, the lower RPM limit 𝜏 is set at 5% of 𝜏 .
5.4.2. Station keeping

The controller design for DP station keeping is based on the conventional PID controller structure as
proposed by Serraris [31]. However, this controller output is a thrust command. Since the controller
structure is based on RPM as control variable, the DP station keeping controller thrust estimate is
transformed to a RPM command by:
1 𝜏∗ ,
𝜏 , = √ −𝑣 (5.7)
0.7𝜋𝐷 𝐷 𝐶 𝜌
where 𝜏∗ , represents the thrust-based controller as described in the previous chapter.
The next step is to combine the output of the roll controller and station keeping controller. This is
done by using the allocation matrices H and H :
𝜏 =H 𝜏 +H 𝜏 (5.8)
𝜏 ≤𝜏≤𝜏
By using this approach, the high-frequency roll controller and the low-frequency DP controller are com-
bined. The principle is schematically visualized in Figure 5.9.
Max
RPM
+ DP
Min
Time
RPM
DP
Time
RPM
Time
Figure 5.9: 3DP control system high-level motion control principle
As can be observed from Figure 5.9, the high- and low-frequency controller signal are combined into
one. In this way, the controller ensures that both the second-order mean wave drift, wind, current forces
and the first-order roll moment are counteracted.
5.4.3. Azimuth angle

Since all thrusters pairs are aligned in the sway direction to achieve maximum roll reduction, compen-
sation of environmental forces in the surge direction is not possible. In beam waves, the wave forces
acting on the vessel in the surge direction are negligible. However, the wind and current can induce sig-
nificant forces in the surge direction, since these forces can differ from the wave direction. To increase
the station keeping capability of the vessel in these conditions, the azimuth angle of the thrusters used
for surge motion control are rotated by an angle as can be seen in Figure 5.2. There are two different
methods to determine the azimuth angle:
• Static method, in which fixed azimuth angles are used.
• Dynamic method, in which the azimuth angles are dependent on a particular variable.
First, the optimum static azimuth angles resulting in maximum roll reduction are determined by con-
ducting systematic simulations. During these simulations, the environmental forces induced by wind
5.4. Controllers 53
and current are varied between 0∘ and 180∘ with a constant interval of 30∘ . The wave forces are al-
ways beam on. The azimuth angles are varied in the interval [0∘ , 20∘ ] for thruster T3 and [0∘ , −30∘ ] for
thruster T5. The optimum azimuth angles per direction are shown in Figure 5.10.
20
T3
T5
Azimuth angle [deg]
15
10
0
0 20 40 60 80 100 120 140 160 180
[deg]
Figure 5.10: Optimum azimuth angle resulting in maximum roll reduction ( . m and . s)
From Figure 5.10 can be observed that the optimum azimuth angles vary significantly per environ-
mental direction. Next to this, one can see that the magnitude of the azimuth angle of T3 is bigger
compared to T5 in the first quadrant (0∘ − 90∘ ) and vice versa in the second quadrant (120∘ − 180∘ ).
This is the result of the fact that thruster T3 is used to counteract positive environmental forces and T5
is used to counteract negative environmental forces in the surge direction. From both observations is
concluded that a static method is not suitable to determine the azimuth angles of the thruster used for
surge motion control. A static method, where the azimuth angles are fixed to a specific value, does
not yield maximum roll reduction and more importantly, is not able to adapt to changing environmental
conditions.
The dynamic method is based on the idea that the azimuth angles need to be controllable in order
to adjust to changing environmental conditions. It was observed that there exist a relation between the
DP controller mean thrust command in the surge direction and the RMS of the roll angle time trace,
see Figure 5.11.
Roll angle RMS [deg] Mean T x command [kN]

1.54 50
1.52 45
-5 -5
1.5
40
Thruster angle T3 [deg]
Thruster angle T3 [deg]
-10 1.48 -10

35
1.46
30
-15 1.44 -15
25
1.42
-20 -20 20
1.4
15
1.38
-25 -25
1.36 10
-30 1.34 -30 5

2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20
Thruster angle T5 [deg] Thruster angle T5 [deg]
Figure 5.11: Azimuth angles for T3 and T5 vs. roll angle RMS and mean command ( . m and . s, ∘)
From Figure 5.11 can be observed that there exist a relation between the roll angle RMS and the mean
𝑇 command of the DP controller. When the mean 𝑇 command is minimum, the roll angle RMS is
also at a low value. Therefore, it is decided to minimize the DP controller 𝑇 command by using a
proportional controller:
𝜏 =𝐾 𝑒 (5.9)
where 𝑒 is defined by:

𝑒 =0−𝑇 (5.10)
By using a proportional azimuth angle controller, the controller behaves similar to a mechanical spring,
adjusting the azimuth angle when the environmental load increases. There is no mean offset or damp-
ing term that has to be included in this controller, since the azimuth angle controller is intended purely
proportional. Therefore, no 𝐾 and 𝐾 terms are included. To verify the behavior of the controller, the
average azimuth angles per environmental direction are visualized in Figure 5.12.
10
T3
T5
8
Azimuth angle [deg]

6
0
0 20 40 60 80 100 120 140 160 180
[deg]
Figure 5.12: Mean azimuth angles for T3 and T5 per environmental direction ( . m and . s)
From Figure 5.12 can be observed that the mean azimuth angles resulting from the proportional az-
imuth controller show behavior as expected. The azimuth angles of T3 are bigger compared to T5
for positive environmental forces in the surge direction. The opposite applies for T5 during negative
environmental forces in the surge direction. By implementing this type of azimuth angle controller, the
control system is able to counteract significant environmental loads in the surge direction and is able
to counteract the change in environmental loads.
The proposed azimuth angle controller is tuned by conducting systematic simulations per environmen-
tal direction and different values for 𝐾 . The results are given in Figure 5.13.
50
40
30
20
=0
10 =90
=180
0
0 0.5 1 1.5
Kp [-]
Figure 5.13: Roll reduction percentage for three environmental directions and tuning coefficients ( . m and . s)
From Figure 5.13 can be observed that the RMS roll reduction converges for increasing 𝐾 coefficient
values. The roll reduction percentages of every environmental direction are averaged and the 𝐾 co-
efficient which result in the highest averaged roll reduction percentage is selected. This resulted in a
value of 𝐾 = 0.9. A typical time trace of the surge, sway and yaw motions of the 3DP model is given
in Appendix D.
5.4.4. Shaft speed

The shaft speed controller is in place to keep the shaft speed equal to the required shaft speed as
determined by the high-level motion controllers. According Leonhard [18], a proportional-integral (PI)
controller is considered the most appropriate type for electrical motor control. Also Smogeli [32] pro-
poses a PI controller for shaft speed control. Therefore, a PI controller is implemented to control the
shaft speed:
𝜏 =𝑄 = 𝐾 𝑒 + 𝐾 ∫ 𝑒 (𝑡)𝑑𝑡 (5.11)
where 𝑄 and 𝑒 represent the commanded motor torque and shaft speed error calculated by:
𝑒 =𝜏−𝑛 (5.12)
where 𝜏 and 𝑛 represent the requested shaft speed by the high-level motion controller and the actual
shaft speed.
5.4. Controllers 55
A first estimate of the tuning parameters 𝐾 and 𝐾 is obtained by using the approach as formulated
by Smogeli [32]:
𝑘 𝑄
𝐾 =𝑘 (5.13)
𝑛
𝐾
𝐾 = (5.14)
𝑇
𝑇 =𝑘 𝑇 (5.15)
where 𝑘 , 𝑄 , 𝑛 , 𝑘 and 𝑇 represent a tuning parameter in the domain [5, 10], the propeller torque,
shaft speed at bollard pull condition, a tuning parameter ≥ 4 and a time constant defined by:
𝑇 =𝑇 +𝑇 (5.16)
where 𝑇 and 𝑇 represent the motor time constant and the shaft speed constant. The model is manu-
ally tuned to achieve desired controller behavior. The tuning coefficients are given in Appendix D. The
performance of the shaft speed controller for the thruster type TH1250 is visualized in Figure 5.14.
120
RPM Command
RPM Response
100
80
RPM [%]
60
40
20
0
640 645 650 655 660 665 670 675 680
Time [s]
Figure 5.14: Commanded RPM signal and actual RPM behavior of T3
As can be observed from Figure 5.14 the shaft speed controller performs as desired. No RPM over-
shoot and oscillatory behavior is observed.
5.5. Conclusions
The fourth sub-research object is defined as:
”Develop a control system model that enables combined DP station keeping and active roll reduc-
tion.”
In this chapter a new control system model is presented that combines DP station keeping and active
roll reduction. The newly proposed control model consist out of a high-level motion controller, which
is a combination of a proportional anti-roll controller and a conventional DP controller, and lower-level
shaft speed and azimuth angle controllers. Allocation matrices are used to transfer the controller com-
mands to the correct thrusters. The motion controller command is subsequently used as input variable
for the shaft speed controller in the dynamic thruster model. The azimuth angle of the thrusters used
for motion control in the surge direction is controlled by a proportional controller to ensure the ability
to counteract environmental forces in the surge direction. Tuning of the newly proposed controllers is
done by conducting systematic simulations. The conventional shaft speed and DP controller are tuned
by obtaining first estimates from literature.

• It is essential to form thruster pairs to achieve combined DP station keeping and active roll reduc-
tion.
• The optimal minimum RPM limit to achieve maximum roll reduction is 5% of the maximum thruster
RPM.
• A proportional controller is the most effective controller for active roll reduction.
• By implementing a proportional azimuth angle controller, the control system is able to counteract
environmental loads in the surge direction.
6
Results
In this chapter, results regarding the performance of the proposed control system model are presented
and compared to a conventional DP system model.
6.1. Introduction
The performance of the proposed control model is analyzed on the basis of multiple aspects. The
aspects that will be treated in this chapter are:
• Roll reduction
• DP footprint
• DP capability
• Thruster power consumption
• Workability
Each aspect will be treated in a dedicated section.
6.2. Roll reduction

In order to analyze the roll reduction obtained by the 3DP control system, systematic simulations in
irregular sea states have been carried out. First, the simulation is carried out without roll control ac-
tivated. Subsequently, active roll control is activated. A snapshot of the roll amplitude time traces is
given in Figure 6.1.
5
DP
4 3DP
2
Roll amplitude [deg]
-1
-2
-3
-4
-5
860 880 900 920 940 960 980 1000 1020 1040 1060
Time [s]
Figure 6.1: Snapshot of the roll angles time traces with and without active roll control ( . m and . s)
57
58 6. Results
As can be derived from Figure 6.1, the control system is able to actively reduce the roll angle ampli-
tudes of the vessel in irregular waves. To quantify the roll reduction, the RMS of both roll angle time
traces is calculated for a range of wave periods at a constant significant wave height and vice versa for
a range significant wave heights at a constant wave period. The results are observable in Figure 6.2.
Constant significant wave heigth (Hs =0.5m) Constant period (T p =8.5s)

1.5 6
DP DP
3DP 3DP
5
1 4
[deg]
[deg]
3
RMS
RMS
0.5 2
0 0
0 5 10 15 20 0 0.5 1 1.5 2 2.5 3
Tp [s] Hs [m]
Figure 6.2: Roll angle RMS as a function of wave period and significant wave height for both control models
From the first plot in Figure 6.2 can be observed that the roll reduction increases significantly in the
lower period range (𝑇 =4.5-7s), this can be explained by the fact that the thruster are less able to coun-
teract the short waves due to their response time, which is in the range of the wave period. In the higher
period range (𝑇 =8-20s) it can be observed that thrusters are able to significantly reduce the roll RMS.
The maximum reduction is achieved when 𝑇 = 8.5s, the natural frequency of the vessel.
The second plot in Figure 6.2 indicates that thruster induced roll reduction decreases when the signifi-
cant wave height increases. This is the result of the fact that the wave induced roll moment becomes
larger in comparison to the maximum thruster induced counteracting roll moment. Therefore, the rela-
tive contribution of the thrusters decreases when the significant wave height is increased.
To quickly obtain an estimate of the relative RMS roll reduction percentage at the roll natural frequency
for different sea states, the RRR results are fitted by an exponential model. The result is given in Figure
6.3.
100
Data
Fit
80
60
40
20
0
0 0.5 1 1.5
Hs [m]
Figure 6.3: Relative RMS roll angle percentage at the roll natural frequency versus significant wave height
The exponential fit model is given by:

.
𝑅𝑅𝑅 = 14.83 ⋅ 𝐻 (6.1)
As can be derived from Figure 6.3, the fit model somewhat underestimates the RMS roll reduction per-
centages during low significant wave heights and slightly overestimates the RRR percentage during
the higher wave heights. However, the model is deemed appropriate to obtain a quick estimate of the
6.3. DP footprint 59
expected roll reduction percentage for a particular significant wave height.
In order to analyze the roll reduction performance of the 3DP model in different conditions, system-
atic simulations with different values of 𝐻 and 𝑇 have been carried out. During these simulations, the
current velocity 𝑉 and wind velocity 𝑉 is kept constant at:
𝑉 = 3.32 m/s
𝑉 = 0.75 m/s
𝜇 = 90∘ (for both wind and current)
The resulting roll reduction during every sea state is visualized in the contour plot given in Figure 6.4.
20
5
18 0.2
16
14 25
0.
0.375 5
Tp [s]
0.2
0.
12
37
5
0.25
5
10 0.37
0.5 0.25
0.5
5
37
8 0.25
0.5
0.
0.375 0.375 0.25

0.25 0.25
6 0.25 0.25 0.25
0.125 0.125 0.125 0.125 0.125
0.125
4
0.2 0.4 0.6 0.8 1 1.2 1.4
Hs [m]
Figure 6.4: Contour plot of roll reduction for different combinations of and
As can be derived from Figure 6.4, the roll reduction is maximum for low sea states around the natural
frequency. The minimum roll reduction is found at the lower periods, since the thrusters are not able to
counteract the short wave periods. Again a trend is visible that the roll reduction decreases when the
significant wave height increases.
6.3. DP footprint
The DP footprint visualizes the maximum excursions of the vessel in the horizontal earth-fixed refer-
ence frame with respect to the DP setpoint. To assess the station keeping performance of the proposed
control model, the DP footprint of a conventional DP system is compared to the DP footprint of the 3DP
control system. The simulated cases are given in Table 6.1.
Case Draft Vw 𝜇 Vc 𝜇 Hs Tp 𝜇 Spectrum

[#] [m] [m/s] [deg] [m/s] [deg] [m] [s] [deg] [-]
1 4.7 0.39 90 0.75 90 0.5 8.5 90 JONSWAP
2 4.7 1.94 90 0.75 90 0.75 8.5 90 JONSWAP
3 4.7 3.32 90 0.75 90 1.0 8.5 90 JONSWAP
Table 6.1: Overview of simulation cases and corresponding environmental conditions
The wind and current velocity is chosen according the environmental conditions table given in IMCA.
The wave period is equal to the natural roll frequency of the vessel. The results are given in Figure 6.5.
From Figure 6.5 can be observed that the DP footprint of both control models increases during an in-
creasing significant wave height. This behavior is as expected, since the environmental loads increase.
Figure 6.5 also indicates that the DP footprint of the DP model is smaller compared to the DP footprint
of the 3DP model. This can be explained by the fact that the thrusters produce oscillating thrust levels
due to the combination of the anti-roll and DP controller. The maximum increase of the DP footprint is
1 meter in the sway direction. This is considered an acceptable magnitude.
60 6. Results
Hs =0.5m Hs =0.75m Hs =1.0m

0.5 0.5 0.5
Setpoint Setpoint Setpoint
0.4 DP 0.4 DP 0.4 DP
3DP 3DP 3DP
0.3 0.3 0.3
0.2 0.2 0.2
0.1 0.1 0.1
x [m]
x [m]
x [m]
0 0 0
-0.1 -0.1 -0.1
-0.2 -0.2 -0.2
-0.3 -0.3 -0.3
-0.4 -0.4 -0.4
-0.5 -0.5 -0.5

3 2 1 0 -1 -2 -3 3 2 1 0 -1 -2 -3 3 2 1 0 -1 -2 -3
y [m] y [m] y [m]
Figure 6.5: DP footprints results for case 1,2 and 3
Next to the station keeping performance of the vessel regarding its horizontal position, also the yaw
motions are of interest. The thrust variations will result in significant yaw moments, which the control
system has to counteract. A normal distribution has been used to visualized the yaw station keeping
performance of the vessel when both control models are applied, see Figure 6.6.
Hs =0.5m Hs =0.75m Hs =1.0m

20 20 20
DP DP DP
18 3DP 18 3DP 18 3DP
16 16 16
14 14 14
12 12 12
) [-]
) [-]
) [-]
10 10 10
6
6
f(
f(
f(
8 8 8
6 6 6
4 4 4
2 2 2
0 0 0
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
Yaw angle amplitude [deg] Yaw angle amplitude [deg] Yaw angle amplitude [deg]
Figure 6.6: Normal distribution probability density function of the yaw motion for case 1,2 and 3
From Figure 6.6 can be observed that DP model has the highest probability density at a yaw angle of
0∘ . This is as expected, since the DP model yaw angle setpoint is 0∘ . Figure 6.6 also shows that the
variance increases proportional to the significant wave height. This can be explained by the fact that
the vessel motions increase as a result of the increased environmental loads during higher significant
wave heights.
The 3DP model results show that the mean of the probability density function is slightly shifted to-
wards a yaw angle of 0.1∘ . This can be explained by the fact that the yaw moments in the 3DP model
are significantly higher compared to the 3DP model. The 𝐾 coefficient term of the DP controller in
the 3DP model should be increased to remove the yaw offset. However, for comparison reasons, it is
decided to use the same DP controller tuning variables in both the DP and the 3DP model.
The plots show that the probability density function remains more or less constant during increased
significant wave heights. This is the result of the fact that the yaw moments induced by the thrusters
are significantly bigger in magnitude compared to the yaw moments resulting from waves, current and
wind.
The yaw angle variance of the 3DP model is bigger compared to the DP model. This is explained
by the oscillating yaw moments induced by the thrusters during active roll reduction. The maximum
increase of the yaw angle is 0.7∘ , which is considered acceptable.
6.4. DP capability 61
6.4. DP capability
Also the station keeping capability of the vessel during active roll reduction in beam waves is of interest.
The station keeping capability is in this study defined by the limiting current speeds at which the vessel
is still able to maintain position. Maintaining position is defined according the positioning limits as de-
fined by DNV-GL [9] for DP capability level 3. The positioning limits herein are defined as a maximum
of 5 meter excursion from the DP setpoint and a maximum of +/-3∘ yaw angle excursion. The current
velocity is increased per environmental direction until the vessel exceeded the specified positioning
limits. The last current velocities at which the vessel was still able to remain within the positioning limits
are visualized in Figure 6.7
180
2.5
150
1.5
vc [m/s]
120
0.5
Hs=0.5m
0 90 H =0.75m
s
Hs=1.0m
60
30
0
[°]
Figure 6.7: Limiting current velocities per environmental direction for case 1,2 and 3, wave force is always beam on
From Figure 6.7 can be observed that the limiting current velocities per environmental direction de-
crease when the significant wave height increases. This is as expected since the environmental loads
on the vessel increase proportional to the wave height. The vessel is able to maintain position during
relative high current velocities of 2.2 m/s and 2.0 m/s in respectively head and following seas. In beam
seas, the station keeping capability is limited by a current velocity of 0.9 m/s for a 𝐻 =1m.
62 6. Results
6.5. Thruster power consumption

The 3DP control model is designed to obtain maximum roll reduction. The achieved roll reduction
comes at a cost. This cost can be expressed in terms of power consumed by the thrusters. The
thruster power consumption is calculated according Stapersma [34] by:
𝑃 = 2𝜋𝑄 𝑛 (6.2)
where 𝑄 and 𝑛 represent the propeller torque and the shaft speed.
In the DP control model, the shaft speed and propeller torque are not calculated. Therefore, a dif-
ferent approach is used to calculate the power consumption by the thrusters.
From Equation 6.2 can be derived that the thruster power consumption has a cubic relation with respect
to the shaft speed, since it is known that propeller torque scales ∼ with shaft speed, see Equation 3.13.
Using this relation, the power consumption of every thruster is calculated according:
.
𝑇
𝑃 =𝑃 ( ) (6.3)
𝑇
where 𝑃 , 𝑇, 𝑇 represent the maximum power consumption of the thruster, actual thrust produced
by the thruster and the maximum thrust magnitude of the thruster.
The sum of the power consumption of every thruster is calculated to obtain the total power consumption
of the control system over time. The mean of the total power consumption of both control models is
calculated and is expressed in percentage of the total installed thruster power. The results are given
in Figure 6.8
50
DP
45
3DP
40
Power consumption [%]
35
30
25
20
15
10
0
Hs =0.5m Hs =0.75m Hs =1.0m
Figure 6.8: Mean power consumption percentage of total installed DP power for case 1,2 and 3
Figure 6.8 clearly indicates that the power consumption of the 3DP control system is significantly higher
compared to the DP system. This is as expected, since the thrusters are constantly ramping up to
maximum thrust. It also logical that the power consumption of the thrusters is around 50% of the total
installed thruster power, since the thrusters are working 50% of their time at thrust levels near the max-
imum level to counteract the wave induced roll moment. The power consumption is around a factor 10
higher compared to the 3DP model.
From Figure 6.8 can also be observed that the power consumption of both control systems increases
during increasing sea state. This is as expected, since the environmental loads increase accordingly.
6.6. Workability 63
6.6. Workability
The final step is to analyze whether application of the 3DP control model results in increased vessel
workability. Workability is defined as the percentage of time the vessel is able to carry out operational
activities like crane operations, tool over-boardings and subsea power cable pull-ins. The ability to
conduct these operations is determined by vessel operability. The operability is defined by operating
criteria, see for example the NORDFORSK criteria in Journée [20].
In this research, the focus is on critical offshore operations, which are sensitive to roll motions. For
example, the over-boarding of a trencher for subsea power cable burial operations. It is assumed that
this type of operations have an operability limit of 1∘ roll RMS. To calculate the workability, the vessel
roll RMS is calculated for a range of sea states defined by different combinations of 𝐻 and 𝑇 . It is
assumed that the following conditions are present during each sea state:
𝑉 = 3.32 m/s
𝑉 = 0.75 m/s
𝜇 = 90∘ (for both wind and current)
The sea states during which the vessel roll RMS exceeds the defined operability limit are marked with
a ’x’. This procedure has been done for both the DP and the 3DP control system model. The results
are given in operability tables, see Table 6.2 and Table 6.3.
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6 x x
6.5 x x x x x x
7 x x x x x x x x x
7.5 x x x x x x x x x x x
8 x x x x x x x x x x x x
8.5 x x x x x x x x x x x x
10 x x x x x x x x x x
Table 6.2: Operability table for the DP control model
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6
6.5 x x x x
7 x x x x x x x
7.5 x x x x x x x x x
8 x x x x x x x x x
8.5 x x x x x x x x x x
10 x x x x x x x x
Table 6.3: Operability table for the 3DP control model

64 6. Results
To determine the workability increase, both the operability tables are applied to a yearly wave scatter
that represent the wave climate at the Borssele offshore wind farm. The wind farm location is visualized
in Figure 6.9
Figure 6.9: Location of the Borssele offshore wind farm (marked with x)
The wave scatter used is given in Table 6.4. The complete wave scatter is given in Appendix E.
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4 0.4 1.0 1.3 1.0 0.7 0.4 0.5 0.7 0.7 0.5 0.4 0.2 0.1 0.0
4.5 0.5 0.9 0.9 0.9 0.4 0.3 0.4 0.4 0.5 0.3 0.4 0.4 0.2 0.1
5 0.8 0.6 0.6 0.3 0.4 0.8 0.5 0.5 0.4 0.3 0.2 0.2 0.1 0.1
5.5 0.3 0.8 0.7 0.4 0.4 0.6 0.4 0.4 0.3 0.2 0.3 0.3 0.2 0.1
6 0.2 1.1 1.1 0.8 0.6 0.3 0.3 0.3 0.4 0.2 0.2 0.1 0.2 0.0
6.5 0.3 1.3 1.3 1.3 1.1 0.7 0.2 0.2 0.3 0.1 0.1 0.2 0.2 0.2
7 0.2 0.8 1.1 0.9 1.0 0.7 0.2 0.1 0.1 0.0 0.1 0.1 0.0 0.0
7.5 0.0 0.5 1.2 1.0 0.9 0.4 0.3 0.2 0.1 0.1 0.0 0.0 0.1 0.0
8 0.0 0.2 0.2 0.5 0.4 0.3 0.2 0.1 0.1 0.1 0.0 0.0 0.0 0.0
8.5 0.0 0.1 0.2 0.2 0.2 0.1 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0
9 0.2 0.1 0.1 0.1 0.2 0.0 0.1 0.0 0.0 0.0 0.0 0.1 0.0 0.0
9.5 0.1 0.0 0.1 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0
10 0.1 0.1 0.0 0.0 0.0 0.0 0.1 0.0 0.1 0.1 0.1 0.0 0.0 0.1
Table 6.4: Yearly wave scatter of the Borssele field, values are in percentage. Light grey represent sea states at which the
operability limit is exceeded for the DP model. Dark gray represent the 3DP model
The wave scatter given in Table 6.4 shows the percentage of time a particular sea state occurs at the
specific offshore location. By removing the percentages of the sea states marked with an ’x’ in the
operability tables and adding up the percentages in the wave scatter, the workability of the operation
is calculated. This process is indicated by the light gray and dark gray colored cells in Table 6.4.
The increased vessel workability is calculated by:
𝑊 =𝑊 −𝑊 (6.4)
where 𝑊 and 𝑊 represent respectively the workability obtained when the 3DP model and the DP
model are applied. This corresponds by summing up the light gray colored cells in Table 6.4. The
calculation results are given in Table 6.5.
𝑊 [%] 89.9
𝑊 [%] 95.1
𝑊 [%] 5.2
Table 6.5: Yearly workability calculation results in beam waves
A workability increase of 5.2% on a yearly basis corresponds to roughly 19 days. So by applying the
6.6. Workability 65
3DP model, the amount of days the vessel is limited by operability criteria decreased by 19 days. This
corresponds to a relative decrease of 48% of the amount of days the vessel is not able to operate. One
should take into account that this workability increase is only valid for operations in beam waves.
To assess the workability increase when the vessel is subjected to waves coming from different envi-
ronmental conditions, simulations have been carried out for waves coming from 0∘ to 180∘ with intervals
of 15∘ . The resulting operability tables for each environmental direction are given in Appendix E. The
yearly workability of the vessel per environmental direction is given in Table 6.6. As can be derived from
𝜇 [∘ ] 0 15 30 45 60 75 90 105 120 135 150 165 180

𝑊 [%] 100 100 100 99.7 95.5 92.3 89.9 91.3 94.6 98.9 100 100 100
𝑊 [%] 100 100 100 100 98.7 97 95.1 96 98.5 99.8 100 100 100
Table 6.6: Yearly workability per environmental direction for the DP and 3DP model for the Borssele field
Table 6.6, the workability of the vessel is only increased for environmental directions between 45∘ -135∘ .
This is due to the fact that the roll response of the vessel during the other environmental conditions
is much smaller in magnitude and is therefore not exceeding the operability limits. The workability in-
crease per direction is visualized in Figure 6.10.
180
6
150
5
Workability increase [%]
120
3
0 90
60
30
0
[°]
Figure 6.10: Yearly workability increase per environmental direction for the Borssele field when the 3DP model is engaged
As can be derived from Figure 6.10, the workability increase is maximum for beam waves and becomes
less for other environmental directions. This can be explained by the fact that during beam waves, the
vessel roll response is maximum. Application of the 3DP control system is therefore most effective
during this condition.
It is assumed that the probability of occurrence of each environmental direction during critical opera-
tions is equal. The yearly workability increase can therefore be calculated by averaging the workability
increase per direction. This results in an average workability increase of 2.0%. This corresponds to a
week of increased workability on a yearly basis for the Borssele field. One should take into account
that this estimate could change when a different wave scatter is used. For example, the Borssele wind
farm is close to shore, so shorter wave periods are present. At a more offshore location, longer wave
periods are more likely to occur, and therefore a higher workability increase can be expected.
7
Conclusions and Recommendations
In this chapter the conclusions are drawn based on the results of the conducted research. Next to this,
recommendations regarding further research are given.
7.1. Conclusions
The fifth sub-research objective is defined as:
”Identify and quantify the possibility and magnitude of thruster induced roll reduction and its effect
on the station keeping performance of the ship.”
From the results presented in Chapter 6 is concluded that:
• It is possible to actively reduce the roll motion amplitude of the Ndurance by active control of the
thrusters.
• The most significant roll reduction is obtained during sea states defined by significant wave
heights in the range of 0.2 to 0.7 meter and wave peak periods around the natural roll period
of the vessel. During these conditions, the achieved roll reduction is around 0.5∘ roll RMS.
• It is possible to combine active roll reduction and DP station keeping.
• The maximum DP footprint increase is 1 meter in the sway direction and 0.7∘ for the yaw motion
of the vessel when the 3DP model is activated.
• In beam seas, the station keeping capability during active roll reduction is limited by a current
velocity of 0.9 m/s during 𝐻 =1m.
• The thruster power consumption increases with a factor 10 when the 3DP control model is en-
gaged.
• The yearly workability of the Ndurance for the Borssele field increases with 5.2% in beam waves.
• The yearly workability increase for the Borssele field averaged over all environmental directions
is 2.0%.
The main research objective is to evaluate the possibility of reducing the ship roll motion during DP
operations by active control of the thrusters. From the conducted research is concluded that this is
possible.
67
68 7. Conclusions and Recommendations
7.2. Recommendations
The sixth and final sub-research objective is defined as:
”Give recommendations regarding possible thruster alternatives, the ship lay-out and DP system con-
figuration that influence the possibility and magnitude of actively controlled roll reduction.”
The recommendations based on the conducted research given below.

• Wear and tear of the thrusters
Probably the most import recommendation is to investigate the wear and tear experienced by the
thruster system during active roll reduction mode. This study confirms that thruster induced roll
reduction is possible. However, no attention has been paid to the wear and tear of the whole
thruster drive train system. Since the system will operate in quite an unorthodox way, it is ex-
pected that the wear and tear of especially the gearbox and transmission will be significant. The
increased maintenance cost, next to the increased fuel consumption, has to be taken into account
when one is looking into the performance of the 3DP control system from a financial point of view.
• Thruster alternatives
Alternative methods that can be used to achieve actively controlled roll reduction during DP op-
erations are active anti-roll tanks as investigated by Habing [12] and active fins as proposed by
Ikeda [15]. These systems however require vessel modifications and the installation of additional
equipment and power. The limited amount of increased operability is probably not worth such an
investment, as stated by Habing [12].
Passive roll reduction can be achieved by fitting bilge keels at the bilge of the Ndurance. This
option is also investigated by using a empirical bilge keel damping model, see Appendix F. From
this research is concluded that fitting bilge keels at the Ndurance results in 4 times less roll re-
duction compared to the 3DP model during low to moderate sea states with wave periods around
the roll natural period of the vessel. So fitting bilge keels at the Ndurance is probably not worth
the investment.
• Ship lay out
The achievable thruster induced roll reduction is probably more significant for ship-shaped ves-
sels than for barge-shaped vessels. This is due to the fact that barge-shaped vessel experience
a bigger roll moment in waves compared to a ship shaped vessel. The ratio between the counter-
acting thruster induced roll moment and the wave roll moment becomes smaller for ship-shaped
vessels. Therefore, the thrusters are able to more significantly reduce the roll motion. To achieve
optimum thruster induced roll reduction it is important that the thrusters are located as far away
as possible from the vessel CoG to obtain a large roll moment arm.
• Model tests
The performance of the constructed 3DP control model has only been investigated by conduct-
ing mathematical simulations. It is recommended to validate the implementation of a quadratic
viscous damping term in the time domain model by conduction systematic roll tests on model
scale for different wave heights and wave periods. On the basis of these results a new method or
model can be constructed that more accurately describes the viscous damping term in the time
domain.
• Operational use
The use of the 3DP control model is probably only beneficial for critical moments during offshore
operations. For example, when the vessel needs to conduct an operation that is at the roll oper-
ability limit and the weather forecast shows worsening weather conditions for the coming weeks.
The RRR fit model can be used to asses the effectiveness of engaging the 3DP system. When
considered effective, the 3DP system can enable the crew to carry out the operation, instead of
leaving the field and start the waiting on weather for multiple weeks. In this way, the 3DP control
model is more used as a back-up instrument instead of an operational mode that is engaged for
longer periods of time.
7.2. Recommendations 69
• New innovations
New innovations like lighter composite propellers and improved electrical drives are currently
introduced in the maritime industry. These innovations can increase the feasibility and perfor-
mance of thruster induced active roll reduction. Also new innovations can be thought off. For
example, one could design a system where the thruster pair ramping up could benefit from the
wasted energy of the thruster pair that is braking or ramping down at the same moment. Maybe
an intelligent energy conversion or transfer system can be developed and applied.
• Modelling approach
A lot of students use software packages like OrcaFlex or AQWA when one wants to conduct time
domain simulations. I would encourage everybody that is planning to do so during their graduation
to try to develop a time domain model by yourself. I believe it strongly increases the knowledge
regarding time domain simulations, vessel behavior and dynamics. Next to this, it enables you
in a lot of cases to conduct simulations much more faster and to post-process the results more
quickly.
A
Appendix
A.1. Vessel particulars
A photo of the Ndurance in the field at the Egmond aan Zee wind farm is given in Figure A.1.
Figure A.1: Ndurance in the field at the Egmond aan Zee wind farm
The vessel particulars of the cable-lay vessel Ndurance are given in Table A.1.
Vessel Ndurance
L [m] 99
B [m] 30
D [m] 4.7
∇ [m ] 11400
𝐼 [kgm ] 1.24E+09
𝐼 [kgm ] 7.33E+09
𝐼 [kgm ] 7.93E+09
𝐶 [N/m] 2.77E+07
𝐶 [Nm/rad] 1.25E+09
𝐶 [Nm/rad] 1.95E+10
𝐶 [Nm/rad] 7.74E+07
KG [m] 4.29
𝐴 [m ] 430
𝐴 [m ] 807
Table A.1: Vessel particulars Ndurance
71
72 A. Appendix
A.2. Ramp function

A ramp function is used when environmental loads are applied to the vessel time domain model. By
doing so, the transient response of vessel time domain model is decreased. Especially for calculation
of the motion RAO’s this is important to avoid, since one does not want a motion offset in the results.
The ramp function is chosen similar as in OrcaFlex and is calculated by:
𝑓 = 𝑢 (6𝑢 − 15𝑢 + 10) (A.1)
where 𝑢 is a parameter in the range [0, 1]. The length of the ramp function can be selected by adjusting
the slope of 𝑢:
1
Δ𝑢 = (A.2)
𝑡
where 𝑡 is the desired endtime of the ramp function. A visualization of the ramp function for an
endtime of 100 seconds is given in Figure A.2.
0.9
0.8
0.7
0.6
f ramp [-]
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60 70 80 90 100
Time [s]
Figure A.2: Ramp function for =100s

B
Appendix
B.1. Thruster and electrical drive data
The thruster parameters are given in Table B.1.
Thruster type TH1250 TH1500

D [m] 1.75 2.00
P/D [-] 1 1.2
Z [-] 4 4
EAR [-] 0.7 0.7
Propeller type Kaplan Kaplan
Nozzle 19A 19A
𝑘 [-] 3.409 6.207
Power [kW] 1000 1250
Thrust in bollard pull [kN] 177 222
𝐼 [kgm ] 10.2 9.2
𝐼 as fraction of 𝐼 [-] 0.385 0.332
𝐼 [kgm ] 0.6 0.6
𝑄 [Nm] 320 270
𝐾 [Nms/rad] 0.4 0.32
Roll moment arm [m] 10.2 7.7
Table B.1: Thruster TH1250 and TH1500 particulars
The electrical drive particulars are given in Table B.2.
Electrical drive AEM WH450-L6 Hyundai HRN3-48E

Power [kW] 1000 1250
Torque [Nm] 8000 6745
Rel. break down torque [-] 2.6 2.5
Nom. RPM [1/m] 1200 1800
Frequency [Hz] 60 60
𝑇 [s] 0.1 0.1
𝐼 [kgm ] 65.0 33.8
Table B.2: Electrical drive particulars
B.2. Four-quadrant Fourier coefficients

The four-quadrant coefficients 𝐴 , 𝐴 , 𝐵 and 𝐵 are given in Table B.3. The application of the coeffi-
cients resulted in a erroneous estimation of the thruster coefficients. It was concluded that Oosterveld
reported some mistaken coefficients in his report. After some trial-and-error procedures the correct
values for the Ka 4-70 propeller with P/D=1.0 and P/D=1.2 in the 19A nozzle were obtained. The cor-
rected values are given in Table B.3 and Table B.4.
73
74 B. Appendix
Order CT CQ
k A(k) B(k) A(k) B(k)
0 -1.09850E-01 0.0000E+00 3.15890E-02 0.00000E+00
1 1.40640E-01 -1.0583E+00 2.44060E-01 -1.17170E+00
2 1.57850E-01 4.7284E-02 -7.38800E-03 5.11550E-02
3 4.55440E-02 1.3126E-01 2.82600E-02 8.90690E-02
4 5.16390E-03 -7.7539E-03 -5.59590E-03 6.56700E-03
5 -2.55600E-03 9.3507E-02 2.65580E-04 1.42040E-01
6 -6.05020E-03 9.2520E-03 1.13680E-02 7.70520E-03
7 6.73680E-03 -1.4828E-02 -4.74010E-02 -3.60910E-02
8 6.85710E-03 -9.6554E-03 -6.56860E-03 4.20360E-03
9 4.72450E-03 9.6216E-03 -7.49900E-03 2.11390E-03
10 2.35910E-03 -7.5453E-04 1.28730E-03* 1.30950E-02
11 8.79120E-03 2.4453E-03 4.65020E-03 3.09610E-02
12 1.19680E-03 -8.7981E-03 -4.66760E-03 -9.94590E-03
13 8.38080E-03 1.8184E-03 3.34380E-03 1.79210E-02
14 -8.20980E-04 -2.0077E-03 2.20460E-03 -8.19170E-03
15 2.73710E-03 -3.3070E-03 7.00340E-03 -7.84280E-04
16 -2.61210E-04 -7.9201E-04 3.91470E-03* 7.26610E-03
17 1.91330E-03 -3.6311E-04 7.37190E-03 -4.73160E-03
18 3.22900E-04 -1.9377E-03* 9.40830E-04 -2.57310E-03
19 1.52230E-03 -1.2135E-03 6.05600E-03 1.11360E-03
20 -1.01510E-03 -3.1678E-04 4.23900E-04 -1.54700E-03
Table B.3: Corrected (marked with an *) coefficients as reported by Oosterveld for the Ka 4-70 (P/D=1.0) propeller in a 19A
nozzle [22]
Order CT CQ
k A(k) B(k) A(k) B(k)
0 -9.08880E-02 0.0000E+00 4.38000E-02 0.0000E+00
1 1.79590E-01 -1.0260E+00 3.52990E-01 -1.2949E+00
2 1.49560E-01 6.1459E-02 -1.09170E-02 5.9030E-02
3 6.56750E-02 1.3715E-01 4.70620E-02 9.3540E-02
4 5.21070E-03 -1.7280E-02 -1.07790E-02 -6.1148E-03
5 -6.82320E-03 9.6579E-02 -1.01930E-02 1.6121E-01
6 -6.28960E-03 5.8809E-03 -8.88240E-04 1.4624E-02
7 1.81780E-02 -2.2587E-02 -3.78930E-02 -5.3549E-02
8 6.06940E-03 -1.4819E-02 -7.03460E-03 -3.1589E-03
9 6.19420E-03 1.0398E-02* -8.01300E-03 1.4382E-02
10 2.64820E-03 -2.9324E-03 7.26220E-03 9.9836E-03
11 1.21370E-02 4.0913E-03 -5.43900E-03 3.8781E-02
12 -3.57050E-03 -4.4436E-03 -2.00380E-03 -4.6749E-03
13 3.29850E-03 -1.2190E-03 3.92810E-03 1.4944E-02
14 -8.86520E-04 -2.2551E-03 -6.52560E-04 -6.3253E-03
15 6.98070E-03 -3.2272E-03 1.54140E-02 2.2275E-03
16 -1.75600E-04 1.7533E-03 3.03560E-03 7.1826E-03
17 2.16430E-03 1.4875E-03 5.90730E-03 1.0229E-03
18 3.53620E-04 4.5353E-05 4.14330E-03 -5.9201E-03
19 2.57720E-03 -8.8702E-04 4.61020E-03 -1.4814E-03
20 -1.82790E-03 -9.4609E-04 -5.74230E-04 -4.3092E-03
Table B.4: Corrected (marked with an *) coefficients as reported by Oosterveld for the Ka 4-70 (P/D=1.2) propeller in a 19A
nozzle [22]
B.3. Thrust deduction factor 75
B.3. Thrust deduction factor

The thrust deduction factor due to thruster-hull interaction is based on the work of Wichers [35]. An
azimuth angle of 90∘ has been assumed. A factor of 0.04 is used, see Figure B.1.
Figure B.1: Thrust deduction percentage due to thruster-hull interaction as a factor of azimuth angle according Wichers [35]
The thruster moment arms in the horizontal plane with respect to the CoG are given in Table B.5
Thruster x [m] y [m]

T1 45 0
T2 36.5 7.5
T3 36.5 -7.5
T4 -45.5 7.5
T5 -45.5 -7.5
Table B.5: Thruster moment arms with respect to vessel CoG

C
Appendix
C.1. Wind and current coefficients
The wind coeffcients are given in Figure C.1
0.5
C coefficient [-]
-0.5 Cx
C
y
Cz
-1
0 50 100 150 200 250 300 350
[deg]
Figure C.1: Ndurance wind coefficients per environmental direction
The current coeffcients are given in Figure C.2
1
C
x
C
y
0.5 C
z
C coefficient [-]
-0.5
-1
0 50 100 150 200 250 300 350
[deg]
Figure C.2: Ndurance current coefficients per environmental direction
C.2. DP model tuning parameters

The tuning parameters used in the controller and the Kalman filter implemented in the DP model system
for comparison with the AnySim simulations are given in Table C.1.
Motion Surge Sway Yaw

𝐾 212.2 235.8 4932.21
𝐾 2309.35 2776.35 34781.5
Q 1.00E-09 2.00E-09 1.00E-09
R 1.30E-04 6.50E-05 1.30E-04
Table C.1: Controller and Kalman filter tuning parameters
The integral tuning coefficient 𝐾 is omitted, since this is also the case in the AnySim simulations.
77
78 C. Appendix
C.3. Forbidden zones

The performance of the allocation algorithm is checked by subjecting the vessel to environmental loads
coming from 0∘ to 360∘ . The calculated thruster angles by the allocation algorithm to counteract the
environmental loads are visualized in Figure C.3.
Figure C.3: Thrust magnitude and azimuth angles for different environmental directions, the forbidden zones are indicated by
the dotted triangles
As can be observed from Figure C.3, the allocation algorithm does not result in azimuth angles within the
forbidden zones. From this result is concluded that the allocation algorithm is performing as intended.
D
Appendix
D.1. Tuning coefficients of 3DP model
The tuning coefficients of the DP controller in the 3DP control system model are given in Table D.1.
DP Controller Surge Sway Yaw

𝐾 106.1 117.9 4932.2
𝐾 0.25 0.30 0.01
𝐾 1154.7 1388.2 34781.5
Q 1.00E-09 2.00E-09 1.00E-09
R 1.30E-04 6.50E-05 1.30E-04
Table D.1: Tuning coefficients of the DP controller and Kalman filter in the 3DP model
The controller coefficients as used in the 3DP control model for the roll controller, azimuth controller
and shaft speed controllers are given in Table D.2.
Roll controller
𝐾 2.50E+06
Azimuth controller
𝐾 0.9
Shaft speed controller (TH1250)
𝐾 33.3
𝐾 16.65
Shaft speed controller (TH1500)
𝐾 19.9
𝐾 9.95
Table D.2: Tuning coefficients of the roll controller, azimuth controller and shaft speed controllers as implemented in the 3DP
control model
79
80 D. Appendix
D.2. Motion time traces

A typical simulation time trace of the horizontal motions surge, sway and when both the DP and the
3DP model are applied are given in Figure D.1.
Figure D.1: Time traces of horizontal motions when both DP and 3DP model are applied ( =0.75m, =8s
E
Appendix
E.1. Yearly wave scatter Borssele
In Table E.1, the complete yearly wave scatter of the Borssele field is given.
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6
1 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1.5 2.3 2.6 1.3 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 0.6 3.6 3.5 1.4 1.1 0.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2.5 1.0 0.8 2.6 3.5 1.7 0.4 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3 1.3 1.8 0.9 1.2 1.5 1.6 0.5 0.3 0.1 0.0 0.0 0.0 0.0 0.0 0.0
3.5 0.5 1.7 0.8 0.7 0.7 0.8 0.8 0.6 0.3 0.2 0.1 0.0 0.0 0.0 0.0
4 0.4 1.0 1.3 1.0 0.7 0.4 0.5 0.7 0.7 0.5 0.4 0.2 0.1 0.0 0.0
4.5 0.5 0.9 0.9 0.9 0.4 0.3 0.4 0.4 0.5 0.3 0.4 0.4 0.2 0.1 0.0
5 0.8 0.6 0.6 0.3 0.4 0.8 0.5 0.5 0.4 0.3 0.2 0.2 0.1 0.1 0.1
5.5 0.3 0.8 0.7 0.4 0.4 0.6 0.4 0.4 0.3 0.2 0.3 0.3 0.2 0.1 0.1
6 0.2 1.1 1.1 0.8 0.6 0.3 0.3 0.3 0.4 0.2 0.2 0.1 0.2 0.0 0.0
6.5 0.3 1.3 1.3 1.3 1.1 0.7 0.2 0.2 0.3 0.1 0.1 0.2 0.2 0.2 0.0
7 0.2 0.8 1.1 0.9 1.0 0.7 0.2 0.1 0.1 0.0 0.1 0.1 0.0 0.0 0.1
7.5 0.0 0.5 1.2 1.0 0.9 0.4 0.3 0.2 0.1 0.1 0.0 0.0 0.1 0.0 0.1
8 0.0 0.2 0.2 0.5 0.4 0.3 0.2 0.1 0.1 0.1 0.0 0.0 0.0 0.0 0.0
8.5 0.0 0.1 0.2 0.2 0.2 0.1 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0
9 0.2 0.1 0.1 0.1 0.2 0.0 0.1 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0
9.5 0.1 0.0 0.1 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.0
10 0.1 0.1 0.0 0.0 0.0 0.0 0.1 0.0 0.1 0.1 0.1 0.0 0.0 0.1 0.0
10.5 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Table E.1: Yearly wave scatter of the Borssele field, values are in percentage
81
82 E. Appendix
E.2. Operability tables

In this section the operability tables for 𝜇 = 45∘ to 𝜇 = 135∘ are given with intervals of 15∘ .
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6
6.5
7
7.5
8 x x x
8.5 x x x x
9 x x x x
9.5 x x
10 x
Table E.2: Operability table ∘. Light grey represents sea states at which the operability limit is exceeded for the DP
model. Dark gray represents the 3DP model
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6
6.5 x
7 x x x x x
7.5 x x x x x x x x
10 x x x x x x x x
E.2. Operability tables 83
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6
6.5 x x x x x
7 x x x x x x x x
8 x x x x x x x x x x x
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6 x x
6.5 x x x x x x
7 x x x x x x x x x
8.5 x x x x x x x x x x x x
84 E. Appendix
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6
6.5 x x
7 x x x x x x
10 x x x x x x x x
Tp [s]/Hs [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
4
4.5
5
5.5
6
6.5
7
7.5 x x x
8 x x x x x x
8.5 x x x x x x x
9 x x x x x x
9.5 x x x x x
10 x x x x
F
Appendix
A commonly used method to passively damp the roll motion of a vessel is fitting bilge keels at the bilge
radius of the vessel. To be able to quantify and compare the roll reduction resulting from fitting bilge
keels at the Ndurance, a short study is carried out. It is assumed that the bilge keels are fitted on both
sides of the vessel over a length of 54 meters, where 𝐶 = 0.99. The bilge keel width is assumed 0.3
meters. A visualization of the bilge keel configuration at the starboard side of the vessel is given in
Figure F.1.
z
0.3m
Figure F.1: Schematic visualization of bilge keel fitted at section
F.1. Empirical bilge keel damping model

According Chakrabarti [5], the viscous bilge keel damping consist out of two components. These com-
ponents are damping terms due to the normal and hull pressure induced by the presence of the bilge
keels:
𝐵 =𝐵 +𝐵 (F.1)
Bilge keel normal pressure damping
The normal component per unit length is defined by:
8
𝐵 = 𝜌𝑟 𝑏 𝜔𝑅 𝑓 𝐶 (F.2)
3𝜋

𝑏 = bilge keel width
The equivalent drag coefficient 𝐶 is defined by:
𝑏
𝐶 = 22.5 + 2.4 (F.3)
𝜋𝑟 𝑅 𝑓
85
86 F. Appendix
where the factor 𝑓 is included to account for the flow speed at the bilge keel. It is defined by:
𝑓 = 1 + 0.3exp(−160(1 − 𝜎)) (F.4)
where 𝜎 represents the area coefficient of a cross-section of the hull. The mean distance from the roll
axis to the bilge keel 𝑟 is defined by:
.
𝑟 = 𝐷((𝐻 − 0.293𝑅 /𝐷) + (1 − 𝑂𝐺/𝐷 − 0.293𝑅 /𝐷) ) (F.5)
where
𝐷 = vessel draft
𝑅 = bilge radius
𝐻 = half the beam-draft ratio
𝑂𝐺 = vertical distance from SWL to COG
Bilge keel hull pressure damping

The pressure damping component due to the hull surface is defined by:
4 𝑏
𝐵 = 𝜌𝑟 𝐷 𝜔𝑅 𝑓 ( − ( − 22.5 − 1.2)𝐴 + 1.2𝐵 ) (F.6)
3𝜋 𝜋𝑟 𝑓𝑅
where 𝐴 and 𝐵 are defined by:
𝐴 = (𝑚 + 𝑚 )𝑚 − 𝑚 (F.7)
𝑚 (1 − 𝑚 ) (2𝑚 − 𝑚 )
𝐵 = + + (𝑚 𝑚 + 𝑚 𝑚 )𝑚 (F.8)
3(𝐻 − 0.215𝑚 ) 6(1 − 0.215𝑚 )
The 𝑚-coeffients are given by:
𝑚 = 𝑅 /𝐷 (F.9)
𝑚 = 𝑂𝐺/𝐷 (F.10)
𝑚3 = 1 − 𝑚 − 𝑚 (F.11)
𝑚 =𝐻 −𝑚 (F.12)
0.414𝐻 + 0.0651𝑚 − (0.382𝐻 + 0.0106)𝑚
𝑚 = (F.13)
(𝐻 − 0.215𝑚 )(1 − 0.215𝑚 )
0.414𝐻 + 0.0651𝑚 − (0.382 + 0.0106𝐻 )𝑚
𝑚 = (F.14)
(𝐻 − 0.215𝑚 )(1 − 0.215𝑚 )
𝑚 = 𝑆 /𝐷 − 0.25𝜋𝑚 (F.15)
𝑚 = 𝑚 + 0.414𝑚 (F.16)
where 𝑆 is defined by:

𝑆 = 0.3𝜋𝑓𝑟 𝑅 + 1.95𝑏 (F.17)
F.2. Bilge keel model implementation

The damping coefficients are calculated by using the roll motion amplitude estimates of the frequency
domain model without the implementation of a viscous damping term. Subsequently, the viscous damp-
ing coefficients are calculated and implemented in the frequency domain model:
[M + A(𝜔)]𝜂̈ + B(𝜔)𝜂̇ + 𝐵 𝜙̇ + C𝜂 = F(𝜔) cos(𝜔𝑡 + 𝜖) (F.18)
where 𝐵 is now defined by:

𝐵 =𝐵 +𝐵 +𝐵 (F.19)
where 𝐵 and 𝐵 represent the eddy making and skin friction damping respectively. A similar iterative
approach as described in Chapter 2 is used to obtain the converged viscous roll damping values.
F.2. Bilge keel model implementation 87
108 Eddy making damping 106 Skin friction damping


15 2.5
Iterations 2 Iterations
10
1.5
1
5
0.5
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
107 Bilge keel normal pressure damping 108 Bilge keel hull pressure damping

8 4
6 Iterations 3 Iterations
4 2
2 1
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
108 Converged total viscous damping Converged roll RAO
3 8
Roll RAO [deg/m]
6
2
1
2
0 0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Figure F.2: Results of iterative calculation of viscous damping RAO’s with bilge keel damping included and final roll motion
RAO in the frequency domain model ( ∘ , T=4.7m, WD=1000m)
The results of the iterative procedure are visualized in Figure F.2.
In the time domain model the viscous damping term is included by a quadratic viscous damping term.
This is considered a more physical approach, since it is known that drag forces have a quadratic re-
lation with velocity, see for example Chakrabarti [5]. The quadratic viscous damping term is defined
by:
𝐵 = 𝑐 |𝜙|̇ 𝜙̇ (F.20)
where 𝑐 is a fitting parameter. The value of this fitting parameter is obtained by fitting the time domain
roll motion RAO to the frequency domain roll motion RAO with a linearized viscous roll damping term
included. The time domain model with a viscous roll damping term included is now defined by:
[M + A ]𝜂̈ + K 𝜂̇ + 𝑐 |𝜙|̇ 𝜙̇ + C𝜂 = 𝐹(𝑡) (F.21)
By adjusting the value of 𝑐 , the time domain roll motion RAO is fitted to the frequency domain roll
motion RAO. The result is visualized in Figure F.3
The value of the tuning parameter is determined as 𝑐 = 3.5 ⋅ 10 . As can be derived from Figure F.3,
the time and frequency domain model results match.
88 F. Appendix
7
Time domain model
6
Roll RAO [deg/m]

4
0
0 0.5 1 1.5 2 2.5
Frequency [rad/s]
Figure F.3: Results of time domain model RAO fitting to frequency domain RAO with bilge keel damping included ( ∘,
T=4.7m, WD=1000m)
F.3. Bilge keel model results

The RMS of the roll motion amplitudes for different wave periods and wave heights are calculated when
the bilge keels are applied. The result is given in Figure F.4. To enable comparison with the 3DP model
and DP model, both results are also displayed in Figure F.4.
Constant significant wave heigth (Hs =0.5m) Constant period (T p =8.5s)

1.5 6
DP DP
3DP 3DP
Bilge keels 5 Bilge keels
1 4
[deg]
[deg]
3
RMS
RMS
0.5 2
0 0
0 5 10 15 20 0 0.5 1 1.5 2 2.5 3
Tp [s] Hs [m]
Figure F.4: Roll angle RMS as a function of wave period and significant wave height for both control models and when bilge
keels are applied
From the results visualized in Figure F.4 can be observed that the relative roll reduction, at the roll
natural frequency, realized by the 3DP is a factor 4 bigger compared to the roll reduction when the
bilge keels are applied. It is also concluded that roll reduction by the bilge keels is more significant for
wave heights greater than 1.5m. This can be explained by the increased roll velocities during these
relative high sea states.
From this study is concluded that fitting bilge keels on the Ndurance reduces the roll motion amplitudes
with a maximum of 10%, whereas the 3DP control model is able to realize a roll reduction of 40%. The
bilge keel induced relative roll reduction increases to 15% for sea states exceeding significant wave
heights of 1.5 meter.
List of Figures
1.1 Schematic visualization of the Ndurance . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2.1 Ship motions in 6 degrees of freedom [20] . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Panel model of the Ndurance in AQWA . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Ndurance Motion RAOs (𝜇 = 45∘ , T=4.7m, WD=1000m) . . . . . . . . . . . . . . . . . . 8
2.4 Strategy to obtain state-space representation . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 Calculated motion IRFs using the OrcaFlex and Journée models, 𝑇 =50s . . . . . . . . 11
2.6 Frequency-dependent added mass for six ship motions and calculated infinity frequency
added mass estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.7 Frequency domain motion IRFs and fitted transfer functions . . . . . . . . . . . . . . . . 14
2.8 Comparison of calculated RAOs by both frequency and time domain model (𝜇 = 45∘ ,
T=4.7m, WD=1000m) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.9 Schematic visualization of eddy making damping and skin friction damping . . . . . . . 17
2.10 Results of iterative calculation of viscous damping RAOs and final roll motion RAO in the
frequency domain model (𝜇 = 90∘ , T=4.7m, WD=1000m) . . . . . . . . . . . . . . . . . 19
2.11 Results of time domain model RAO fitting to frequency domain RAO (𝜇 = 90∘ , T=4.7m,
WD=1000m) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.12 Viscous roll damping trend visualization in both frequency and time domain model . . . 20
3.1 Schematic structure of thruster system lay-out . . . . . . . . . . . . . . . . . . . . . . . 23

3.2 Schematic visualization of oscillating thruster inflow velocities due to the rolling motion
of the vessel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.3 Four-quadrant diagram of Ka 4-70 thruster with different P/D ratios as found in Ooster-
veld [22] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 Model results for different advance velocities compared to manufacturer data for thruster
type 1250 and 1500 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.5 Transient response for thruster type 1250 and 1500 models . . . . . . . . . . . . . . . . 28
4.1 Definition of reference frames and environmental force direction . . . . . . . . . . . . . 31

4.2 Schematic visualization of the thruster system lay-out of the Ndurance . . . . . . . . . . 32
4.3 Schematic overview of DP system model . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.4 Block diagram of closed-loop PID controller . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.1 Thruster pair configuration 3DP control model for roll reduction . . . . . . . . . . . . . . 46
5.2 Thruster configuration 3DP control model for station keeping (sway, surge and yaw motion) 47
5.3 3DP control model structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.4 Selected search space for controller tuning coefficients in range [0.01 10000], 𝐻 = 0.5m
and 𝑇 = 8.5s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.5 Roll reduction for different tuning coefficients and maximum sensitivity of every tuning
coefficient (𝐻 = 0.5m and 𝑇 = 8.5s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.6 Roll reduction percentage for different 𝐾 values (𝑇 = 8.5s) . . . . . . . . . . . . . . . . 50
5.7 Response time of thruster TH1250 and TH1500 for different minimum limit RPMs . . . . 50
5.8 Roll reduction percentage for different 𝜏 as percentage of 𝜏 (𝑇 = 8.5s) . . . . . 51
5.9 3DP control system high-level motion control principle . . . . . . . . . . . . . . . . . . . 52
5.10 Optimum azimuth angle resulting in maximum roll reduction (𝐻 = 0.5m and 𝑇 = 8.5s) . 53
5.11 Azimuth angles for T3 and T5 vs. roll angle RMS and mean 𝑇 command (𝐻 = 0.5m
and 𝑇 = 8.5s, 𝜇 = 60∘ ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.12 Mean azimuth angles for T3 and T5 per environmental direction (𝐻 = 0.5m and 𝑇 = 8.5s) 54
89
90 List of Figures
5.13 Roll reduction percentage for three environmental directions and tuning coefficients 𝐾
(𝐻 = 0.5m and 𝑇 = 8.5s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.14 Commanded RPM signal and actual RPM behavior of T3 . . . . . . . . . . . . . . . . . 55
6.1 Snapshot of the roll angles time traces with and without active roll control (𝐻 = 0.5m
and 𝑇 = 8.5s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.2 Roll angle RMS as a function of wave period and significant wave height for both control
models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.3 Relative RMS roll angle percentage at the roll natural frequency versus significant wave
height . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.4 Contour plot of roll reduction for different combinations of 𝑇 and 𝐻 . . . . . . . . . . . 59
6.5 DP footprints results for case 1,2 and 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.6 Normal distribution probability density function of the yaw motion for case 1,2 and 3 . . 60
6.7 Limiting current velocities per environmental direction for case 1,2 and 3, wave force is
always beam on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.8 Mean power consumption percentage of total installed DP power for case 1,2 and 3 . . 62
6.9 Location of the Borssele offshore wind farm (marked with x) . . . . . . . . . . . . . . . . 64
6.10 Yearly workability increase per environmental direction for the Borssele field when the
3DP model is engaged . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
A.1 Ndurance in the field at the Egmond aan Zee wind farm . . . . . . . . . . . . . . . . . . 71
A.2 Ramp function for 𝑡 =100s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
B.1 Thrust deduction percentage due to thruster-hull interaction as a factor of azimuth angle
according Wichers [35] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
C.1 Ndurance wind coefficients per environmental direction . . . . . . . . . . . . . . . . . . 77

C.2 Ndurance current coefficients per environmental direction . . . . . . . . . . . . . . . . . 77
C.3 Thrust magnitude and azimuth angles for different environmental directions, the forbid-
den zones are indicated by the dotted triangles . . . . . . . . . . . . . . . . . . . . . . . 78
D.1 Time traces of horizontal motions when both DP and 3DP model are applied (𝐻 =0.75m,
𝑇 =8s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
F.1 Schematic visualization of bilge keel fitted at section . . . . . . . . . . . . . . . . . . . . 85

F.2 Results of iterative calculation of viscous damping RAO’s with bilge keel damping in-
cluded and final roll motion RAO in the frequency domain model (𝜇 = 90∘ , T=4.7m,
WD=1000m) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
F.3 Results of time domain model RAO fitting to frequency domain RAO with bilge keel damp-
ing included (𝜇 = 90∘ , T=4.7m, WD=1000m) . . . . . . . . . . . . . . . . . . . . . . . . . 88
F.4 Roll angle RMS as a function of wave period and significant wave height for both control
models and when bilge keels are applied . . . . . . . . . . . . . . . . . . . . . . . . . . 88
List of Tables
2.1 Description of ship motions and symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Convolution term properties and corresponding implications on the fitted transfer function
(TF) according Perez [25] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 RMAE values of indicated ship motions . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.1 Fitting coefficients according MacPherson [19] . . . . . . . . . . . . . . . . . . . . . . . 25

3.2 Four quadrant operation of a propeller [34] . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 RMAE values for thruster type TH1250 and TH1500 during different advance velocities 28
4.1 Thruster azimuth angle forbidden zones . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4.2 Overview of simulation cases and corresponding environmental conditions for model
verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.3 Comparison of simulation results for case 1 . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.6 RMAE values for model results compared to MARIN results . . . . . . . . . . . . . . . . 42
5.1 Differences between a conventional DP model and the 3DP model . . . . . . . . . . . . 45
6.1 Overview of simulation cases and corresponding environmental conditions . . . . . . . 59

6.2 Operability table for the DP control model . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.3 Operability table for the 3DP control model . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.4 Yearly wave scatter of the Borssele field, values are in percentage. Light grey repre-
sent sea states at which the operability limit is exceeded for the DP model. Dark gray
represent the 3DP model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.5 Yearly workability calculation results in beam waves . . . . . . . . . . . . . . . . . . . . 64
6.6 Yearly workability per environmental direction for the DP and 3DP model for the Borssele
field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
A.1 Vessel particulars Ndurance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
B.1 Thruster TH1250 and TH1500 particulars . . . . . . . . . . . . . . . . . . . . . . . . . . 73

B.2 Electrical drive particulars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
B.3 Corrected (marked with an *) coefficients as reported by Oosterveld for the Ka 4-70
(P/D=1.0) propeller in a 19A nozzle [22] . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
B.4 Corrected (marked with an *) coefficients as reported by Oosterveld for the Ka 4-70
(P/D=1.2) propeller in a 19A nozzle [22] . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
B.5 Thruster moment arms with respect to vessel CoG . . . . . . . . . . . . . . . . . . . . . 75
C.1 Controller and Kalman filter tuning parameters . . . . . . . . . . . . . . . . . . . . . . . 77
D.1 Tuning coefficients of the DP controller and Kalman filter in the 3DP model . . . . . . . 79
D.2 Tuning coefficients of the roll controller, azimuth controller and shaft speed controllers
as implemented in the 3DP control model . . . . . . . . . . . . . . . . . . . . . . . . . . 79
E.1 Yearly wave scatter of the Borssele field, values are in percentage . . . . . . . . . . . . 81
E.2 Operability table 𝜇 = 45∘ . Light grey represents sea states at which the operability limit
is exceeded for the DP model. Dark gray represents the 3DP model . . . . . . . . . . . 82
91
92 List of Tables
Bibliography
[1] Jose A. Armesto, Raul Guanche, Fernando del Jesus, Arantza Iturrioz, and Inigo J. Losada. Com-
parative analysis of the methods to compute the radiation term in cummins’ equation. Ocean
Engineering and Marine Energy, 1:377–393, 2015.
[2] J. Bargmeyer. Wetenswaardigheden over frequentieomvormers. Danfoss, Nordborg, 2016.
[3] R.H. Byrd et al. A trust region method based on interior point techniques for nonlinear program-
ming. Mathematical Programming, 89(1):149–185, 2000.
[4] O. Cadet. Introduction to kalman filter and its use in dynamic positioning systems. In MTS Dynamic
Positioning Conference September 16-17, 2003.
[5] Subrata Chakrabarti. Empirical calculation of roll damping for ships and barges. Ocean Engineer-
ing, 28:915–932, 2001.
[6] J.G. Cooke. Incorporating Thruster Dynamics in the Control of an Underwater Vehicle. MIT, 1989.
[7] W.E. Cummins. The impulse response function and ship motions. In Symposium Ship Theory
Hamburg, Germany, 1962.
[8] G. de Backer. Hydrodynamic design optimization of wave energy converters consisting of heaving
point absorbers. PhD thesis, Ghent University, 2009.
[9] DNV-GL. St-0111 assessment of station keeping capability of dynamic positioning vessels. DNV-
GL Standard, 2016.
[10] T. Duarte. SS Fitting - Theory and user manual. NWTC, 2012.
[11] T.I. Fossen. Handbook of Marine Craft Hydrodynamics and Motion Control. Wiley & Sons Ltd.,
Sussex, 2011.
[12] R. Habing. Roll damping to improve the operational limits of existing barge-shaped vessels. TU
Delft, Delft, 2016.
[13] H. Hatecke. Robust identification of parametric radiation force models via impulse response fitting.
Proceedings in Applied Mathematics and Mechanics, 15(1):593–594.
[14] Y. Himeno. Prediction of ship roll damping - state of the art. University of Michigan, 1981.
[15] Y. Ikeda, T. Fujiwara, and T.Katayama. Roll damping of a sharp-cornered barge and roll control
by a new-type stabilizer. In Proceedings of the Third International Offshore and Polar Engineering
Conference, Singapore, 1993.
[16] J.M.J. Journée. Hydromechanic coefficients for calculating time domain motions of cutter suction
dredgers by cummins equation, 2000.
[17] Philipp Koschorrek et al. Dynamic positioning with active roll reduction using voith schneider
propeller. In 10th IFAC Conference on Manoeuvring and Control of Marine Craft, 2015.
[18] W. Leonhard. Control of Electrical Drives 2nd edition. Springer-Verlag, Berlin, 1996.
[19] Donald M. MacPherson, Vincent R. Puleo, and Matthew B. Packard. Estimation of entrained water
added mass properties for vibration analysis. SNAME New England section, 2007.
[20] J.M.J. Journée & W.W. Massie. Offshore Hydromechanics. Delft University of Technology, Delft,
2001.
93
94 Bibliography
[21] J.N. Newman. Marine Hydrodynamics. MIT, Massachusetts, 1977.

[22] M.W.C. Oosterveld. Wake Adapted ducted propellers. PhD thesis, Delft University of Technology,
Delft, 1970.
[23] Orcina. Orcaflex documentation. https://www.orcina.com/SoftwareProducts/OrcaFlex/Documentation,

2018.
[24] T. Perez. Parametric time-domain models based on frequency domain data, 2007.
[25] T. Perez. Practical aspects of frequency-domain identification of dynamic models of marine struc-
tures from hydrodynamic data. Ocean Engineering, 38(1):426–435, 2011.
[26] Ivan Pinéda. The European offshore wind industry - key trends and statistics 2015. EWEA,
Brussel, 2015.
[27] P. Rici et al. Time-domain models and wave energy converters performance assessment. In 27th
ASME Conference on Ocean, Offshore and Arctic Engineering OMAE2008, 2008.
[28] Sigbørn Eng Rudaa. Use of tunnel thrusters and azimuthing thrusters for roll damping of ships,
2016.
[29] H.E. Saunders. Hydrodynamics in Ship Design. The Society of Naval Architects and Marine
Engineers, New York, 1957.
[30] H. Schwanecke. Gedanken zur frage der hydrodynamisch erregten schwingungen des propellers
und der wellenleitung. Jahrbuch der Schiffbauteschnischen Gesellschaft 57. Band, 1963.
[31] Jorrit Jan Serraris. Time domain analysis for dp simulations. In 28th ASME Conference on Ocean,
Offshore and Arctic Engineering OMAE2009, 2009.
[32] Øyvind Notland Smogeli. Control of marine propellers. PhD thesis, NTNU, Trondheim, 2006.
[33] Asgeir J. Sørensen and Jann Peter Strand. Positioning of small-waterplane-area marine construc-
tions with roll and pitch damping. Control Engineering Practice, 8:205–213, 2000.
[34] H. Klein Woud & D. Stapersma. Design of Propulsion and Electric Power Generation Systems.
IMAREST, London, 2002.
[35] Johan Wichers, Stephen Bultema, and Richard Matten. Hydrodynamic research on and optimizing
dynamic positioning system of a deep water drilling vessel. In Offshore Technology Conference,
1996.
[36] Shengwen Xu et al. Mitigating roll-pitch motion by a novel controller in dynamic positioning system
for marine vessels. Ships and offshore structures, 12:1136–1144, 2017.

Thesis RGdeJong

Uploaded by

Copyright:

Available Formats

Thesis RGdeJong

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Thesis RGdeJong

Uploaded by

Copyright:

Available Formats

What is the research aiming to evaluate?

What is the research aiming to evaluate?

What models have been constructed and applied in the research?

What models have been constructed and applied in the research?

3DP

A control system model for combined DP

Student number: 4195647

An electronic version of this thesis is available at http://repository.tudelft.nl/

DOF - Degrees of Freedom

IRF - Impulse Response Function

QTF - Quadratic Transfer Function

RAO - Response Amplitude Operator

RMS - Root Mean Square

RPM - Rounds Per Minute

RRR - Relative Roll Reduction

VFD - Variable Frequency Drive

𝜂 vessel motion vector [m] or [rad]

A infinite frequency added mass matrix [kg] or [kgm ]

𝐶 four-quadrant thrust coefficient [-]

5 3DP System Model 45

Figure 1.1: Schematic visualization of the Ndurance

1.2. Problem statement

1.3. Research objective

1.4. Research approach

1.5. Report structure

2.1.1. Ship motions

Motion Name Symbol Unit

2.2. Frequency domain model

[M + A(𝜔)]𝜂̈ + B(𝜔)𝜂̇ + C𝜂 = F(𝜔) cos(𝜔𝑡 + 𝜖) (2.1)

M = ship inertia matrix

Figure 2.2: Panel model of the Ndurance in AQWA

𝑆 (𝜔) = motion response spectrum

Surge RAO Sway RAO

Surge amplitude [m/m]

Sway amplitude [m/m]

Heave RAO Roll RAO

Roll amplitude [deg/m]

Yaw amplitude [deg/m]

Figure 2.3: Ndurance Motion RAOs ( ∘, T=4.7m, WD=1000m)

2.3. Time domain model

[M + A ]𝜂̈ + ∫ K(𝑡 − 𝜏)𝜂𝑑𝜏

A = infinite frequency added mass matrix

function matrix are defined by:

2.3.1. The convolution term

Time domain IRF(t) SS(t)

Frequency domain IRF(𝜔) TF(𝜔)

Figure 2.4: Strategy to obtain state-space representation

2.3.2. IRF in the time domain

The model used in OrcaFlex [23] is formulated by:

According to Journée [16], the discretized model of the IRF is:

105 IRF Surge 106 IRF Sway

106 IRF Heave 107 IRF Roll

109 IRF Pitch 109 IRF Yaw

2.3.3. IRF in the frequency domain

K(𝜔) = B(𝜔) + 𝑖𝜔[A(𝜔) − A ] (2.10)

105 Surge added mass 106 Sway added mass

10 AQWA data AQWA data

Added mass [kg]

107 Heave added mass 108 Roll added mass