A control strategy for combined DP station keeping and active roll reduction

Abstract

Dynamic positioning (DP) systems are used for station keeping during offshore operations. The safety and operability of several offshore operations can be increased when the roll motion is actively controlled, especially in beam seas. We propose a novel control strategy for combined roll motion control and station keeping, using no additional hardware than the installed DP thrusters. The control strategy is applied to an offshore construction vessel and the performance is demonstrated by time domain simulations. The DP footprint is compared to a conventional dynamic positioning control model. The proposed control model enables active roll reduction while the station keeping performance remains unaffected. The code has been made open source and is available on https://github.com/pwellens/3dp.git.

1. Introduction

The use of vessels equipped with a Dynamic Positioning (DP) system has become the standard for offshore operations. Numerous FPSOs, drill, cable-laying, pipe-laying, heavy-lift and offshore supply vessels are equipped with a DP system to actively control the horizontal motions of the vessel. However, the safety and operability of several DP operations, like tool overboardings and lifting operations, can be increased when also the roll motion is actively controlled. This is especially the case when the vessel is operating in beam waves. Typical active roll reduction systems such as rudder roll damping and anti-roll fins are not effective during DP operations since the speed of the vessel is near zero. The use of anti-roll tanks is effective at near zero speed, but many vessels lack such a system.

Jürgens and Palm [9] showed that it is possible to actively decrease the roll motion by using Voith Schneider propellers (VSP). The fast thrust generation of a VSP makes it highly suitable for both positioning and roll reduction purposes. Koschorrek et al. [10] proposed a method to combine DP and roll reduction using VSP. However, the major part of DP vessels is equipped with conventional thrusters, such as azimuthing thrusters and tunnel thrusters. Sørensen and Strand [21] and Xu et al. [23] showed that it is possible to actively damp the unintentional low-frequency roll-pitch motion induced by the DP system of a semi-submersible.

Rudaa et al. [17] showed that conventional thrusters, such as azimuthing and tunnel thrusters, can also be used to significantly reduce the wave-frequency roll motion of a offshore supply vessel by controlling both the shaft speed and pitch angle of controllable pitch thrusters. However, the control model was developed for roll reduction purposes only. The effect of thruster induced wave-frequency roll reduction on the DP station keeping performance and the possibility of merging both control models was left uninvestigated.

As our main contribution, therefore, we present a control strategy for combined DP station keeping and thruster induced wave-frequency roll reduction is presented. The model is called 3DP, because it adds roll reduction to dynamic position in the two-dimensional horizontal plane. In our strategy, the thrusters are used to counteract the wave-frequency roll moment induced by waves. Naturally, the power consumption of the thrusters increases significantly during roll reduction mode. Therefore, the 3DP model is not envisaged as an operational mode that is executed for long periods of time, but rather as a back-up instrument to enable critical operations, like tool overboardings or subsea cable pull-ins, that are near the operability limit. By using the active roll reduction mode in these situations, the operability of the vessel can be increased at no additional operational cost but increased fuel consumption. Maintenance cost as a result of increased wear may increase, but is not quantified in this article.

The 3DP model is applied to an offshore construction vessel. The performance of the model is investigated by numerical analysis. The vessel motions are calculated with a time domain model based on frequency domain vessel data. The vessel time domain model is coupled with a dynamic thruster model to incorporate the transient response of the thrusters. Subsequently, a combined control strategy is proposed and the performance regarding station keeping, thruster power consumption and roll reduction is compared to a conventional DP control system. Application of the 3DP control model system actively decreases the roll motion. The effect on the station keeping performance is shown to be limited. The code has been made open source and is available on https://github.com/pwellens/3dp.git.

2. Design approach

To achieve both active roll reduction and DP station keeping, a new approach regarding thrust allocation and controller structure is necessary. The main idea behind the developed control strategy is presented in the design approach.

2.1. Thrust allocation

Rudaa et al. [17] propose to use thruster pairs to counteract the wave induced roll moment. When the thruster pairs are fixed and pointing in opposing direction, this yields three main advantages:

Yaw stability, when the thruster pair consists 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.

Reduced thruster-thruster interaction, since the thruster wakes are pointed away from each other on average

An example of a thruster pair configuration is indicated by dashed lines in Fig. 1, in which the bow thruster is indicated by T1, and the azimuthing thrusters by T2 until T5. The position of the centre of gravity is indicated by CoG.

Fig. 1.

Schematic thruster pair configuration for roll reduction purposes. Dashed lines indicate a thruster pair.

The thruster pairs work in counterphase to achieve maximum 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, see Fig. 1, are already pointing in the direction of the incoming wave loads, these thrusters are used for station keeping in the sway direction.

The yaw motion is controlled by the bow thruster, 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, the azimuth angle of the thrusters used for surge control are controlled by an azimuth controller. Thruster T3 is used for compensation of positive surge forces and thruster T5 is used to compensate negative surge forces. A schematic visualization of the thrusters used for station keeping in sway, surge and yaw direction is given in Fig. 2.

Fig. 2.

Thruster configuration for station keeping in the horizontal plane (left: sway, middle:surge and right: yaw motion). Dashed lines indicate which thruster (pair) controls which motion.

2.2. Control system structure

The vessel motions are controlled by a hierarchical control system. The system consists of high-level motion controllers and low-level shaft speed ( ${\mathit{\tau }}_{\mathit{s}}$ , ${\mathit{\tau }}_{\mathrm{DP}}$ and ${\mathit{\tau }}_{\mathit{\varphi }}$ , where ϕ indicates the roll motion) and azimuth angle controllers ( ${\mathit{\tau }}_{\mathrm{azi}}$ ). The control structure is visualized in Fig. 3

Fig. 3.

Schematic 3DP control system structure.

A shaft speed controller is implemented to control the shaft speed of the dynamic thruster model. The high-level motion controllers are merged by superimposing the shaft speed controller command for DP with the shaft speed controller command for active roll reduction. The principle hereof is visualized in Fig. 4.

Fig. 4.

Schematic visualization of high-level motion control merging principle.

As illustrated in Fig. 4, merging the roll reduction controller and the DP controller results in a decrease of the roll reduction capacity by increasing the minimum number of shaft revolutions per minute (RPM).

3. Mathematical model

3.1. Kinematics

The horizontal vessel motions surge, sway and yaw as calculated by the model are based on a body-fixed reference frame $[\mathit{u}\mathrm{,}\mathit{v}\mathrm{,}\mathit{r}]$ and a earth-fixed reference frame $[\mathit{x}\mathrm{,}\mathit{y}\mathrm{,}\mathit{\psi }]$ . The vessel motions, velocities and accelerations are calculated in the body-fixed frame and are subsequently translated to the earth-fixed frame by using the rotation matrix: $$\begin{array}{}\text{(1)}& \mathbf{R}\mathrm{(}\mathit{\psi }\mathrm{)}=\left[\begin{array}{ccc}\hfill cos\mathrm{(}\mathit{\psi }\mathrm{)}\hfill & \hfill -sin\mathrm{(}\mathit{\psi }\mathrm{)}\hfill & \hfill 0\hfill \\ \hfill sin\mathrm{(}\mathit{\psi }\mathrm{)}\hfill & \hfill cos\mathrm{(}\mathit{\psi }\mathrm{)}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right]\end{array}$$ The reference frames are illustrated in Fig. 5. Also the orientation of the environmental loads μ is indicated.

Fig. 5.

Definition of reference frames and environmental load direction μ.

3.2. Time domain model

