Proceedings of the IFAC Conference on

Proceedings of the IFAC Conference on Manoeuvring and Control of Marine Craft (MCMC'2012),
Genova, September 19-21, 2012
HOW TO INCORPORATE WIND, WAVES AND
OCEAN CURRENTS IN THE MARINE CRAFT
EQUATIONS OF MOTION
Thor I. Fossen
Department of Engineering Cybernetics,
Norwegian University of Science and Technology,
NO-7491 Trondheim, NORWAY
Abstract: This paper demysti…es how ocean currents together with wind and wave
loads in‡uence the marine craft equations of motion. In the literature there exists
great confusion of the use of absolute and relative velocity terms when modeling
rigid-body and hydrodynamic forces. The article is useful for engineers who want
to simulate and predict the motions of marine craft exposed to wind, wave and
ocean currents as well as control engineers evaluating the performance of marine
craft control systems. The results are also very useful for testing and tuning of
integral action time constants for compensation of ocean current and 2nd-order
wave-induced drift forces.
Keywords: Marine systems, wind, waves, ocean currents, equations of motion
INTRODUCTION
The equations of motion for underwater vehicles,
ships, ocean structures and high-speed craft are
usually derived using Newtonian and Lagrangian
mechanics (Fossen, 1994, 2011). The resulting
models are nonlinear mass-damper-springs which
include rigid-body, hydrostatic and hydrodynamic
generalized forces. The motions are coupled in six
degrees of freedom (DOF). Marine craft are exposed to environmental forces due to wind, waves
and ocean currents, which act like forcing on the
mass-damper-spring system. In hydrodynamics it
is common to assume linear superposition such
that forcing due to wind and waves can be treated
as generalized forces that can be directly added
to the nonlinear equations of motion. However,
generalized forces due to ocean currents do not
follow the law of linear superposition and there
exists much diversity and misunderstandings in
the existing literature on how to include the e¤ects
of ocean currents in the equations of motion.
This paper addresses the e¤ect of ocean currents
on marine craft in a tutorial perspective by using
the concept of relative velocity, that is the velocity
of the craft with respect to the ocean current, to
e¤ectively describe current-induced forces. Di¤erent properties and representations of the marine
craft equations of motion are discussed and it is
shown how the current velocity enters the equations. Nonlinear models for time-domain simulations and control systems design are presented
using the compact vectorial notation of Fossen
(1994, 2011).
1. MARINE CRAFT EQUATIONS OF
MOTION
In Fossen (1991) it was shown that the coupled
6-DOF equations of a marine craft could be expressed as:
_ = J( )
(1)
M _ + C( ) + D( ) + g( ) =
(2)
>
where = [N; E; D; ; ; ] and = [u; v; w; p; q; r]>
are the generalized position and velocity vectors,
respectively. The rotation and angular velocity
transformation matrices between body coordinates (BODY) and the North-East-Down (NED)
geographical reference frames are denoted Rnb and
T , respectively. The other quantities follow the
notation of (Fossen, 1994, 2011):
Rnb
Euler angle velocity
transformation matrix
System inertia matrix
M = MRB + MA
(including added mass)
Coriolis-centripetal matrix
C( )=CRB ( )+CA ( )
(including added mass)
Damping matrix
D( )
Vector of gravitational/
g( )
buoyancy forces
Vector of control inputs
J( ) =
06
6
06
T
6
The subscripts RB and A are used for the rigidbody and added mass terms, respectively. The
rigid-body system inertia matrix MRB satis…es:
_ RB = 06
M
MRB = M>
RB > 0;
MRB =
6
mS(rbg )
mI3 3
mS(rbg )
Ib
where
S( ) =
2
0
S> ( ) = 4
3
2
0
3
2
1
0
1
3
2
5;
=4
1
2
3
3
5
(3)
is the cross-product operator de…ned such that
a := S( )a.
The matrix CRB in (2) represents the Coriolis
b
vector term ! bb=n vb=n
and the centripetal vector
b
b
b
term ! b=n (! b=n rg ). Contrary to the representation of MRB , it is possible to …nd a large
number of representations for the matrix CRB :
We will present some useful representations below.
Theorem 1. (Coriolis-Centripetal Matrix from M).
Consider the 6 6 constant system inertia matrix:
M = M> =
M11 M12
M21 M22
>0
(4)
where M21 = M>
12 . Then the Coriolis-centripetal
matrix can always be parameterized such that
C( ) = C> ( ) by choosing:
C( ) =
03
S(M11
1
3
+ M12
S(M11
S(M21
with
1
:= [u; v; w]> and
2
2)
+ M12
1 + M22
1
2)
2)
:= [p; q; r]> .
PROOF. Sagatun and Fossen (1991).
(5)
The rigid-body Coriolis and centripetal matrix
CRB ( ) can always be parametrized such that it
is skew-symmetrical:
CRB ( ) =
C>
RB ( );
8 2 R6
(6)
The skew-symmetric property is very useful when
designing nonlinear motion control systems since
the quadratic form > CRB ( )
0: This is
exploited in energy-based designs where Lyapunov functions play a key role. There exists several parametrizations (Fossen and Fjellstad, 1995)
that satisfy (6).
Lagrangian parametrization: Application of Theorem 1 with M = MRB yields the following
expression:
CRB ( ) =
mS(
1)
03 3
mS(S(
b
2 )rg )
mS( 1 ) mS(S( 2 )rbg )
(7)
mS(S( 1 )rbg ) S(Ib 2 )
Velocity-independent parametrization: By using
the cross-product property S( 1 ) 2 = S( 2 ) 1 ,
f12g
f11g
it is possible to move S( 1 ) 2 from CRB to CRB
in (7). This gives an expression for CRB ( ) that
is independent of linear velocity 1 :
mS( 2 )
mS(rbg )S( 2 )
CRB ( ) =
mS( 2 )S(rbg )
S(Ib 2 )
(8)
Remark 1. Formula (8) is useful when irrotational
ocean currents enter the equations of motion since
CRB ( ) does not depend on the linear velocity 1
(uses only angular velocity 2 and lever arm rbg ).
This is discussed in Section 3.2.
Gravitational and buoyancy forces for surface
vessels show a linear behavior g( ) = G where
2
3
0 0 0
0
0 0
60 0 0
0
0 07
6
7
6 0 0 Zz 0
Z 07
6
7
G=6
(9)
K
0 07
60 0 0
7
4 0 0 Mz 0
5
M 0
0 0 0
0
0 0
For underwater vehicles
2
6
6
6
g( ) = 6
6
6
4
(W B) sin ( )
(W B) cos ( ) sin ( )
(W B) cos ( ) cos ( )
(y g W y b B) cos ( ) cos ( )
(z g W z b B) sin ( )
(xg W xb B) cos ( ) sin ( )
3
+ (z g W
+ (xg W
(y g W
7
7
z b B) cos ( ) sin ( ) 7
7 (10)
xb B) cos ( ) cos ( ) 5
y b B) sin ( )
which both are functions of the relative velocities:
W in d C o e ffic ien ts
urw = u
vrw = v
0 .8
0 .6
0 .8
0 .4
0
0 .4
- 0.2
0 .2
- 0.4
- 0.6
0
CX
0
30
60
90
120
150
Angle of wind γ (deg) relative bow
180
0
w
30
60
90
120
150
Angle of wind γ (deg) relative bow
180
w
0.15
0 .8
CN
0 .1
0 .6
CK
0.05
0 .4
0
0 .2
- 0.0 5
0
- 0.1
0
30
60
90
120
150
Angle of wind γ (deg) relative bow
180
(15)
(16)
CY
0 .6
0 .2
uw
vw
The nondimensional wind coe¢ cients CX ; CY ;
CZ ; CK ; CM and CN are usually computed using
h = 10 m as reference height while a is the air
density. The frontal and lateral projected wind
areas are denoted by AF w and ALw while Loa
is the length over all. The mean heights of the
areas AF w and ALw are denoted by HLw and HFw ,
respectively. Figure 1 shows four wind coe¢ cients
for a typical research vessel (Blendermann, 1994).
Wind coe¢ cients for other vessels are found in
Fossen (2011) and references therein.
2.2 Wave-induced forces and moments
0
w
30
60
90
120
150
Angle of wind γ (deg) relative bow
w
Fig. 1. Wind coe¢ cients CX ; CY ; CK and CN for
a research vessel.
where the gravitational and buoyancy forces act
through the centers of gravity (CG) and buoyancy
(CB ) de…ned by the vectors rbg := [xg ; yg ; zg ]>
and rbb := [xb ; yb ; zb ]> ; respectively.
180
The …rst- and second-order wave forces for varying
wave directions i and wave frequencies ! k are
fdofg
fdofg
denoted ~wave1 (! k ; i ) and ~wave2 (! k ; i ) where
dof 2 f1; 2; 3; 4; 5; 6g. The normalized force response amplitude operators (RAOs) are complex
variables given by (WAMIT Inc., 2010):
For control systems design it is common to assume
the principle of superposition when considering
wind and wave-induced forces such that (2) takes
the following form:
M _ + C( ) + D( ) + g( ) =
+
(11)
where wind 2 R6 and wave 2 R6 represent the
generalized forces due to wind and waves.
wind
+
wave
i)
fdofg
i)
~wave2 (! k ;
gA2k
=
ej\~w a v e 1 (!k ;
fdofg
i)
fdofg
i)
ej\~w a v e 2 (!k ;
The output from the hydrodynamic code is usually an ASCII …le containing RAOs in table format. Let us denote the imaginary and real parts
of the force RAOs by; Imwave1 fdofg(k; i) and
Rewave1 fdofg(k; i). The amplitudes and phases for
di¤erent frequencies ! k and wave directions i for
the …rst-order wave-induced forces can be computed according to the formulae:
fdofg
Fwave1 (! k ;
i)
= Imwave1 fdofg(k; i)2
+ Rewave1 fdofg(k; i)2
i)
2.1 Wind forces and moments
are functions of relative wind speed Vrw and angle
of attack rw according to:
p
2
(13)
Vrw = u2rw + vrw
atan2(vrw ; urw )
(14)
rw =
fdofg
~wave1 (! k ;
gAk
fdofg
\Fwave1 (! k ;
Wind can be de…ned as the movement of air
relative to the surface of the Earth. For a marine
craft moving at a forward speed, the wind forces
and moments:
2
3
CX ( rw )AF w
6 CY ( rw )ALw 7
6
7
6 CZ ( rw )AF w 7
1
2 6
7 (12)
V
wind =
7
2 a rw 6
6 CK ( rw )ALw HLw 7
4 CM ( rw )AF w HFw 5
CN ( rw )ALw Loa
i)
i) =
Fwave2 (! k ;
2. SUPERPOSITION OF WIND AND
WAVE-INDUCED FORCES
fdofg
fdofg
Fwave1 (! k ;
1=2
(17)
= atan2 (Imwave1 fdofg(k; i);
Rewave1 fdofg(k; i))
(18)
The amplitudes and phases for the second-order
mean forces are:
fdofg
Fwave2 (! k ;
fdofg
\Fwave2 (! k ;
i)
i)
= Rewave2 fdofg(k; i)
(19)
=0
(20)
Since the …rst- and second-order wave forces are
represented in terms of the complex variables
fdofg
fdofg
Fwave1 (! k ; i ) and Fwave2 (! k ; i ), the responses
for sinusoidal excitations can be computed using
di¤erent wave spectra. A frequently used family
of wave spectra is:
S(!) = A!
5
exp( B!
4
)
(21)
where di¤erent values for A and B are used. These
values depend on geographical location and wind
speed (see Chapter 8.2, Fossen 2011).
Fig. 2. Computation of …rst- and second-order
wave-induced forces from force RAOs.
When computing the wave-induced forces, linear
superposition is employed as illustrated in Fig. 2.
The relationship between the a wave spectrum
S(! k ) and the wave amplitude Ak for a wave
component k is (Faltinsen, 1990):
1 2
A = S(! k ) !
(22)
2 k
where ! is a constant di¤erence between the
frequencies. Let the wave-induced forces in 6 DOF
be denoted by the vectors:
h
i>
f1g
f2g
f6g
=
;
;
:::;
wave1
wave1
wave1
wave1
h
i>
f1g
f2g
f6g
;
;
:::;
wave2 =
wave2
wave2
wave2
For the no spreading case, the wave direction
constant such that:
fdofg
wave1
=
N
X
De…nition 1. (Irrotational ‡uid). The generalized
ocean current velocity of an irrotational ‡uid is:
c
)+
k
fdofg
g Fwave2 (! k ; )
k=1
k)
(24)
where
! e (U; ! k ; ) = ! k
= [uc ; vc ; wc ; 0; 0; 0]>
| {z }
! 2k
U cos( )
g
(27)
vcb
where vcb = [uc ; vc ; wc ]> is the linear velocity.
vcn = Rnb (
fdofg
\Fwave1 (! k ;
(26)
c
where c 2 R is the velocity of the ocean current
expressed in BODY.
fdofg
A2k cos (! e (U; ! k ; )t +
=
The ocean current linear velocity vector satis…es:
(23)
=
r
6
g Fwave1 (! k ; )
Ak cos ! e (U; ! k ; )t +
fdofg
wave2
The forces on a marine craft due to ocean currents
can be accounted for by replacing the generalized
velocity vector in the hydrodynamic terms with
relative velocities:
=
k=1
N
X
to turn the major currents to the East in the
northern hemisphere and West in the southern
hemisphere. Finally, the major ocean circulations
will also have a tidal component arising from planetary interactions and gravity. In coastal regions
and fjords, tidal components can reach very high
speeds, in fact speeds of 2–3 m/s or more have
been measured.
b
nb )vc
(28)
>
where
are the Euler angles
nb = [ ; ; ]
between BODY and NED, and Rnb ( nb ) 2 SO(3)
is the corresponding rotation matrix.
De…nition 2. (Irrotational constant ocean current).
An irrotational constant ocean current in NED
satis…es:
_ nb ( nb )vcb + Rnb ( nb )v_ cb := 0
v_ cn = R
(29)
where
(25)
is the encounter frequency. The assumption that
= constant can be relaxed to model spreading of
the main wave propagation direction (see Chapter
8.3, Fossen 2011).
3. EQUATIONS OF MOTION INCLUDING
OCEAN CURRENTS
Ocean currents are horizontal and vertical circulation systems of ocean waters produced by gravity,
wind friction and water density variation in different parts of the ocean. Besides wind-generated
currents, the heat exchange at the sea surface
together with salinity changes, develop an additional sea current component, usually referred to
as thermohaline currents.
The oceans are conveniently divided into two water spheres, the cold and warm water sphere. Since
the Earth is rotating, the Coriolis force will try
_ nb (
R
nb )
= Rnb (
b
nb )S(! b=n )
(30)
This implies that the ocean current linear velocity
vector in BODY coordinates is given by:
v_ cb =
S(! bb=n )vcb
(31)
3.1 Equations of Motion including Ocean Currents
In order to simulate irrotational ocean currents
and their e¤ect on marine craft motion, the following model can be applied:
MRB _ + CRB ( ) + g( )
{z
}
|
rigid-b o dy and hydrostatic term s
+ MA _r + CA ( r )
|
{z
r
+ D(
r) r
hydro dynam ic term s
=
where
r
=
wind
vb vbc
! bb=n
+
wave
}
+
(32)
(33)
is the relative velocity vector. Notice that the
rigid-body kinetics is independent of the ocean
current.
3.2 Equations of Relative Velocity
It is possible to simplify (32) by exploiting the
structure of CRB ( r ):
Theorem 2. If the Coriolis and centripetal matrix
CRB ( r ) is parametrized independent of linear
velocity 1 = [u; v; w]> , for instance by using
(8), and the ocean current is irrotational and
constant (De…nition 2), the rigid-body kinetics
satis…es (Hegrenæs, 2010):
MRB _ +CRB ( )
MRB _ r +CRB (
r) r
(34)
PROOF. Since the Coriolis and centripetal matrix represented by (8) is independent of linear velocity 1 = [u; v; w]> , it follows that CRB ( r ) =
CRB ( ). The property:
Fig. 3. Current speed Vc , current direction c and
current angle of attack c relative bow.
wind tunnels. The resulting forces are measured
on the model, which is restrained from moving.
(36)
In many textbooks and papers, for instance Blendermann (1994), wind and current coe¢ cients are
de…ned relative to the bow using a counter clockwise rotation c (see Figure 3). The current forces
on a marine craft at rest can be expressed in terms
of the area-based current coe¢ cients CX ; CY and
CN as:
1
AF c CX ( c )Vc2
(39)
Xcurrent =
2
1
Ycurrent =
ALc CY ( c )Vc2
(40)
2
1
Ncurrent =
ALc Loa CN ( c )Vc2
(41)
2
where Vc is the speed of the ocean current. The
frontal and lateral projected currents areas are
denoted AF c and ALc ; respectively while Loa is
the length over all and is the density of water.
Theorem 2 when applied to (32) gives the di¤erential equations:
For vehicles at rest and motions limited to surge,
sway and yaw, ocean currents are linearly superimposed according to:
MRB _ c + CRB (
r) c
=0
(35)
is proven by expanding the matrices MRB and
CRB ( r ); and corresponding acceleration and velocity vectors according to:
+
"
mI3 3
mS(rbg )
mS(rbg )
Ib
mS(! bb=n )
mS(rbg )S(! bb=n )
S(! bb=n )vcb
03 1
#
b
mS(! b=n )S(rbg )
vcb
b
03 1
S(Ib ! b=n )
=0
Finally, it follows that:
MRB _ + CRB ( ) = MRB [ _ r + _ c ]
+ CRB ( r )[ r + c ]
= MRB _ r + CRB ( r )
M _ r +C(
r)
r
r +D(
r +g( ) =
r)
+
vcn
0
_ = J( )
r
(37)
wind +
wave +
(38)
where M = MRB + MA and C( r ) = CRB ( r ) +
CA ( r ). Notice that only r and not
is used
in (38) if compared to (32). The model (37)–(38)
includes the bias v_ cn = 0 at the kinematic level
while (32) models drift due to ocean currents at
the kinetic level using r =
c.
3.3 Equations of Motion for Zero Speed
For low-speed applications such as DP, ocean
currents and damping can be modeled by three
current coe¢ cients CX ; CY and CN . These can
be experimentally obtained using scale models in
M _ + C( ) + D( ) + g( ) =
current + wind +
where
current
wave
+
= [Xcurrent ; Ycurrent ; Ncurrent ]> :
The current coe¢ cients can also be used at forward speed U > 0 and related to the surge resistance, cross-‡ow drag and the Munk moment
used in maneuvering theory by using the concept
of relative velocity (see Chapter 7.3, Fossen 2011).
4. OCEAN CURRENT SIMULATION
MODELS
Let the ocean current speed be denoted by Vc
while its direction relative to the moving craft is
expressed in terms of two angles: angle of attack
c and sideslip angle
c as shown in Figure 4.
For computer simulations the ocean current speed
and direction can be generated by using …rst-order
Gauss–Markov processes:
V_ c +
_c +
_ +
c
1 Vc
= w1
(42)
2 c
= w2
(43)
3 c
= w3
(44)
where wi (i = 1; 2; 3) are zero-mean Gaussian
white noise processes and i
0 (i = 1; 2; 3)
are constants. If 1 = 2 = 3 = 0, the models
reduce to a random walks, corresponding to time
integration of white noise. A saturating element
is usually used in the integration process to limit
the current speed to:
Vmin
Vc (t)
Vmax
(45)
The direction of the current can also be …xed by
specifying constant values for c and c :
3-D Irrotational Ocean Current Model: A 3-D
irrotational ocean current model is obtained by
transforming the ocean current speed Vc and directions ( c ; c ) from FLOW axes to NED velocities:
>
vcn = R>
y; c Rz;
2
2
3
Vc
4 0 5
c
0
3
Vc cos( c ) cos( c )
5
= 4 Vc sin( c )
Vc sin( c ) cos( c )
(46)
This expression can be transformed from NED to
BODY using the Euler angle rotation matrix Rnb .
Consequently,
2
3
2
3
uc
Vc cos( c ) cos( c )
>4
4 vc 5 = Rnb (
5 (49)
Vc sin( c )
nb )
wc
Vc sin( c ) cos( c )
2-D Irrotational Ocean Current Model: For motions in the horizontal plane, the 3-D equations
(49) reduce to:
c
= 0 and
uc = Vc cos(
c
)
(50)
vc = Vc sin(
c
)
(51)
= 0. Consequently,
p
Vc = u2c + vc2
A tutorial on how to include models for wind,
waves and ocean currents for marine craft has
been presented. The concept of equations of relative motions is used to model the e¤ect of ocean
currents while wind and wave-induced forces are
added under the assumption of linear superposition. The article is intended for control engineers
who want to simulate and predict the motions of
marine craft exposed to wind, wave and ocean
currents and use time-series to evaluate the performance of control systems.
REFERENCES
where the principal rotations Ry; c and Rz; c
are recognized as:
3
2
cos( ) 0 sin( )
0
1
0 5
(47)
Ry; = 4
sin( ) 0 cos( )
2
3
cos( ) sin( ) 0
(48)
Rz; = 4 sin( ) cos( ) 0 5
0
0
1
for
Fig. 4. Angle of attack and sideslip angle for a
marine craft.
5. CONCLUSIONS
=
(52)
W. Blendermann. Parameter Identi…cation of
Wind Loads on Ships. J. Wind Eng. Ind.
Aerodyn, JWEIA-51:339–351, 1994.
O. M. Faltinsen. Sea Loads on Ships and O¤ shore
Structures. Cambridge University Press, 1990.
T. I. Fossen. Nonlinear Modeling and Control
of Underwater Vehicles. PhD thesis, Department of Engineering Cybernetics, Norwegian
University of Science and Technology, Trondheim, Norway, June 1991.
T. I. Fossen. Guidance and Control of Ocean
Vehicles. John Wiley and Sons Ltd., 1994. ISBN
0-471-94113-1.
T. I. Fossen. Handbook of Marine Craft Hydrodynamics and Motion Control. John Wiley and
Sons Ltd., 2011. ISBN 978-1-1199-9149-6.
T. I. Fossen and O. E. Fjellstad. Nonlinear Modelling of Marine Vehicles in 6 Degrees of Freedom. International Journal of Mathematical
Modelling of Systems, JMMS-1(1):17–28, 1995.
Øyvind Hegrenæs. Autonomous Navigation for
Underwater Vehicles. PhD thesis, Dept. of
Engineering Cybernetics, Norwegian University
of Science and Technology, Trondheim, Norway,
2010.
S. I. Sagatun and T.I. Fossen. Lagrangian Formulation of Underwater Vehicles’Dynamics. In
Proceedings of the IEEE International Conference on Systems, Man and Cybernetics, pages
1029–1034, Charlottesville, VA, October 1991.
WAMIT Inc.
WAMIT User Manual.
<www.wamit.com>, 2010.