In order to calculate the vessel’s (linearized) motion response when subjected to (non-linear) forces and moments induced by the thrusters and the environment, a 6 degree of freedom time domain model is used. The time domain model is based on Cummins [3]: $$\begin{array}{}\text{(2)}& [\mathbf{M}+\mathbf{A}\mathrm{(}\infty \mathrm{)}]\ddot{\mathit{\eta }}+{{\displaystyle \int }}_0^{\infty }\mathbf{K}\mathrm{(}\mathit{t}-\mathit{\tau }\mathrm{)}\dot{\mathit{\eta }}\mathit{d}\mathit{\tau }+\mathbf{C}\mathit{\eta }={\mathbf{F}}_{\mathrm{ex}}\mathrm{(}\mathit{t}\mathrm{)}+{\mathbf{F}}_{\mathit{t}}\mathrm{(}\mathit{t}\mathrm{)}\end{array}$$ where $$\begin{array}{}\text{(3)}& \mathbf{K}\mathrm{(}\mathit{\tau }\mathrm{)}={\displaystyle \frac{2}{\mathit{\pi }}}{{\displaystyle \int }}_0^{\infty }\mathbf{B}\mathrm{(}\mathit{\omega }\mathrm{)}cos\mathrm{(}\mathit{\omega }\mathit{\tau }\mathrm{)}\mathit{d}\mathit{\omega }\mathrm{,}\text{(4)}& \mathbf{A}\mathrm{(}\infty \mathrm{)}=\mathbf{A}\mathrm{(}\mathit{\omega }\mathrm{)}+{\displaystyle \frac{1}{\mathit{\omega }}}{{\displaystyle \int }}_0^{\infty }\mathbf{K}\mathrm{(}\mathit{\tau }\mathrm{)}sin\mathrm{(}\mathit{\omega }\mathit{\tau }\mathrm{)}\mathit{d}\mathit{\tau }\end{array}$$ where M, B, A, ω, $\mathbf{A}\mathrm{(}\infty \mathrm{)}$ , η, $\mathbf{K}\mathrm{(}\mathit{t}\mathrm{)}$ , C, ${\mathbf{F}}_{\mathrm{ex}}\mathrm{(}\mathit{t}\mathrm{)}$ and ${\mathbf{F}}_{\mathit{t}}\mathrm{(}\mathit{t}\mathrm{)}$ are the vessel inertia matrix, hydrodynamic damping matrix, added mass matrix, wave frequency, infinite added mass matrix, motion vector, retardation or impulse response function (IRF) matrix, restoring coefficient matrix, the environmental forces and moments and the forces and moments induced by the thrusters. The coefficient matrices in the time domain model can be derived from panel method results in the frequency domain.

There exist different methods to evaluate the convolution term in the Cummins’ equation, see Armesto et al. [1]. A state-space method is used, due to its favourable computational performance. It is possible to approximate the convolution operation in Cummins’ equation by a state-space model: $$\begin{array}{}\text{(5)}& \mathit{\kappa }={{\displaystyle \int }}_0^{\mathit{t}}\mathbf{K}\mathrm{(}\mathit{t}-\mathit{\tau }\mathrm{)}\dot{\mathit{\eta }}\mathit{d}\mathit{\tau }\simeq \left|\begin{array}{cc}\hfill & \dot{\mathit{x}}={\mathbf{A}}_{\mathrm{ss}}\mathit{x}+{\mathbf{B}}_{\mathrm{ss}}\dot{\mathit{\eta }}\hfill \\ \hfill & \mathit{\kappa }={\mathbf{C}}_{\mathrm{ss}}\mathit{x}\hfill \end{array}\right.\end{array}$$ where $\dot{\mathit{x}}$ , ${\mathbf{A}}_{\mathrm{ss}}$ , ${\mathbf{B}}_{\mathrm{ss}}$ , ${\mathbf{C}}_{\mathrm{ss}}$ represent the state vector, system matrix, input matrix and the output matrix. The strategy to obtain the state-space representation κ is visualized in Fig. 6, in which IRF represents an impulse response function, TF represents a transfer function and SS represents a state-space model.

Fig. 6.

Schematic visualization of the strategy used to obtain a time domain state-space representation.

The impulse response function is calculated numerically in the time domain for every motion and coupling term according to Journée [8]: $$\begin{array}{rcl}IRF\mathrm{(}\mathit{t}\mathrm{)}& =& {\displaystyle \frac{2}{\mathit{\pi }{\mathit{t}}^2}}{\displaystyle \sum }_{\mathit{n}=1}^{\mathit{N}}\mathrm{\left(}{\displaystyle \frac{\mathrm{\Delta }{\mathit{B}}_{\mathit{n}}}{\mathrm{\Delta }\mathit{\omega }}}\right[cos\mathrm{(}{\mathit{\omega }}_{\mathit{n}}\mathit{t}\mathrm{)}-cos\mathrm{(}{\mathit{\omega }}_{\mathit{n}-1}\mathit{t}\mathrm{)}\left]\mathrm{\right)}\\ \text{(6)}& & & +{\displaystyle \frac{2}{\mathit{\pi }\mathit{t}}}{\mathit{B}}_{\mathit{N}}sin\mathrm{(}{\mathit{\omega }}_{\mathit{N}}\mathit{t}\mathrm{)}\end{array}$$ where ${\mathit{B}}_{\mathit{n}}$ and N represent the hydrodynamic damping term at frequency n and the total number of evaluated frequencies.

Subsequently, the impulse response function is transferred to the frequency domain according to Duarte and Sarmento [5]: $$\begin{array}{}\text{(7)}& \mathbf{IRF}\mathrm{(}\mathit{\omega }\mathrm{)}=\mathbf{B}\mathrm{(}\mathit{\omega }\mathrm{)}+\mathit{i}\mathit{\omega }[\mathbf{A}\mathrm{(}\mathit{\omega }\mathrm{)}-\mathbf{A}\mathrm{(}\infty \mathrm{)}]\end{array}$$ Subsequently, a transfer function TF(ω) is fitted through the retardation function by using a least squares method and a weight factor as described in Duarte and Sarmento [5]. By using this approach, the impulse response function of each motion and coupling term is fitted by a transfer function (TF) and converted to a state-space representation SS(t). The time domain model with state-space representation of the convolution term is defined by: $$\begin{array}{}\text{(8)}& [\mathbf{M}+\mathbf{A}\mathrm{(}\infty \mathrm{)}]\ddot{\mathit{\eta }}+{\mathbf{B}}_{\mathrm{SS}}\dot{\mathit{\eta }}+\mathbf{C}\mathit{\eta }={\mathbf{F}}_{\mathrm{ex}}\mathrm{(}\mathit{t}\mathrm{)}+{\mathbf{F}}_{\mathit{t}}\mathrm{(}\mathit{t}\mathrm{)}\end{array}$$ where ${\mathbf{B}}_{\mathrm{SS}}$ represents the state-space representation matrix. Results of the time domain model with state-space representation of the convolution term are compared to frequency domain results further below in Section 5. The frequency domain model is defined by Newman [13]: $$\begin{array}{}\text{(9)}& \mathrm{(}-{\mathit{\omega }}^2[\mathbf{M}+\mathbf{A}\mathrm{(}\mathit{\omega }\mathrm{)}]-\mathit{\omega }\mathbf{B}\mathrm{(}\mathit{\omega }\mathrm{)}+\mathbf{C}\mathrm{)}\mathit{\eta }\mathrm{(}\mathit{\omega }\mathrm{)}=\mathbf{F}\mathrm{(}\mathit{\omega }\mathrm{)}\end{array}$$ where $\mathbf{F}\mathrm{(}\mathit{\omega }\mathrm{)}$ represents the wave force as a function of wave frequency.

3.3. Viscous roll damping

Viscous effects and energy dissipation are neglected when diffraction analysis based on linear potential theory is used. Since the motions of a ship other than roll are dominated by potential damping, diffraction analysis is, in general, sufficiently accurate at predicting the motion response. For the roll motion, however, viscous effects are dominant. Vortex shedding at the bilge induces damping of the ship.

According to Chakrabarti [2], the viscous roll damping term can be expressed by a third-order polynomial: $$\begin{array}{}\text{(10)}& \mathit{B}\mathrm{(}\dot{{\mathit{\eta }}_4}\mathrm{)}={\mathit{B}}_1\dot{{\mathit{\eta }}_4}+{\mathit{B}}_2|\dot{{\mathit{\eta }}_4}|\dot{{\mathit{\eta }}_4}+{\mathit{B}}_3{\dot{{\mathit{\eta }}_4}}^3\end{array}$$ where ${\mathit{B}}_1$ , ${\mathit{B}}_2$ and ${\mathit{B}}_3$ represent the linear, quadratic and cubic damping coefficients. It is not straightforward to evaluate the higher-order terms in a frequency domain model. Therefore, the viscous damping terms are included in the frequency domain by assuming an approximated equivalent linear damping coefficient: $$\begin{array}{}\text{(11)}& \mathit{B}\mathrm{(}\dot{{\mathit{\eta }}_4}\mathrm{)}={\mathit{B}}_{\mathrm{eq}}\dot{{\mathit{\eta }}_4}=\mathrm{(}{\mathit{B}}_{\mathit{f}}+{\mathit{B}}_{\mathit{e}}+{\mathit{B}}_{\mathit{w}}+{\mathit{B}}_{\mathit{L}}+{\mathit{B}}_{\mathrm{BK}}\mathrm{)}\dot{{\mathit{\eta }}_4}\end{array}$$ where ${\mathit{B}}_{\mathit{f}}$ , ${\mathit{B}}_{\mathit{e}}$ , ${\mathit{B}}_{\mathit{w}}$ , ${\mathit{B}}_{\mathit{L}}$ and ${\mathit{B}}_{\mathrm{BK}}$ represent the hull skin friction damping, hull eddy shedding damping, radiation damping, lift force damping and bilge keel damping coefficients. The radiation damping is accurately computed by diffraction analysis and since the vessel is operating at zero-speed during DP operations, the lift force damping is assumed zero. The equivalent viscous damping coefficient is therefore reduced to: $$\begin{array}{}\text{(12)}& {\mathit{B}}_{\mathrm{eq}}=\mathrm{(}{\mathit{B}}_{\mathit{f}}+{\mathit{B}}_{\mathit{e}}+{\mathit{B}}_{\mathrm{BK}}\mathrm{)}\end{array}$$ The empirical model as proposed by Ikeda et al. [7] is used to calculate the eddy shedding damping. The hull skin friction and bilge keel damping coefficients are obtained by using the empirical models in Himeno [6]. The empirical models for zero-speed condition are used. The damping coefficient models included in ${\mathit{B}}_{\mathrm{eq}}$ are a function of both the vessel’s roll amplitude and wave frequency. The magnitude of the damping coefficient has a direct effect on the resulting roll amplitude, since an additional damping term is added. Therefore, first the equivalent viscous damping is calculated according to the roll amplitudes without inclusion of the viscous damping term. The calculated viscous damping term per wave frequency is subsequently included in (9), in which ${\mathit{B}}_{\mathrm{eq}}$ for roll is added to B. The viscous damping term is recalculated according to the newly calculated roll motion amplitudes. This process is repeated until both the roll motion amplitudes and the viscous damping term have converged.

In the time domain model it is possible to incorporate a quadratic damping term ${\mathit{B}}_{\mathrm{visc}}|\dot{{\mathit{\eta }}_4}|\dot{{\mathit{\eta }}_4}$ , with ${\mathit{B}}_{\mathrm{visc}}$ a tuning parameter that is found by iterating the time domain roll motion Response Amplitude Operator (RAO) until it matches the frequency domain roll motion RAO for a unit wave amplitude. The time domain model with the viscous damping term included in the $\mathrm{(}4\mathrm{,}4\mathrm{)}$ position of ${\mathbf{B}}_{\mathrm{visc}}$ is now expressed by: $$\begin{array}{}\text{(13)}& [\mathbf{M}+\mathbf{A}\mathrm{(}\mathbf{\infty }\mathrm{)}]\ddot{\mathit{\eta }}+{\mathbf{B}}_{\mathrm{ss}}\dot{\mathit{\eta }}+{\mathbf{B}}_{\mathrm{visc}}|\dot{\mathit{\eta }}|\dot{\mathit{\eta }}+\mathbf{C}\mathit{\eta }={\mathbf{F}}_{\mathrm{ex}}\mathrm{(}\mathit{t}\mathrm{)}+{\mathbf{F}}_{\mathit{t}}\mathrm{(}\mathit{t}\mathrm{)}\end{array}$$

3.4. Environmental forces

The vessel motions are induced by waves, wind and current. The wave force consists of a first- and a second-order component. The first-order wave force (indicated by superscript $\mathrm{(}1\mathrm{)}$ ) is expressed in the time domain as: $$\begin{array}{}\text{(14)}& {\mathit{F}}^{\mathrm{(}1\mathrm{)}}\mathrm{(}\mathit{t}\mathrm{)}={\displaystyle \sum }_{\mathit{n}=1}^{\mathit{N}}\sqrt{2·{\mathit{S}}_{\mathit{\eta }\mathrm{,}\mathit{F}}\mathrm{(}\mathit{\omega }\mathrm{)}\mathrm{\Delta }\mathit{\omega }}cos\mathrm{(}{\mathit{\omega }}_{\mathit{n}}\mathit{t}+{\mathit{ϵ}}_{\mathit{n}}+{\mathit{ϵ}}_{\mathit{\eta }\mathrm{,}\mathit{F}}\mathrm{)}\end{array}$$ where N, ${\mathit{ϵ}}_{\mathit{n}}$ and ${\mathit{ϵ}}_{\mathit{\eta }\mathrm{,}\mathit{F}}$ represent the total number of frequencies evaluated by diffraction analysis, random phase angle chosen in the interval $[02\mathit{\pi }]$ and phase shift of the wave load. The wave force RAO ${\mathit{S}}_{\mathit{\eta }\mathrm{,}\mathit{F}}\mathrm{(}\mathit{\omega }\mathrm{)}$ is defined by: $$\begin{array}{}\text{(15)}& {\mathit{S}}_{\mathit{\eta }\mathrm{,}\mathit{F}}\mathrm{(}\mathit{\omega }\mathrm{)}={|{\displaystyle \frac{{\mathit{F}}_{\mathit{\eta }}}{\mathit{\xi }}}\mathrm{(}\mathit{\omega }\mathrm{)}|}^2{\mathit{S}}_{\mathit{\xi }}\mathrm{(}\mathit{\omega }\mathrm{)}\end{array}$$ where ξ, $\frac{{\mathit{F}}_{\mathit{\eta }}}{\mathit{\xi }}$ and ${\mathit{S}}_{\mathit{\xi }}$ represent the free surface amplitude, the wave force RAO and the wave amplitude spectrum. The wave amplitude spectrum is characterized by a significant wave height ${\mathit{H}}_{\mathit{s}}$ and a peak period ${\mathit{T}}_{\mathit{p}}$ .

The low-frequency and mean wave drift force (indicated by superscript $\mathrm{(}2\mathrm{)}$ ) in irregular waves are calculated in the time domain according to Newman [13] by: $$\begin{array}{rcl}{\mathit{F}}^{\mathrm{(}2\mathrm{)}}\mathrm{(}\mathit{t}\mathrm{)}& =& {\displaystyle \sum }_{\mathit{i}}{\displaystyle \sum }_{\mathit{j}}{\mathit{\xi }}_{\mathit{i}}{\mathit{\xi }}_{\mathit{j}}{\mathit{P}}_{\mathit{i}\mathit{j}}^{-}cos\mathrm{(}\mathrm{(}{\mathit{\omega }}_{\mathit{i}}-{\mathit{\omega }}_{\mathit{j}}\mathrm{)}\mathit{t}+\mathrm{(}{\mathit{ϵ}}_{\mathit{i}}-{\mathit{ϵ}}_{\mathit{j}}\mathrm{)}\mathrm{)}\\ \text{(16)}& & & +{\displaystyle \sum }_{\mathit{i}}{\displaystyle \sum }_{\mathit{j}}{\mathit{\xi }}_{\mathit{i}}{\mathit{\xi }}_{\mathit{j}}{\mathit{Q}}_{\mathit{i}\mathit{j}}^{-}sin\mathrm{(}\mathrm{(}{\mathit{ϵ}}_{\mathit{i}}-{\mathit{ϵ}}_{\mathit{j}}\mathrm{)}\mathit{t}+\mathrm{(}{\mathit{ϵ}}_{\mathit{i}}-{\mathit{ϵ}}_{\mathit{j}}\mathrm{)}\mathrm{)}\end{array}$$ where ${\mathit{P}}^{-}$ , ${\mathit{Q}}^{-}$ and ϵ represent the in-phase quadratic transfer function (QTF), out-of-phase QTF and phase angle. The QTFs ${\mathit{Q}}^{-}$ , ${\mathit{P}}^{-}$ and corresponding phase angles are obtained by diffraction analysis using the near-field approach, see Pinkster [16]. The wind forces and moment for surge, sway and yaw are calculated according to Serraris [19]: $$\begin{array}{}\text{(17)}& {\mathit{F}}_{\mathit{x}\mathrm{,}\mathrm{wind}}={\mathit{C}}_{\mathit{x}\mathit{w}}\mathrm{(}\mathit{\mu }\mathrm{)}{\displaystyle \frac{1}{2}}{\mathit{\rho }}_{\mathit{a}}{\mathit{v}}_{\mathit{w}}^2{\mathit{A}}_{\mathit{T}}\text{(18)}& {\mathit{F}}_{\mathit{y}\mathrm{,}\mathrm{wind}}={\mathit{C}}_{\mathit{y}\mathit{w}}\mathrm{(}\mathit{\mu }\mathrm{)}{\displaystyle \frac{1}{2}}{\mathit{\rho }}_{\mathit{a}}{\mathit{v}}_{\mathit{w}}^2{\mathit{A}}_{\mathit{L}}\text{(19)}& {\mathit{M}}_{\mathit{z}\mathrm{,}\mathrm{wind}}={\mathit{C}}_{\mathit{z}\mathit{w}}\mathrm{(}\mathit{\mu }\mathrm{)}{\displaystyle \frac{1}{2}}{\mathit{\rho }}_{\mathit{a}}{\mathit{v}}_{\mathit{w}}^2{\mathit{A}}_{\mathit{L}}\mathit{L}\end{array}$$ where ${\mathit{C}}_{\mathit{w}}$ , μ, ${\mathit{\rho }}_{\mathit{a}}$ , ${\mathit{v}}_{\mathit{w}}$ , ${\mathit{A}}_{\mathit{T}}$ , ${\mathit{A}}_{\mathit{L}}$ and L represent the wind force coefficient, environmental direction, air density, wind speed, transverse wind area, lateral wind area and vessel length.

A similar approach is used to calculate the forces and moment due to the current: $$\begin{array}{}\text{(20)}& {\mathit{F}}_{\mathit{x}\mathrm{,}\mathrm{current}}={\mathit{C}}_{\mathit{x}\mathit{c}}\mathrm{(}\mathit{\mu }\mathrm{)}{\displaystyle \frac{1}{2}}{\mathit{\rho }}_{\mathit{w}}{\mathit{v}}_{\mathit{c}}^2\mathit{L}\mathit{T}\text{(21)}& {\mathit{F}}_{\mathit{y}\mathrm{,}\mathrm{current}}={\mathit{C}}_{\mathit{y}\mathit{c}}\mathrm{(}\mathit{\mu }\mathrm{)}{\displaystyle \frac{1}{2}}{\mathit{\rho }}_{\mathit{w}}{\mathit{v}}_{\mathit{c}}^2\mathit{L}\mathit{T}\text{(22)}& {\mathit{M}}_{\mathit{z}\mathrm{,}\mathrm{current}}={\mathit{C}}_{\mathit{z}\mathit{c}}\mathrm{(}\mathit{\mu }\mathrm{)}{\displaystyle \frac{1}{2}}{\mathit{\rho }}_{\mathit{w}}{\mathit{v}}_{\mathit{c}}^2{\mathit{L}}^2\mathit{T}\end{array}$$ where ${\mathit{C}}_{\mathit{c}}$ , ${\mathit{\rho }}_{\mathit{w}}$ , ${\mathit{v}}_{\mathit{c}}$ and T represent the current force coefficient, water density, current velocity and vessel draft. The wind and current force coefficients can be obtained by conducting wind tunnel experiments or by estimation using empirical models, see for example Nienhuis [14]. The total environmental forces and moments acting on the vessel model are expressed by: $$\begin{array}{}\text{(23)}& {\mathbf{F}}_{\mathrm{ex}}\mathrm{(}\mathit{t}\mathrm{)}={\mathit{F}}^{\mathrm{(}1\mathrm{)}}\mathrm{(}\mathit{t}\mathrm{)}+{\mathit{F}}^{\mathrm{(}2\mathrm{)}}\mathrm{(}\mathit{t}\mathrm{)}+{\mathit{F}}_{\mathrm{wind}}+{\mathit{F}}_{\mathrm{current}}\end{array}$$

3.5. Dynamic 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. Since fixed-pitch propellers are assumed, the thruster dynamics are governed by the inertia of the thruster system, propeller torque demand, torque produced by the electrical motor and shaft friction. This can be formulated according the simulation model proposed by Smogeli [20]: $$\begin{array}{}\text{(24)}& {\mathit{I}}_{\mathit{s}}\dot{\mathit{n}}={\mathit{Q}}_{\mathit{m}}-{\mathit{Q}}_{\mathit{p}}\mathrm{(}\mathit{\beta }\mathrm{,}\mathit{n}\mathrm{)}-{\mathit{Q}}_{\mathit{f}}\mathrm{(}\mathit{n}\mathrm{)}\text{(25)}& {\dot{\mathit{Q}}}_{\mathit{m}}={\displaystyle \frac{1}{{\mathit{T}}_{\mathit{m}}}}\mathrm{(}{\mathit{Q}}_{\mathit{c}\mathit{m}}-{\mathit{Q}}_{\mathit{m}}\mathrm{)}\end{array}$$ where ${\mathit{I}}_{\mathit{s}}$ , n, ${\mathit{Q}}_{\mathit{m}}$ , ${\mathit{Q}}_{\mathit{p}}$ , ${\mathit{Q}}_{\mathit{f}}$ , ${\mathit{T}}_{\mathit{m}}$ and ${\mathit{Q}}_{\mathit{c}\mathit{m}}$ represent the rotational inertia of the thruster system, propeller rotational speed, torque delivered by the motor, torque demand by the propeller, shaft friction, motor time constant and commanded motor torque. Note that he commanded motor torque needs to stay within the confines of the maximum torque a specific motor can deliver. In our case, we limited the torque at the rated torque of a typical motor.

The hydrodynamic pitch angle β is defined by: $$\begin{array}{}\text{(26)}& \mathit{\beta }={\displaystyle \frac{{\mathit{v}}_{\mathit{a}}}{0.7\mathit{\pi }\mathit{n}\mathit{D}}}\end{array}$$ where ${\mathit{v}}_{\mathit{a}}$ is the advance velocity and D is the propeller diameter. The advance velocity ${\mathit{v}}_{\mathit{a}}$ is based on current and thruster velocity.

The shaft friction term is defined by [20]: $$\begin{array}{}\text{(27)}& {\mathit{Q}}_{\mathit{f}}\mathrm{(}\mathit{n}\mathrm{)}=sign\mathrm{(}\mathit{n}\mathrm{)}{\mathit{Q}}_{\mathit{s}}+{\mathit{K}}_{\mathit{n}}\mathit{n}\end{array}$$ where ${\mathit{Q}}_{\mathit{s}}$ and ${\mathit{K}}_{\mathit{n}}$ represent the static shaft friction and the linear shaft friction coefficient.

Smogeli [20] acknowledges that an added inertia term due to the hydrodynamic forces in phase with the propeller rotational acceleration exists, but chooses to neglect it. An involved method based on lifting lines for the hydrodynamic properties of open water propellors is available in Krüger and Abels [11]. An empirical estimate of the inertia term of open water propellors can be obtained from MacPherson et al. [12] and Schwanecke [18]. In reality, DP propellors are ducted and below the ship. Methods or data for those circumstances are – to our knowledge – not readily available. Not knowing how close either method is to our specific situation, we chose to base our inertia term on the lesser involved empirical methods of MacPherson et al. [12] and Schwanecke [18] as follows.

In the latter two references, the added rotational inertia term of the entrained water is defined by: $$\begin{array}{}\text{(28)}& {\mathit{I}}_{\mathit{E}}={\mathit{C}}_{\mathrm{IE}}\mathit{\rho }{\mathit{D}}^5\end{array}$$ where ρ is the water density and ${\mathit{C}}_{\mathrm{IE}}$ is a fitting parameter defined by Schwanecke [18]: $$\begin{array}{}\text{(29)}& {\mathit{C}}_{\mathrm{IE}}={\displaystyle \frac{0.0703{\mathrm{(}\mathit{P}\mathrm{/}\mathit{D}\mathrm{)}}^2{\mathrm{EAR}}^2}{\mathit{\pi }\mathit{Z}}}\end{array}$$ where $\mathit{P}\mathrm{/}\mathit{D}$ , $\mathrm{EAR}$ and Z represent the propeller pitch-diameter ratio, expanded blade area ratio and the total number of propeller blades. MacPherson et al. [12] propose the following expression for the fitting parameter: $$\begin{array}{}\text{(30)}& {\mathit{C}}_{\mathrm{IE}}={\mathit{C}}_1\mathrm{EAR}\mathrm{(}\mathit{P}\mathrm{/}\mathit{D}\mathrm{)}-{\mathit{C}}_2\end{array}$$ where ${\mathit{C}}_1$ and ${\mathit{C}}_2$ are fitting parameters given in Table 1.

Table 1

${\mathit{C}}_{\mathit{IE}}$ fitting coefficients according MacPherson et al. [12]

Parameter	$\mathit{Z}=3$	$\mathit{Z}=4$	$\mathit{Z}=5$	$\mathit{Z}=6$
${\mathit{C}}_1$	0.00477	0.00394	0.00359	0.00344
${\mathit{C}}_2$	0.00093	0.00087	0.00080	0.00076

Since both empirical models are based on model test results carried out with different propeller types, the average of the model results is used as a representative value for the inertia term.

The thrusters experience oscillating inflow velocities as a result of the roll motion of the vessel. Therefore, the four-quadrant model as developed by Oosterveld [15] is used to calculate the thruster torque demand and thrust production: $$\begin{array}{}\text{(31)}& {\mathit{Q}}_{\mathit{p}}\mathrm{(}\mathit{\beta }\mathrm{,}\mathit{n}\mathrm{)}={\mathit{C}}_{\mathit{Q}}\mathrm{(}\mathit{\beta }\mathrm{)}\mathit{\rho }\mathrm{(}{\mathit{v}}_{\mathit{a}}^2+{\mathrm{(}0.7\mathit{\pi }\mathit{n}\mathit{D}\mathrm{)}}^2\mathrm{)}{\displaystyle \frac{\mathit{\pi }}{8}}{\mathit{D}}^3\text{(32)}& {\mathit{T}}_{\mathit{p}}\mathrm{(}\mathit{\beta }\mathrm{,}\mathit{n}\mathrm{)}={\mathit{C}}_{\mathit{T}}\mathrm{(}\mathit{\beta }\mathrm{)}\mathit{\rho }\mathrm{(}{\mathit{v}}_{\mathit{a}}^2+{\mathrm{(}0.7\mathit{\pi }\mathit{n}\mathit{D}\mathrm{)}}^2\mathrm{)}{\displaystyle \frac{\mathit{\pi }}{8}}{\mathit{D}}^2\text{(33)}& {\mathit{C}}_{\mathit{Q}}\mathrm{(}\mathit{\beta }\mathrm{)}\approx {\displaystyle \sum }_{\mathit{k}=0}^{20}[{\mathit{A}}_{\mathit{Q}}\mathrm{(}\mathit{k}\mathrm{)}cos\mathrm{(}\mathit{\beta }\mathit{k}\mathrm{)}+{\mathit{B}}_{\mathit{Q}}\mathrm{(}\mathit{k}\mathrm{)}sin\mathrm{(}\mathit{\beta }\mathit{k}\mathrm{)}]\text{(34)}& {\mathit{C}}_{\mathit{T}}\mathrm{(}\mathit{\beta }\mathrm{)}\approx {\displaystyle \sum }_{\mathit{k}=0}^{20}[{\mathit{A}}_{\mathit{T}}\mathrm{(}\mathit{k}\mathrm{)}cos\mathrm{(}\mathit{\beta }\mathit{k}\mathrm{)}+{\mathit{B}}_{\mathit{T}}\mathrm{(}\mathit{k}\mathrm{)}sin\mathrm{(}\mathit{\beta }\mathit{k}\mathrm{)}]\end{array}$$ where the Fourier coefficients ${\mathit{A}}_{\mathit{Q}}$ , ${\mathit{B}}_{\mathit{Q}}$ , ${\mathit{A}}_{\mathit{T}}$ and ${\mathit{B}}_{\mathit{T}}$ are based on experimental data resulting from ducted propeller tests obtained by Oosterveld [15]. A thrust deduction factor $\mathit{t}=0.04$ is included to take thruster-hull interaction into account, see Wichers et al. [22].

The roll moment of the thrusters is calculated by multiplying the produced thrust in the sway direction with the thruster moment arm. The thruster moment arm is defined as the vertical distance from the centre of the thruster to the centre of gravity (CoG) of the vessel.

3.6. 3DP control strategy

The high level motion controller is a combination of a conventional DP controller and roll controller: $$\begin{array}{}\text{(35)}& \mathit{\tau }={\mathbf{H}}_1{\mathit{\tau }}_{\mathrm{DP}}+{\mathbf{H}}_2{\mathit{\tau }}_{\mathit{\varphi }}\end{array}$$ where $$\begin{array}{}\text{(36)}& {\mathit{\tau }}_{\mathrm{DP}}={\mathit{K}}_{\mathit{p}}{\mathit{e}}_{\mathrm{DP}}+{\mathit{K}}_{\mathit{i}}{{\displaystyle \int }}_0^{\mathit{t}}{\mathit{e}}_{\mathrm{DP}}\mathrm{(}\mathit{t}\mathrm{)}\mathit{d}\mathrm{(}\mathit{t}\mathrm{)}+{\mathit{K}}_{\mathit{d}}{\dot{\mathit{e}}}_{\mathrm{DP}}\text{(37)}& {\mathit{\tau }}_{\mathit{\varphi }}={\mathit{K}}_{\mathit{p}}\dot{\mathit{\varphi }}\end{array}$$ where ${\mathit{K}}_{\mathit{p}}$ , ${\mathit{K}}_{\mathit{i}}$ , ${\mathit{K}}_{\mathit{d}}$ and $\dot{\mathit{\varphi }}$ represent the proportional gain, integral gain, the derivative gain and the roll velocity. The error values ${\mathit{e}}_{\mathrm{DP}}$ are defined by: $$\begin{array}{}\text{(38)}& {\mathit{e}}_{\mathrm{DP}}=\left[\begin{array}{c}\hfill {\mathit{x}}_{\mathit{s}}-\mathit{x}\hfill \\ \hfill {\mathit{y}}_{\mathit{s}}-\mathit{y}\hfill \\ \hfill {\mathit{\psi }}_{\mathit{s}}-\mathit{\psi }\hfill \end{array}\right]\end{array}$$ where ${\mathit{x}}_{\mathit{s}}$ , ${\mathit{y}}_{\mathit{s}}$ and ${\mathit{\psi }}_{\mathit{s}}$ are the surge, sway and yaw setpoint. ${\mathbf{H}}_1$ and ${\mathbf{H}}_2$ represent allocation matrices. These matrices allocate the required thrust to the available thrusters by a corresponding RPM command. The thrusters are limited to a maximum and a minimum RPM to avoid overloading and to increase thruster reaction time, therefore: $$\begin{array}{}\text{(39)}& {\mathit{n}}_{min}⩽\mathit{\tau }⩽{\mathit{n}}_{max}\end{array}$$ The commanded RPM τ is used subsequently as setpoint input for the lower-level shaft speed controller.

The thruster shaft speed is controlled by a PI controller as proposed by Smogeli [20]: $$\begin{array}{}\text{(40)}& {\mathit{\tau }}_{\mathit{s}}={\mathit{Q}}_{\mathit{c}\mathit{m}}={\mathit{K}}_{\mathit{p}}{\mathit{e}}_{\mathit{s}}+{\mathit{K}}_{\mathit{i}}{{\displaystyle \int }}_0^{\mathit{t}}{\mathit{e}}_{\mathit{s}}\mathrm{(}\mathit{t}\mathrm{)}\mathit{d}\mathrm{(}\mathit{t}\mathrm{)}\end{array}$$ The shaft speed error term ${\mathit{e}}_{\mathit{s}}$ is calculated by: $$\begin{array}{}\text{(41)}& {\mathit{e}}_{\mathit{s}}=\mathrm{(}\mathit{\tau }-{\mathit{n}}_{\mathit{s}}\mathrm{)}\end{array}$$ where ${\mathit{n}}_{\mathit{s}}$ represents the actual shaft speed.

In beam waves, the wave forces acting on the vessel in the surge direction are negligible. However, wind and current can induce significant forces in the surge direction when their direction differs from the wave direction. To increase the station keeping capability of the vessel, the azimuth angle of the thrusters allocated to control the surge motion is controlled by a proportional controller: $$\begin{array}{}\text{(42)}& {\mathit{\tau }}_{\mathrm{azi}}={\mathit{K}}_{\mathit{p}}{\mathit{T}}_{\mathit{x}}\end{array}$$ where ${\mathit{T}}_{\mathit{x}}$ is the thrust command in the surge direction from the DP controller. By using a proportional azimuth angle controller, the controller behaves like a mechanical spring, adjusting the azimuth angle when the environmental load in the surge direction increases. Another advantage of using a proportional azimuth angle controller is the fact that the thrusters are aligned to achieve maximum roll reduction when no environmental loads in the surge direction are present. The maximum azimuth angle is set at ${45}^{\circ }$ to ensure stability of the system.

3.7. Controller tuning

Tuning of the controllers is important to ensure stable vessel behaviour. A first estimate of the proportional, integral and derivative gains of the high-level DP controller ${\mathit{\tau }}_{\mathrm{DP}}$ is obtained by the rules of thumb proposed by Serraris [19].

The first estimates of the proportional and integral gain terms of the low-level shaft speed controller ${\mathit{\tau }}_{\mathit{s}}$ are based on parameter estimates as presented in the work of Smogeli [20]. The controller gains of both controllers are tuned manually until the desired stable system behaviour was obtained.

The high-level anti-roll controller and the low-level azimuth angle controller are tuned by conducting systematic numerical experiments. In this way both controllers are tuned such that roll reduction is at maximum.

3.8. DP model

The station keeping performance of the control model is compared to a conventional DP control system to analyze the effect of active roll reduction onto the DP footprint. The conventional DP system consists of a Kalman filter, PID controller and an allocation algorithm [19]. The coefficients of the Kalman filter and the PID controller of the conventional DP system are kept the same as in our proposed roll reduction control model to ensure equal comparison.

4. Numerical experiment description

4.1. Vessel particulars

The performance of the control model is analyzed by applying the 3DP model to a barge-shaped offshore construction vessel. The vessel particulars are given in Table 2.

Table 2

Vessel particulars offshore construction vessel

Parameter	Value
L [m]	99
B [m]	30
T [m]	4.7
∇ [ ${\text{m}}^3$ ]	11 683
${\mathit{T}}_{\mathit{\varphi }}$ [s]	8.5

In Table 2, L is the length of the vessel, B the width, T is the draft, ∇ the displacement and ${\mathit{T}}_{\mathit{\varphi }}$ is the natural roll period. Barge-shaped vessels are particularly sensitive to roll motions due to their large $\mathit{B}\mathrm{/}\mathit{T}$ -ratio.

The vessel is equipped with 5 thrusters in total, of which 4 azimuthing thrusters and 1 bow thruster. The thruster lay-out is given in Fig. 1. The DP control point is located in the vessel’s CoG. The DP setpoint is defined as $[0\mathrm{,}0\mathrm{,}0]$ .

The thruster particulars are given in Table 3.

Table 3

Thruster particulars of offshore construction vessel

Thruster	Power [kW]	Thrust [kN]	Type	D [m]
T1	500	75	tunnel	1.0
T2	1000	177	azimuth	1.75
T3	1000	177	azimuth	1.75
T4	1250	222	azimuth	2.0
T5	1250	222	azimuth	2.0

4.2. Environmental conditions

The vessel is subjected to several environmental conditions to assess the performance of the control model for roll reduction together with station keeping. The simulated environmental conditions are given in Table 4.

Table 4

Simulated environmental conditions

Case	${\mathit{H}}_{\mathit{s}}$ [m]	${\mathit{T}}_{\mathit{p}}$ [s]	${\mathit{v}}_{\mathit{w}}$ [m/s]	${\mathit{v}}_{\mathit{c}}$ [m/s]	${\mathit{M}}_{\mathit{t}}\mathrm{/}{\mathit{M}}_{\mathit{env}}$ [–]
1	0.5	8.5	0.39	0.75	1.24
2	0.75	8.5	1.94		0.83
3	1.0	8.5	3.32		0.63

It is assumed that the wind and current forces act co-linearly with the waves ( $\mathit{\mu }={90}^{\circ }$ ). The wave peak period is chosen equal to the natural roll frequency of the vessel. A JONSWAP spectrum with peak enhancement factor $\mathit{\gamma }=3.3$ is used to calculate the wave loads. The parameter ${\mathit{M}}_{\mathit{t}}\mathrm{/}{\mathit{M}}_{\mathrm{env}}$ represents a dimensionless ratio, in which ${\mathit{M}}_{\mathit{t}}$ and ${\mathit{M}}_{\mathrm{env}}$ represent the maximum roll moment induced by the thrusters and the mean of the absolute first-order wave roll moment. The total simulation time is set to 0.5 hour, since critical situation during typical operations take less time than that.

4.3. Allocation

The thrusters that are used for station keeping and roll reduction are indicated by the allocation matrices ${\mathbf{H}}_1$ and ${\mathbf{H}}_2$ . The allocation vector ${\mathbf{H}}_1$ is used to indicate which thrusters are used to control the surge, sway and yaw motions: $$\begin{array}{}\text{(43)}& {\mathbf{H}}_1=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill \frac{{\mathit{L}}_2}{{\mathit{L}}_2+{\mathit{L}}_4}\hfill & \hfill 0\hfill \\ \hfill {\mathit{X}}_3\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{{\mathit{L}}_4}{{\mathit{L}}_2+{\mathit{L}}_4}\hfill & \hfill 0\hfill \\ \hfill {\mathit{X}}_5\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\end{array}$$ where ${\mathit{L}}_2$ and ${\mathit{L}}_4$ represent the thruster moment arms of respectively thruster T2 and thruster T4, see Fig. 2. ${\mathit{X}}_3$ and ${\mathit{X}}_5$ represent two Boolean variables defined by: $$\begin{array}{ll}\hfill & {\mathit{X}}_3=\left\{\begin{array}{cc}1\hfill & \text{if }{\mathit{T}}_{\mathit{x}}⩽0\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\\ \hfill & {\mathit{X}}_5=\left\{\begin{array}{cc}1\hfill & \text{if }{\mathit{T}}_{\mathit{x}}\mathrm{>}0\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\\ \text{(44)}& \hfill & {\mathbf{H}}_2=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 1\hfill \\ \hfill 1\hfill \\ \hfill 1\hfill \\ \hfill 1\hfill \end{array}\right]\end{array}$$ where ${\mathit{T}}_{\mathit{x}}$ represent the required thrust in the surge direction to maintain position as calculated by the surge DP controller. Allocation matrix ${\mathbf{H}}_1$ indicates that thruster T1, see Fig. 1, is used for yaw control. The required thrust in the sway direction is divided over thruster T2 and T4. The division is weighed according to the thruster yaw moment arms. By doing so, the yaw moment of both thrusters 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. Allocation matrix ${\mathbf{H}}_2$ indicates that all azimuthing thrusters are used for roll reduction purposes.

5. Results and analysis

5.1. Time domain model validation

The motion response amplitude operators (RAOs) for quartering regular waves as calculated by the time domain model, without a viscous roll damping term included, are compared with frequency domain RAOs to confirm the validity of the state-space modelling approach, see Fig. 7.

Fig. 7.

Comparison of calculated RAOs by both frequency and time domain model ( $\mathit{\mu }={45}^{\circ }$ ).

From Fig. 7, we find that both models are in agreement, because the state-space representation of the convolution term yields the same results as the frequency domain model results.

5.2. Viscous roll damping

The time domain model parameter ${\mathit{B}}_{\mathrm{visc}}$ is obtained by tuning the time domain model roll motion RAO in beam waves with the frequency domain model RAO. The resulting RAOs are given in Fig. 8.

Fig. 8.

Comparison of roll RAO calculated by both frequency and time domain model, including viscous roll damping term ( $\mathit{\mu }={90}^{\circ }$ ).

It can be confirmed from the RAO visualized in Fig. 8 that the time domain tuning method results in a similar roll motion RAO for both models. By iterating, the quadratic tuning parameter in the time domain model was determined to be 6% of the critical roll damping (based on moment of inertia and added moment of inertia).

One should be aware of the fact that this roll motion RAO only matches for a unit wave amplitude. When a lower wave amplitude is used, the time domain model will underestimate the viscous damping in comparison with the frequency domain model due to the inclusion of a quadratic viscous damping term. The opposite applies when a higher wave amplitude is applied. Since the maximum simulated significant wave height is 1 m, this approach will result in less viscous damping compared to frequency domain analysis. This is considered to be a conservative approach, since the thrusters will have to counteract a bigger roll moment and will therefore reach their saturation limit more early compared to the frequency domain model.

5.3. Tuning

Systematic numerical simulations for the different environmental conditions were carried out to tune the anti-roll and azimuth controllers. The relative RMS roll reduction is calculated for different controller gains. The result of the anti-roll controller tuning procedure is given in Fig. 9. Figure 9 indicates that the roll RMS reduction percentage converges towards the higher values for ${\mathit{K}}_{\mathit{p}}$ of the roll reduction controller.

Fig. 9.

Tuning of the controller gain for roll reduction: roll RMS reduction percentage for different values ${\mathit{K}}_{\mathit{p}}$ of the anti-roll controller ( ${\mathit{T}}_{\mathit{p}}=8.5\text{ s}$ ).

To tune the azimuth controller, the relative RMS roll reduction is calculated for different directions of the environmental load. The tuning coefficients are increased systemically per direction. The results are given in Fig. 10

Fig. 10.

Tuning of the azimuth controller: roll RMS reduction percentage for different ${\mathit{K}}_{\mathit{p}}$ values of the azimuth controller and the direction of the environmental load ( ${\mathit{H}}_{\mathit{s}}=0.5\text{ m}$ , ${\mathit{T}}_{\mathit{p}}=8.5\text{ s}$ ).

From Fig. 10, it can be observed that the maximum roll reduction for all investigated directions of the environmental load occurs at a value of ${\mathit{K}}_{\mathit{p}}=0.9$ for the azimuth controller gain.

The tuning coefficients of the DP controller in the 3DP control strategy and the Kalman filter tuning parameters Q and R are given in Table 5.

Table 5

Tuning coefficients of the DP controller and Kalman filter in the 3DP model

DP Controller	Surge	Sway	Yaw
${\mathit{K}}_{\mathit{p}}$	106.1	117.9	4932.2
${\mathit{K}}_{\mathit{i}}$	0.25	0.30	0.01
${\mathit{K}}_{\mathit{d}}$	1154.7	1388.2	34 781.5
Q	$1.00\text{E}-09$	$2.00\text{E}-09$	$1.00\text{E}-09$
R	$1.30\text{E}-04$	$6.50\text{E}-05$	$1.30\text{E}-04$

The controller coefficients as used in the 3DP control model for the roll controller, azimuth controller and shaft speed controllers are summarized in Table 6. Note that a different ship would require the described procedure to be executed again to obtain new numbers for adequate performance.

Table 6

Tuning coefficients of the roll controller, azimuth controller and shaft speed controllers as implemented in the 3DP control model

Roll controller
${\mathit{K}}_{\mathit{p}}$	2.50E+06
Azimuth controller
${\mathit{K}}_{\mathit{p}}$	0.9
Shaft speed controller
${\mathit{K}}_{\mathit{p}}$	33.3
${\mathit{K}}_{\mathit{i}}$	16.65
Shaft speed controller
${\mathit{K}}_{\mathit{p}}$	19.9
${\mathit{K}}_{\mathit{i}}$	9.95

5.4. Roll reduction

To analyze the performance of the proposed roll reduction control model, the roll angle amplitudes are compared to the roll angle amplitudes obtained when roll is not actively reduced in the conventional DP model. Both models are subjected to equal environmental forces. A snapshot of the vessel roll angles as a function of time for both models is visualized in Fig. 11.

Fig. 11.

Snapshot of the roll angles time trace for simulation case 1.

As can observed in Fig. 11, the roll reduction varies in magnitude every period. As a way to visualize the realized roll reduction, the roll angle amplitudes are translated to a normal distribution. The results are given in Fig. 12.

Fig. 12.

Normal distribution of roll for simulation case 1, 2 and 3 (from top to bottom).

Shown in Fig. 12 is that the effect of active roll reduction is most significant for simulation case 1. The variance of the roll amplitude decreases and the probability density of the mean amplitude of $0^{\circ }$ increases when the 3DP model is applied. The effect of active roll reduction decreases when the significant wave height increases. To also quantify the results, the root mean square (RMS), maximum value and absolute mean of the roll angle time traces are calculated, see Table 7.

Table 7 confirms the observed behavior in Fig. 12. The roll reduction percentages decrease when the significant wave height increases. This is because the roll moment realized by the thrusters is not able to compensate for the first-order wave roll moment as the significant wave height increases. The ratio between the thruster roll moment and first-order wave roll moment ${\mathit{M}}_{\mathit{t}}\mathrm{/}{\mathit{M}}_{\mathrm{env}}$ can be seen in Table 4. The roll reduction also decreases in higher sea states because the thrusters require more power for DP station keeping purposes when the significant wave height increases.

5.5. Station keeping

The vessel’s station keeping performance can be assessed by analyzing the DP footprint of the vessel. The DP footprint visualizes the maximum excursions of the vessel in the horizontal earth-fixed reference frame with respect to the DP setpoint. To analyze the performance of the proposed control model, the DP footprint of both the proposed model and the conventional DP model is visualized in Fig. 13. Figure 13 shows that the increase of the 3DP footprint is small compared to the conventional DP footprint in Case 1 and Case 2. For Case 3, with the highest significant wave height, the footprint in y-direction is increased by 1 m, which is 1% of the vessel length.

Table 7

Roll reduction results for selected simulation cases

Case		DP Model	3DP model	Reduction
1	RMS [ ${}^{\circ }$ ]	1.11	0.69	38%
	Max. [ ${}^{\circ }$ ]	3.55	2.51	29%
	Abs. mean [ ${}^{\circ }$ ]	0.88	0.54	39%
2	RMS [ ${}^{\circ }$ ]	1.56	1.31	16%
	Max. [ ${}^{\circ }$ ]	4.86	4.33	11%
	Abs. mean [ ${}^{\circ }$ ]	1.28	1.04	19%
3	RMS [ ${}^{\circ }$ ]	2.09	1.84	12%
	Max. [ ${}^{\circ }$ ]	5.59	4.98	11%
	Abs mean [ ${}^{\circ }$ ]	1.68	1.42	15%

Fig. 13.

DP footprint of proposed control model and DP model for simulation case 1, 2 and 3 (from top to bottom).

The yaw motion behavior is interesting. Thrust variations result in significant yaw moments which the vessel’s DP controller has to counteract. A normal distribution has been used to visualize the yaw station keeping performance of the vessel, see Fig. 14.

Fig. 14.

Normal distribution of yaw angles of proposed control model and DP model for simulation case 1, 2 and 3 (from top to bottom).

From Fig. 14 we find that the DP model has the highest probability density at a yaw angle of $0^{\circ }$ . This is as expected, since the DP model yaw angle setpoint is $0^{\circ }$ . Figure 14 also shows that the variance increases proportionally to the significant wave height. This is due to the fact that the vessel motions increase for higher environmental loads.

The 3DP model results show that the highest probability density is shifted slightly towards a yaw angle of $0.1^{\circ }$ . The yaw moments induced by the thrusters in the 3DP model are significantly higher compared to the DP model. The ${\mathit{K}}_{\mathit{i}}$ coefficient of the DP controller in the 3DP model should have been increased to remove the yaw offset. However, for a fair comparison, it was decided to use the same controller coefficients in both the DP and the 3DP model.

The plots show that the probability density function remains nearly constant in larger significant wave heights. This is because the yaw moments induced by the thrusters are significantly larger than the yaw moments resulting from waves, current and wind.

The yaw angle variance of the 3DP model is also 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 variance is $0.7^{\circ }$ , which is considered acceptable.

5.6. Station keeping capability

Next to the vessel DP footprint, also the station keeping capability is of interest. The station keeping capability is defined by the limiting current speeds at which the vessel is still able to maintain position. It is give by the positioning limits in DNV-GL [4] for DP capability level 3. Those limits are a maximum of 5 meter excursion and a maximum of $\pm 3^{\circ }$ yaw angle excursion.

In our investigation, the current velocity is increased per environmental condition until the vessel exceeds the positioning limits. The highest current velocities at which the vessel was still able to maintain position are visualized in Fig. 15.

Fig. 15.

Limiting current speeds per environmental direction for simulation case 1, 2 and 3.

From Fig. 15, we find that the vessel is able to maintain position in current velocities up to 2.3 m/s in head seas and 2.0 m/s following seas. In beam seas the station keeping capability is limited by a current velocity of 1.1 m/s, 1.0 m/s and 0.9 m/s for simulation case 1, 2 and case 3, respectively.

5.7. Power consumption

The cost of active roll reduction is expressed in terms of increased power consumption. 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 the 3DP and the DP control model is calculated and is expressed in a percentage of the total installed thruster power. The results are given in Fig. 16.

Fig. 16.

Mean power consumption percentage of total installed DP power for case 1, 2 and 3.

As expected, Fig. 16 indicates that the power consumption of the 3DP control system is higher compared to the DP system, and becomes higher for higher sea states. It was also to be expected 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 high thrust levels to counteract the roll moment. The power consumption of the 3DP model is around a factor 10 higher than the conventional DP model during brief critical moments in operations.

6. Conclusions

A control strategy for combined active roll reduction and DP station keeping in beam waves is presented. The performance of the strategy is investigated by conducting numerical analyses. We found that:

• The proposed control strategy is able to actively reduce the roll motion of the vessel.

• The effectiveness of the roll reduction reduces with increasing sea state.

• The 3DP control model is able to combine both active roll reduction and DP station keeping.

• The DP footprint of the vessel with active roll reduction shows a small increase with respect to the DP footprint of a conventional DP system. For the case with highest significant wave height under evaluation, the increase of the footprint was 1 meter, corresponding to 1% of the vessel length.

• The thruster power consumption increases with a factor 10 when the 3DP control model is engaged; this is considered acceptable for short periods of time during crucial parts of an offshore operation.

Active roll reduction in combination with station keeping works with the existing installed power. Our control model can be a viable instrument for extending the operability in special circumstances, such as a sudden change of weather, with no additional capital expenditure.