Guidance and Control of Ocean Vehicles-Thor I Fossen

MA = -

[

Xii. 0

o

00]

Yu Yr. Nu

Ni

Find an expression for the fluid kinetic energy TA and use Kirchhoffs' equations to derive a skew-sy=etric matrix CA (v) for this system. 2.10 Assume that M = M T > 0 and D(v) > O'r/ v'f O. Show that:

(a)

M~(TJ) = M;;(TJ) > 0 'r/ TJ E]R6

(b)

D~(v,TJ)

> 0 'r/ v

E ]R6,

TJ

E]R6

2.11 Compute Yu, N r, Yu and N, for a surface ship with main dimensions D = 8 (m) and L = 100 (m) at a wave circular frequenj;y w = LO (rad/s) by applying strip theory. 2.12 Derive (2.98) and (2.99) from (2.72) and (2.87) by using the formulas:

a x (b x c) (a x b) x c where S(a)S(b)

S(a)S(b)c S(S(a)b)c

=I S(S(a)b). Another useful formula is the Jacobi identity: a x (b x c)

+b x

(c x a) + c x (a x b) = 0

which can be expressed in terms of the skew-symmetric operator S(-) E 88(3) according to:

S(a)S(b)c +S(b)S(c)a + S(c)S(a)b = 0 Fina]]y show that (2.100) and (2.101) can be derived from (2.98) and (2.99). 2.13 Derive the nonlinear body-fixed vector representation for a marine vehicle moving in 6 DOF by applying the Quasi-Lagrangian approach. All terms should be expressed by matrices and vectors.

- - - - - - - - _ .. ---- - - - - - ---

Chapter 3 Environmental Disturbances

In the previous chapter a general model structure for marine vehicles was derived. In this chapter we will look further into details on the modeling aspects in terms of environmental disturbance models. Moreover, the following type of environmental disturbances will be considered: e

Waves (wind generated)

• Wind • Ocean currents In general these disturbances Will] be both additive and multiplicative to the dynamic equations of motion. However, in this chapter we will assume that the principle of superposition can be applied. For most marine control applications this is a good approximation.

3.1

The Principle of Superposition

The previous chapter has shown that the nonlinear dynamic equations of motion could be written:

In the analysis below it is convenient to write the damping matrix as a sum of the radiation-induced potential damping matrix Dp(v) and a viscous damping matrix Dy(v) = Ds(v) + Dw(v) + DM(v) containing the remaining damping terms" Hence we can write: D(v) = Dp(v)

+ Dy(v)

(3.2)

Based on this model we will apply the principle of superposition to derive the linear and llOnlinear equations of motion in terms of environmental disturbances.

Environmental Disturbances

58

Linear Equations of Motion Linearization of the Coriolis and centripetal forces C RB(V)V and C A(V)V about zero angular velocity (p = q = T = 0) implies that the Coriolis and centripetal terms can be removed from the above expressions, that is C RB(V)V = C A(V)V = 0. If we also linearize D(v)v about zero angular velocity, and u = uo, v = vo and w = wo, we can write (3.1) as:

[M RB +M A]v+ [N p where N

N

p

p,

+ N v ]V+ G1] =

(33)

TE +T

N v and G are three constant matrices given by:

= 8[D p (v)v]

8v

v=vo

,

8[D v (v)v] v N = 8v

:

8g(1]) G= 81]

(3.4)

1

Linear Equations of Motion Including the Environmental Disturbances Furthermore, the principle of superposition suggests that the environmental disturbances can be added to the right-hand side of (3.3) to yield: lYI RB

v+N v v +

lYI A v + N p v

~

--..".

+ G 1]

= T wave

'....

radiation-induced forces I

+ Twind + Tcurrent +T ...

(3.5)

I

environmental forces

In the previous chapter the mdiation-induced forces were referred to as subproblem one, Section 2.4. In this chapter sub-problem two is considered. Moreover, we want to find the forces on the body when the body is restrained from oscillating, and there are incident regular waves. These forces are recognized as the Fmude-KTiloff and diffmction forces. Generally, the forces of sub-problem two are computed by integrating the pressure induced by the undisturbed waves and the pressure created by the vehicle when the waves are reflected from the vehicle over the wet body surface (Faltinsen 1990). Since this procedure is mathematically involved and not to well suited to control systems design, we will restrict our treatment to the following approximate solution for the Froude-Kriloff and diffraction forces. Approximate Solution for the Froude-Kriloff and Diffraction Forces If the body is totally submerged, has a small volume and the whole body surface is wetted, a special solution to sub-problem two exists. By small volume we mean that a characteristic cross-sectional dimension of the body is small relative to the wavelength A. For a vertical cylinder small volume means that A > 5D, where D is the cylinder diameter. ROVs are usually within this limit. Let the fluid velocity vector be defined by Vc = tuc, vc, Wc, 0, 0, ojT where the last three fluid motion components are zero (assuming irrotational fluid) we can write (Faltinsen 1990):

I

1•

3.1 The Principle of Superposition

59

M FK VC + MA li c + N p 1/c.---.........--. + .N v 1/c (3.6) Froude-Kriloff diffraction forces viscous forces where M FK may be interpreted as the Froude-Kriloff inertia matrix, that is the inertia matrix of the displaced fluid. Moreover, let 'V be the volume of the displaced fluid and p the fluid density, hence the mass of the displaced fluid can be written: Tcurrent

=

'---..----'....

".

(3.7)

m=p'V The moments and products of the inertia of the displaced fluid are:

Ix = Jv (y2 + z2)pd'V l y = Jv (x 2 + z2)pd'V

I.

=

J"

l xy = Jv xypd'V Ix: = Iv xzpd'V l y • = J" yzpd'V

(x 2 + y2)pd'V

(3.8)

We can now establish the concept of displaced fluid inertia for a completely submerged body by defining the FK-inertia matrix M FK = M~K > 0 similar to MRB (see Section 2.3.1). Moreover,

MFJ(=

[ mI 3x3 mS(rB)

-m~(rB) ] 10

-

~

[. 0 0 0

ffiz B -mYB

0

m 0

0

0 0 ffi

-mZB

mYB

0

-ffiXB

ffiXB

0

mZB

-ffiZB

0

mYB

-ffiXB

I. -7:r:", -I=,

-I:r:v

Iv -7'1::'

-rn"0 ] mXE

-7=: -I!J;J:. I,

(3,9)

where rB = [XB, YB, zB]T is the center of buoyancy. Linear Equations of Relative Motion If we assume that M FK = M RB, that is the vehicle is neutrally buoyant and the mass is homogeneously distributed, the linear equations of motion can be combined to give: [MRB

...

..+ MAl. V,

M

+ [N p + N v ] 1/, + G 1] = "--'-..,...-----'

T wave

+ Twind + T

(3.10)

N

where 1/, = 1/ - 1/ c can be interpreted as the relative velocity vector In this case 1/c should contain the contribution from the currents. Nonlinear Equations of Relative Motion For underwater vehicles, an extension to the nonlinear case could be to write (Fossen 1991): JIII V,

+ C(1/,.) 1/, + D(v,) 1/, + g( 1]) = T wave + T wind + T

(3.11)

Envi ronm ental Distu rbanc es

60

and simple Often this appro ximat ion is also used for ships due to its intuiti ve In the next way of treati ng slowly-varying CUITents in terms of relative velocity, sections we will discuss mathe matic al models for vc, 'T wav~ and 'T wind.

3.2

Win d-G ene rate d Wa ves

ets appea ring The process of wave generation due to wind starts with small wavel allows short on the water surface. This increases the drag force which in turn break and waves to grow. These short waves continue to grow until they finally starts with their energy is dissipated. It is observed that a develo ping sea or storm frequency. A high frequencies creating a spectr um with peak at a relative high developed sea. storm which has been blowing for a long time is said to create a fully swell is being After the wind has stopp ed blowing, low frequency decaying sea or frequencyl. If formed. These long waves form a wave spect rum with a low peak , a wave the swell from one storm intera cts with the waves from anoth er storm we will only spectr um with two peak frequencies may be observed. For simplicity -generated consider wave spectr a with one peak frequency, see Figure 3L Wind components. waves are usually represented as a sum of a large numb er of wave wave spectr al The wave ampli tude Ai of wave comp onent i is relate d to the density furlction S(Wi) as (Newman 1977): , (3,12) Af = '2S(Wi) /::,w nt differwhere Wi is the wave frequency of wave comp onent i and /::,w is a consta ence between successive frequencies, 5(0))

,

~I chosen as a random Figur e 3.1: Figure showing wave spectrum with one peak Wi is frequency in the frequency interval Aw.,

k . Hence, Let the wave numbe T of one single wave comp onent be denoted by i

k, = 271"

,

)"

,

(3.13)

cy, lThe peak frequency of the wave spectrum is often referred to as the modal fr~quen

3.2 Wind-Generated Waves

61

where Ai is the wave lengtp, see Figure 3.2. The wave elevation ((x, t) of a longcrested inegular sea propagating along the positive x-a.xis can be written as a sum of wave components (Newman 1977): '. 1; (x,t)

H t

1 1

, ;.

T

1; (X,t)

H x

Figure 3.2: Characteristics of a wave traveling with speed c = AfT = w/k. In the figure A = wave length, H = wave height, A = wave amplitude, T = wave period and ( = wave elevation.

N

((x, t) -

:L A; cos (Wit -

kix

+
i=l

~ 1 ki Ai + L.."2 0

cos 2(Wi t - kix

+
(3.14)

i=l

where
. w; = k g tanh(k d) i

i

(3 . 15)

Here d is used to denote the water depth. This relationship is often referred to as the dispersion relation. For infinite water depth, that is dI Ai > 1/2, the dispersion relation reduces to = ki g since tanh(ki d) -> 1 as d/ Ai -> co. Unfortunately, Expression (3.14) repeats itself after a time 27f/ 6.uJ This suggests that a large number of wave components should be used, typically N = 1000. However, this problem can be circumvented by simply choosing Wi as a random frequency in the frequency interval b.w.

w;

Environmental Disturbances

62

Linear wave theory or Airy, theory represents a 1st-order approximation of the wave elevation ((x, t). This corresponds to the first term Ai in Formula (3.14). Furthermore, 2nd-order theory implies that an additional term 1/2 ki is included in the expression for (x, t). This is done by applying a so-called Stoke's expansion to solve the wave theory problem up to second order. Similarly, forces caused by these terms are usually referred to as 1st-or-deT and 2nd-or-derwave forces, respectively. 2nd-order theory is usually sufficient to describe the response of most marine vehicles in a seaway. 1st-order wave disturbances will describe the oscillatoric motion of the vehicle while the 2nd-order term represents the wave drift forces. The next section shows how to compute S(Wi) and thus Ai in (3.12).

Ar

3.2.1

,

·,i

Standard Wave Spectra

The earliest spectral formulation is due to Neumann (1952) who proposed the one-par'ameter- spectrum:

SeW) =

c w- 6

exp

(-2l w- 2 V- 2 )

(m 2s)

(3.16)

where C is an empirical constant, V is the wind speed and 9 is the acceleration of gravity. More recently, the tendency has gone towards another class of spectra motivated by the early work of Phillips (1958) who showed that the high frequency form of the sea spectrum was asymptotically limited by (w» 1): (3.17)

Sew)

where c< is a positive constant. We will concentmte our discussion on this type of spectra. Table 3.1: Description of wind, p. 162 of Price and Bishop (1974). Reproduced by permission of Chapman and Hall, Ltd. Beaufort number

o 1 2 3 4 5 6

7 8 9 10 11 12

Description of wind Calm Light air Light breeze Gentle breeze Moderate breeze Fresh breeze Strong breeze Moderate gale Fresh gale Strong gale Whole gale Storm Hurricane

Wind speed (knots) 0-1

2-3 4-7 8-11 12-16 17-21

22-27 28-33 34-40

41-48 49-56 57-65 More than 65

'i

I

3.2 Wind-Generated Waves

63

Bretschneider Spectrum ,

A more sophisticated spectrum than the Neumann spectrum has been proposed by Bretschneider (1959). The two-parameter Bretschneider spectrum is written:

1.25 W6 8(1.<1 ) = -4w5 H,2 exp (-1.25 (wo/w) 4)

(m 2s)

(318)

where Wo is the modal frequency and H, is the significant wave height2 (mean of the one-third highest waves). This spectrum was developed for the North Atlantic, for unidirectional seas, infinite depth, no swell and unlimited fetch. Pierson-Moskowitz Spectrum

Independently of this work Pierson and Moskowitz (1963) developed a wave spectral formulation for fully developed wind-generated seas hom analyses of wave spectra in the North Atlantic Ocean. The Pierson-Moskowitz (PM) spectrum is written: 5(0) (m 1s)

10,----,---,---,.--,----,---,----,---,__--,--, 9 8 7

6

5 4

3 2

ool-oJ.-,_-l_JLLL~HS~l~3~·m::::~:~~·~·~=~2I::;oE.8~~0~9~~ .

0,2

0.3

0.,4

0..5

0,,6

0,7

'

frequency (rndls)

Figure 3.3: Figure showing the PM-spectrum for different values of H,.

(319) where

l

(3.20)

(~r

(3.21)

A -

8.1· 10- 3

B -

0.74

'In some textbooks tbe significant wave height is denoted by HI!,'

Environmental Disturbances

64

Here V is the wind speed at a height of 1g.4 m over the sea SUIface and 9 is the gravity constant. By assuming that the waves can be represented by Gaussian random processes and that S(w) is narrow-banded, the PM,.spectmm can be reformulated in terms of significant wave height, that is:

A -

8.1· 10- 3

B _

0.0323

l

(3.22)

(-L) H,

2

= 3.11 H;

(3.23)

Table 3.2: Description of sea, p. 147 of Price and Bishop (1974). Reproduced by permission of Chapman and Hall, Ltd. Notice that the percentage probability for sea state code 0, 1 and 2 is summarized.

Sea state code 0 1 2 3 4 5 6 7 8 9

Description ofsea Calm (glassy) Calm (rippled) Smooth (wavelets) Slight Moderate Rongh Very rough High Very high Phenomenal

Percentage probability Northern North wide Atlantic North Atlantic

Wave height observed (m)

World

0 0-01 0.1-0.5 0.5-L25 L25-2.5 25-40 4.. 0-6.0 6.0-90 \ 9.0-140 1 Over 14.0

1L2486

83103

6.0616

3L6851 40.1944 12.8005 3.0253 09263 0.1190 0.0009

28.1996 42.0273 15.4435 4.2938 L4968 0.2263 0.0016

2L5683 409915 2L2383 7m01 2.. 6931 0.4346 0.0035 J

1

A description of significant wave height with percentage probability is given i n l Table 3.2. This implies that the wind speed V and significant wave height H, will be related through: ~

1 j

H,

=

V 0.21 -

j

2

(3 . 24)

9

This relationship is plotted in Figure 3.4. The modal fTequency Wo for the PMspectmm is found by requiring that:

(dZ) t;wo

= 0

(3.25)

\

I 1I I

•

Straightforward computation yields:

Wo To

J

-

t:

(3.26)

2 £ 1f (4B

(3.27)

i

II II

. I

3.2 Wind-Generated Waves

65

Hs (m)

25 '. 20 15 10 5 0 0

10

20

30

40

50

60 V (knots)

Figure 3.4: Significant wave height H, = 0.21 y2 / g (m) versus wind speed Y (knots) (1 knot = 0.51 m/s).

where To is the modal period. Substituting the values for A and B into (3 . 26), yields: Wo = 0.88

g

V = 0.40

V(9 If,

(3.28)

\

Hence, the ma.'dmum value of S(w) is:

5A Srn"Aw) = S(wo) = - - exp (-5/4) 4Bwo

(3.29)

The Bretschneider spectnim is described by two parameters H, and Wo and is t)lus referred to as a two-parameter spectnJm. Notice that if Wo is chosen as 0.40 H, the Bretschneider spectrum reduces to the one-parameter PMspectrum.

vg/

Wave Spectrum Moments The different wave spectra can be classified by means of so-called wave spectrum moments to illustrate some of the statistical properties of their parameterization. The spectrum moments are defined as:

rn.

=

f" w·

S(w) dw

(k = O...N)

(3 . 30)

For k = 0, we obtain: mo =

A S(UJ) dw = o 4B

l

OO

(3.31)

Environmental Disturbances

66

This simply states that the instantaneous wave elevation is Gaussian-distributed with zero mean and variance a 2 = A/4B. Hence, y'7iiO can be interpreted as the RMS-value of the spectrum. Furthermore, we obtain: '.

mr

-

m2

-

A

0306 B3j,j

V1rA

4-jB

(3.32) (3.33)

For the PM-spectrum the average wave period is defined as: (3.34)

while the aver'age zeTO-cTOssings period is defined as: (3.35)

Furthermore, by assuming that the wave height is Rayleigh distributed it can be shown that (Price and Bishop 1974): H, = 4y'7iiO

(3.36)

I

Modified Pierson-Moskowitz (MPM) Spectrum For prediction of responses of marine vehicles and offshore structures in open sea, the International Ship and Offshore Structures Congress, 2nd ISSC (1964), and the International Towing Tank Conference, 12th ITTC (1969b) and 15th ITTC (1978) have recommended the use of a modified version of the PM-spectrum, that is:

Sew)

411"3 H 2

= T;ws' exp

(-1611"3) T;w 4

(3.37)

This representation of the PM-spectrum has two parameters H, and Tz . Alternatively, we can substitute:

T z = 0.710 To = 0.921 Tr

(3.38)

if To and T r are more convenient to use. This representation of S(w) should only be used for a fully developed sea with infinite depth, no swell and unlimited fetch. For non-fully developed seas the following spectrum has been proposed by the ITTC.

I

I

.-I!'a

3.2 Wind-Generated Waves

67

JONSWAP SpecJrum In 1968 and 1969 an extensive measurement program w\lS carried out in the North Sea, between the island Sylt in Germany and Iceland. The measurement program is lmown as the Joint North Sea Wave Project (JONSWAP) and the results from these investigations have been adopted as an ITTC standard by the 17th ITTC (1984) Since the JONSWAP spectrum is used to describe non-fully developed seas, the spectral de'nsity function will be more peaked than the those for the fully developed spectra. The proposed spectral formulation is representative for wind-generated waves under the assumption of finite water depth and limited fetch The spectral density function is written:

S(w)

= 155 T4IH;u)5 exp

(-944) T4 4 (-y) y I

w

(3.39)

where Hasselmann et aL (1973) suggest that I = 3.3 and:

Y . _ -exp [(0.191WTI-1)2] M

(3.40)

for w:S; 5.. 24/T1 for W > 5.24/T1

(3.41)

y2a

where 007

a- { 009

This formulation can be used with other characteristic periods like To and Tz by substituting: I

T1 = 0.834 To

= 1.073 Tz

(3.42)

The peak value of the JONSWAP spectrum can be related to the PM-spectrum by the ratio: "(=

S(WO)JONSWAP

(3.43)

S(WO)PM

This value is usually between 1 and 7. Example 3.1 (Experimental Wave Spectrum Results) A full-scale experiment was performed west of Bergen in order to measure the motion components in surge, sway, heave, roll, pitch and yaw for a moving supply vessel (Far Scandia). Both time-series and power sp~ctral density functions P(w) are shown in Figure 3.5. The motion components corresponding to the figures are;

u v w

=

-

surge acceleration (m/s 2 ) sway acceleration (m/s 2 ) heave acceleration (m/s 2 )

t/J

=

f)

-

'if;

=

roll angle (deg) pitch angle (deg) yaw angle (deg)

68

Environmental Disturbances

u( t)

.~

0

•

~

i•

•

Fuu(w)

2

i ~

:1

150

~ j

100

j

!

1

I

50

..

-1

.,

-2 0

50

100

0 0.05

0.1

015

02

Hz

sec

v( t)

i ~

-~

F",,(w)

~

5

.~

600 400

0

200 -5 0

50

100

0 005

01

sec

0.15

0.2

0.15

0.2

Hz

w(t)

Fww(w)

4 300 2

~.

AJ

0

200 100

-2 -4 0

50

100

sec

01

Hz

F"",,(w)

q;(t)

X

40

I

10'

2.5

20

2

15

0 -20

-40 0

0

0.05

0..5

50

sec

100

0 0.05

0.15

0.1

Hz

0.2

3.2 Wind-Generated Waves

69

, e(t)

POO(W)

15

6000

...

10 5 0 ·5

4000

~

·10 0

I~ 50

sec

2000

100

0 0.05

0.1

0.15

02

Hz

P,p,p (w)

'f;( t) 10

600

5

400 ·5 200 ·10 ·15 0

50

sec

100

0

0.1

02

03

Hz

004

05

Figure 3.5: Experimental time-series and wave spectra in 6 DOF for a moving ship in sea state code 8 (measured with a Seatex MRU sensor unit) . Notice that the yaw signal is highly affected by the feedback signal from the autopilot. Reproduced by permission of ABB Industry in Oslo.

3.2.2

Linear Approximations to the Wave Spectra

A linear approximation to the PM spectral density function 5("1) can be found by writing the output y(s) from the wave model as:

y(s) = h(s) w(s)

(3.44)

where w(s) is a zero-mean Gaussian white noise process with power spectrum:'

Pww(w) = LO

(345)

and h( s) is a transfer function to be determined. Hence, the power spectral density (PSD) function for y(s) can be computed as:

Environmental Disturbances

70

The ultimate goal is to design an approximation Pyy(w) to Sew), for instance by means of linear regression, such that Pyy(w) reflects the energy distribution of S(w) in the actual frequency rangeo We will in the forthcoming discuss some linear approximations well suited for this purpose. 2nd-Order Wave Transfer Function Approximation

Linear wave model approximations are usually preferred by ship control systems engineers, owing to their simplicity and applicability. The first applications were reported in 1976 by Balchen, Jenssen and Srelid (1976) who proposed to model the high-frequency motion of a dynamically positioned ship in surge, sway and yaw by three harmonic oscillators without damping Later Srelid, Jenssen and Balchen (1983) introduced a damping term in the wave model to better fit the shape of the PM-spectrum. This model is written:

h( ) s =

Kw s

S2

+ 2 ( Wo s + wJ

(3.47)

where it is convenient to define the gain constant according to:

Kw = 2(woow

(3.48)

Here Ow is a constant describing the wave intensity, ( is a damping coefficient while Wo is the dominating wave frequency. Hence, substituting s = jw yields:

' )_ h(JW -

j 21((wo ow)w (wo - w2) + j 2 (wo w ?

(3.49)

This in turn implies that: Ih(jw) [ =

2 (( Wo

V(wJ From (3.46) we recall that:

W2)2

Ow)

w

+ 4 (( Wo W)2

..

2 Pyyw)= hJw [ = (WJ-W 2)2+4((wow)2 (

(3050)

[ (

.

) 2

4 (( Wo ow)2w

(

3.51

)

This expression for Pyy(w) is shown in Figure 3.6 From this it is seen that the maximum value of [Pyy(w)1 is obtained for w = Wo, that is:

State-Space Model

A lineal' state-space model can be obtained from (3A7) by transforming this expression to the time-domain by:

yet) + 2 (wo yet) + w6 yet) = Kw wet)

(3.53)

I

I

3.2 Wind-Generated Waves Defining Xhl = Xh2 and Xh2 = space model can be written:

71

as state variables, this implies that the state-

Yh

XH

-

AHxH+EHwH

(3.54)

YH

-

CHXH

(355)

where W H is a zero-mean white noise process. Writing this expression in component form, yields: (3.56)

Yh

[0 1] [ :~~ ]

=

(3.57)

This model is highly applicable for control systems design due to its simplicity. Applications will be discussed in later chapters..

\

10

.... I

I

~)

5

o o

.

05

1.5

....

.

;

...

3.5

2 3 Frequency (rad/g)

45

4

5

1h(jw)'

20

'" "0

"ca"

]\:(1)

...

;

0

0

-20 -40 10-2

..

..

i

,,

10-'

.. i·.·

..

"'--.

.

' ","-

100

.

. ''-.

. 10'

10 2

Frequency (rad/g)

Figure 3.6: Power spectral density Pyy(w) and amplitude Ih(jw)1 as a function of frequency for the linear wave spectrum (wo = 1.0, ( = 0.1 and cr~ = 10).

-

Environmental Disturbances

72

Higher-Order Wave Transfer Function Approximations An alternative wave transfer function based on five parameters has been proposed by Grimble, Patton and Wise (1980a) and Fung and Grimble (1983) This model is written:

h(s) =

2 w K s S4 + als 3 + a2s2 + a3S + a4

(3 . 58)

where ai (i = LA) are four parameters. Hence, four differential equations are required to describe the wave modeL Moreover, 1 0 0 1 0 0 -a3 -a2

Yh = [0

(3.59)

0 1 0] [

~~~

(3.60)

],

XM

The number of unknown parameters can be reduced by assuming that the denominator can be factorized according to:

h(s) =

2 Kws (s2+2(wos+WJ)2

(3.61)

More recently, Triantafyllou, Bodson and Athans (1983) have shown by applying a rational approximation to the Bretschneider spectrum that a satisfactory approximation of the high-frequency ship motion can be obtained by using the transfer function: K

h(s) = (s2

2

+ 2 (:: s + WJ)3

(3 . 62)

which only has three unknown parameters (, Wo and Kw. The advantage of the higher order models to the simple 2nd-order system (347) is that they will represent a more precise approximation to the wave spectrum. The disadvantage, of course, is higher model complexity and often more parameters to determine. 3.2.3

Frequency of Encounter

For a ship moving with forward speed U, the wave frequency Wo will be modified according to: 2

.

w.(U, Wo,!}) = Wo - Wo U cos f3 9

where w~

=

k 9 (assuming deep water) and:

(3.63)

.~

;

.; i , i

3.2 Wind-Generated Waves

We

73

= encounter frequency (rad/s)

U

= wave frequency (lad/s) = acceleration of gravity (m/s 2 ) = total speed of ship (m/s)

fJ

= the angle between the heading and the direction of the wave (rad)

Wo

9

Notice that the encounter frequency can be negative for large values of U. The definition of the encounter angle fJ is shown in Figure 3 7 This suggests that the wave spectrum for a heading-controlled ship moving at speed U > 0 should be modified to incorporate the frequency of encounter. For instance, we can rewrite (3.47) as:

h( ) .s =

Kw s

S2

+ 2(

te',

S

(

+ w;

3 64

)

However, it should be noted that the wave frequency of a dynamic positioned ship can be perfectly described by uJ e = te'o since U is close to or equal to zero Beam sea

Quartering sea

Bow scn

p = 30°. Following sea

Head sea

Figure 3.7: Definition of ship's heading (encounter) angle (Reid et at 1984).

3.2.4

Wave-Induced Forces and Moments

In order to simulate the motion of ocean vehicles, in the presence of irregular waves we will consider the effect of Ist- and 2nd-order wave disturbances. Superposition in Terms of Ist- and 2nd-Order Wave Disturbances The responses of an ocean vehicle in a seaway are usually computed by applying the principle of superpQsiiion.. Assume that the 1st-order wave disturbances can

I 74

Environmental Disturbances

be described by the damped oscillator (3 . 54)-(357) . Alternatively, Y H can be computed by using the spectral densit.y function S(w). Furthermore, we assume t.hat 2nd-order wave drift forces in the X-, y- and z-directions can be modelled by t.hree slowly-varying parameters: d = [dJ, d 2 , d 3 jT' Hence, 1st-order wave disturbances (oscillatoric motion) 2nd-order wave drift

AH

YH

-

CHXH

d

-

Wd

XH

XH

+ EH wH (3.65)

Here Wd is a vect.or of zero mean Gaussian white noise processes . Moreover, the principle of superposition suggest.s that. the vehicle dynamics and the 2nd-order wave disturbances can be combined t.o yield:

M v

+ C(v) v + D(v) ij

v

+ 9(TI) = d + T

= J(TJ) v

(3 . 66) (3.67)

The measurement equation is modified to include the 1st-order wave induced motion, that is: I

(3.68)

I

where the low-frequency position and attitude components usually are given by Y L = TJ Notice that we have included 1the wave drift forces in the dynamic equation of motion (process noise) while the oscillatoric motion is added to the model output (colored measurement noise). In many practical operations like ship steering and positioning this simple model is sufficient A more intuitive and physical approach would be to model 1st-order wave forces and moments as process noise as well This can be done by applying the following model description

,J

Y=YL+YH

JI

"i i

1, I

1st-Order Wave Forces and Moment on a Block-Shaped Ship

j

Consider the expression (i(X, t) in (3.14) for the wave elevation The wave slope Si for wave component i is defined as:

Si (X, t )

t). = d(i(x, dx = Ai ki sm(wit -

kix

+ rpi) + O(A i2)

(3 . 69)

The wave elevation and wave slope can be expressed in terms of We for a moving ship. For simplicity, we assume that x = 0 and that higher order terms can be neglected. Hence,

(i(t) = (i(O, t)

:= Ai

COS(Weit

+
This implies that the wave slope-can be computed according to:

(3.70)

75

3.2 Wind -Gen erate d Wave s

5i(t)

= 5i(O, t) = A; ki sin(weit +
(3.71)

nent i Based Here Wei is the encou nter frequency corresponding to wave compo on these expressions we can derive the forces and moments:

(1970) makes induced by a regular sea on a block-shaped ship. To do this Zuidweg from water the following assumptions: (1) the forces and mome nts only result bed by the pressure acting on the wetted surface, (2) the wave field is not distur a fluctuating ship and (3) the influence of the waves is accounted for by assuming e itself is pressure distrib ution below the water surface, whereas the water surfac suggests assumed to be undist urbed . Moreover, the principle of superp osition that: M

v + O(v) v + D(v) v + 9(71) =

'wave

+,

(3.73)

length L, For a ship where the wetted part is a rectan gular parallelepiped with der wave breadt h B and draft T, we obtain the following formulas for the 1st-or disturbances (Kalls trom 1979): N

Xwave(t) =

L'p 9 BLTcos (3 Si(t) i=l N

Ywave(t)

-

L -p 9 BLTsin(3 Si(t)

i=l N

Nwave(t)

=

2 pg BL(L L.2.. i;1 24

B 2) sin 2(3 S[(t)

(3.74)

(3.75)

(3.76)

if the ship where fJ is defined in Figure 37. These equations will only hold hull can the is small compared to the wavelength and the water surface across removing the be approximated as a plane surface. An altern ative approach, that (k i L), assumption of the plane water surface, can be derived by assuming and X (k i B) and (k i T) are smalL This results in the same expressions for wave Ywave while N wave is modified to:

Nwave(t) =

:f, 214 pg BL(L

2

-

B 2 ) T k[ sin 2fJ (i(t)

(377)

1=1

eg (1970). More detailed analyses of wave forces and mome nts are found in Zuidw algorithm to In order to imple ment the above formulas we can use the following compute the wave elevation (i and slope Si-

76

Environmental Disturbances

Algorithm 3.1 (Wave Elevation and Wave Slope) L Divide the spectral density function Sew) into N intervals with length !:"w, see Figtrre 3, L 2, Pick a random frequency Wi in each of the frequency intervals and compute S(Wi)" 3, Compute the wave amplitude Ai = !2S(Wi) !:"w and the wave number ki = for (i = LN), 4, Compute

3.3

Si and

wf/g

by applying Formulas (3,70) and (3,71),

Si

Wind

Wind forces and moments on a vessel can usually be described in terms of a mean wind speed in combination with a wind spectrum describing the variation of the wind speed (gusting), We will first describe some standard wind spectra for this purpose and then relate the wind speed and direction to the forces and moments acting on the vehicle,

3.3.1

Standard Wind Spectra

One the most used spectral formulations for wind gust is the Davenport (1961) spectrum:

S

- k w(w) [1

i 916700w

+ (191 wjVw(10)J2l 4 / 3

(3,78)

where

w

= = =

0,05 (turbulence factor) average wind speed at a level of 10 m above the water surface (knots) frequency of the wind oscillations (rad/s)

Another attractive spectral formulation is the so-called Harrls (1971) spectrum which is written: S ( )- k

ww

-

[1

5286 Vw(lO)

+ (286wjVw(10))2]5/6

( 3,,79)

These spectra are based on land-based measurements More recently Ochi and Shin (1988) presented a spectral formulation relying on wind speed measurements carried out at sea, This spectrum is written in non-dimensional form according to:

~1

583 f.

sU.)

for

420 ~O,70

(1+/F') 11 838

f.

,

(1+ 12,3 )".6

o::; f.

< 0003

for 0,003 ::; for

f. > 0,1

f. ::; 0,1

(3,80)

, }

77

3.3 Wind where

f.

-

8(f.) =

f

=

OlO 8(f)

-

10 f

/vw (10)

f· 8(f)/O lO' V~(1.o) frequency of oscillation (Hz) surface drag coefficient, see Ochi and Shin (1988) spectral density

(1972), Simiu Other useful spectr al formulations are Hino (1971), Kaima l et aL pp. 15-18 and and Leigh (1983) and Karee m (1985); see the 10th ISSe (1988) references therein. Linea r Appro ximat ion to the Harris Spect rum der approxThe above wind spectr a are nonlin ear appro ximat ions. A linear 1st-or imatio n for the Harris spectr um is:

h(s) =

K 1+

T

(3 . 81)

s

which implies that:

Sw(w)

~ IhCiwW = 1 +~~T)2

(3.82)

Hence, we can choose the time and gain consta nt accord ing to:

T = J286/ Vw (10)

(3.83)

Wind Veloc ity Profil e e we can use In order to determ ine the local velocity z (m) above the sea surfac the bound ary-la yer profile (see Bretsc hneide r 1969):

Vw(z) = Vw(lO) . (z/10)'/7

(3.84)

e where Vw(lO) is the relative wind velocity 10 (m) above the sea surfac 3.3.2

Wind Force s and Mom ents

n a slowlyAs menti oned in the previous section the total wind speed will contai nent (wind varying compo nent (average wind speed) and a high-frequency compo e vessel are gust). The result ant wind forces and mome nt acting on a surfac angle 'YR (deg) usually defined in terms of relative wind speed VR (knots ) and according to:

1 'YR = tan- (vR/UR) VR = JUh + Vk where the compo nents of VR in the x- and y-dire ctions are:

(3.. 85)

78

Environmental Disturbances

UR

-

VR

-

VW COS(fR) - u +U c Vw sin('/R) - v + Vc

(3..86)

(3.87)

Here (u,v) and (U c,vc) are the ship and current velocity components while iR =

7/Jw - ,p is the angle of relative wind of the ship bow, see Figure 38

r-----

x

I1

1jf

j

r--

I

i1

v.

L.-

.

Y

Figure 3.8: Definition of wind speed and direction.

We can simulate time-series for Vw and component according to: Xl

-

X2

-

VW

Wl

-~

,1pw by adding a mean and a turbulent I

X3 (X2

Xl +X2

-Kwz)

X4

-

7/Jw -

W3

-~

(X4

-le W4)

X3 +X4

where Wi (i = LA) are zero-mean Gaussian white noise processes and T and K are the time and gain constants of the Harris spectrum, for instance. For most ships the wind gust cannot be compensated for by the control system since the dynamics of the ship is too slow compared with the gusts. However, slowly-varying wind forces can be fed forward to the controller by measuring the average wind speed and direction. This requires the wind force and moment coefficients to be known with sufficient accuracy. We will now describe two attractive methods for computation of the wind force and moment vector: (3.88) acting on a surface ship.

3.3 Wind

79

Wind Resistance of Merchant Ships (Isherwood 1972). Isherwood (1972) suggested that one write the wind forces (surge and sway) and moment (yaw) according to:

Xwind Ywind Nwind

-

~ OXh'R) Pw vi AT ~ Oy( IR) Pw vi A L "21 CN ( IR) Pw VR2 A L L

(N)

(3 . 89)

(N)

(3.90)

(Nm)

(3.91 )

where Cx and Cy are the force coefficients and CN is the moment coefficient, and where Pw is the density of air in kgjm 3 , AT and AI, are the transverse and lateral projected areas in m 2 and L is the overall length of the ship in m. Notice that VR is given in knots Based on these equations measured data were analyzed by multiple regression techniques in terms of the following 8 parameters:

= B = AL = AT = Ass = S = L

C M

= =

length overall beam lateral projected area transverse projected area lateral projected area of superstructure length of perimeter of lateral projection of model excluding waterline and slender bodies such as masts and ventilators distance from bow of centroid of lateral projected area number of distinct groups of masts or kingposts seen in lateral projection; kingposts close against the bridge front are not included

Moreover, Isherwood found that the data were best fitted to the following three equations:

Cx Cy CN

-

2AI, 2AT L S 0 Ao + Al V + A 2 B2 + A 3 B + A4 L + As L 2A L 2A L S 0 Bo + B I + B 2 -B2-T + B 3 -B + B4~L + B s-L 12 2A L 2AT L S C Co + 0 1 12 + O2 B2 + 0 3 B + C4 L + Cs L

+ A6 M

(3.92)

Ass

+ B6 -A

(3.93)

L

(3 . 94)

where A; and B i (i =:' 0....6) and OJ (j = 0...5) are tabulated below together with the residual standard errors (S.E.).

---

,

Environmental Disturbances

80

'I'able 3.3: Surge induced wind force parameters (Ishelwood 1972). ,. iR(deg) 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

Ao

Al

A2

A,

2.152 1.714 1.818 1.965 2.333 1.726 0.913 0457 0341 0.355 0.601 0651 0.564 -0.142 -0677 -0.723 -2.148 -2707 -2.529

-5.00 -3.33 -3,97 -4.81 -5.99 -6.54 -4.68 -2.88 -091

0.243 0145 0.211 0.243 0.247 0.189

-0.164 -0.121 -0143 -0.154 -0.190 -0173 -0.104 -0.068 -0.031

129 254 358 364 3.14 2.56 3.97 3.76

0047 0.069 0064 0.081 0126 0.128

-0.175 -0.174

A,

S.E. 0086 0104 0.033 0096 0.041 0117 0.115 0.042 0.109 0048 0348 0052 0.082 0482 0043 0.077 0.346 0.032 0.090 0.018 0.. 094 -0.247 -0020 -0.372 0096 -0.031 0.090 -0.582 -0.024 -0748 0100 -0.028 -0.700 0105 -0.. 032 0123 -0529 -0.475 -0.032 0.128 -0027 0.123 127 0.115 1.81 0.112 1.55 Mean S.E. 0.103

As

A.I

'!'able 3.4: Sway induced wind force parameters (Ishelwood 1972).. iR (deg) 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170

Bo 0.096 0176 0225 0.329 1.164 1.163 0.916 0844 0.889 0.799 0.797 0.996 1.014 0.784 0.536 0.251 0.125

BI

B2

0.22 0.71 138 1.82 1.26 0121 0.96 0.101 053 0.069 0.55 0.082 - 0.138 0.155 0.151 0.184 0191 0.166 0.. 176 0.106 0.046

B, 0.023 0.043

-0.029 -0.022 -0.012

B,

Bs

-029 -059 -0.242 -0.95 -0177 -0.88 -0.65 -0.54 -0.66 -055 -0.55 -0.212 -0.66 -0280 -0.69 -0.209 -0.53 -0163

Bo

S.E. 0.015 0.023 0.030 0.054 0055 0.049 0.. 047 0.046 0051 0.. 050 0049 0.047 0.34 0.44 0.051 0.060 0.38 0.055 0.27 0036 0.022 0.044 Mean S.E.

I

'I

, ;1

.,,

',"

~:';

·i" "<;

a

~?

3.3 Wind

81

Table 3.5: Yaw induced wind moment para.meters (Isherwood 1972). IR (deg)

10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170

Co 0.0596 0.1106 0.2258 0.2017 0.1759 0.1925 0.2133 0.1827 0.2627 0.2102 0.1567 0.0801 -O..Ol89 0.0256 0.0552 0.0881 0.0851

Cl 0061 0204 0245 0457 0573 0480 0315 0254

C2

C3

00067 0.0118 0.0115 00081 00053 -0.0195 -0.0258 -0.0311 -0.0488 -0.0422 -0.0381 -00306 -0.0122

00101 0.0100 0.0109 00091 0.0025

C..

.

C,

-0074 -0170 -0380 -0472 -0523 -0546 -0526 -0443 -0508 -0492 0.0335 -0457 0.0497 -0396 0.0740 -0420 0.1128 -0463 0.0889 -0476 00689 -0415 0.0366 -0.220 Mean S.E.

S.B. 00048 0.0074 Oal05 0.0137 0.0149 0.0133 O.Ol25 00123 0.0141 00146 00163 0.0179 00166 00162 0.0141 0.0105 0.0057 0.0127

Wind Resistance of Very Large Crude Carriers (OCIMF 1977).

Wind loads on very large crude carriers (VLCCs), that is vessels in the 150000 to 500 000 (dwt) class, can be cociputed by applying the following approach.

X wind Ywind N wind

-

;~

vi AT CYwC"!R) ;~ vi A L CNwC"!R) ;~ vi ALL

CXw(''/R)

(N)

(3 . 95)

(N)

(3.96)

(Nm)

(3.97)

Here the non-dimensional force and moment coefficients C Xw , CYw and C Nw aTe given as a function of iR in Figures 3.9-3.1L Pa is the density of air in kg/m 3 , see Appendb: F, while 7.6 is a conversion factor. For ships that are not too asymmetrical with respect to the xz- and yz-planes, we can approximate:

CXwC"!R) ~ CYwC"!R) ~ CNwC"!R) ~

Cxw COSC"!R) CYw sinC"!R) CNw sin(2iR)

(3.98) (3.99) (3.100)

Figures 3.9-3.11 indicate that CXw E {-La, -0.8}, CYw E {-ID, -07} and E {-0.. 2, -0.05}. However, the figures also indicate that these approximations should be used with care.

CNw

82

Envir onme ntal Distu rbanc es

CXW

O.B~-~-~-~~-~- ....-~- -~- ~"~ ; ....

.... : --

O.B

.

. , :"'l' ...

, "." 0.. 2

o -02 -0.4 -0.6

Ballas ted tanker. :,

"

,

FuIIy loaded t~nkel

-OB .1'-_-'-_...L.._~_.....,

o

M

~

W

00

_ _~_..L ..._~ _-'

e_--'

100

lM

1~

lW

100

'Yr (deg) Figur e 3.9: Longitudinal wind force coefficient Cxw as a functio n of relative wind angle of attack 7R (OCIM F 1977).

, -0.1 -02 -03

-0.4

'i. , .. -.E ~.

.\\

.. '

."

..

.

\

;

..\ : '

';'"

"~"

~

.~

~

~

~

,/ .

I

":

"/' .

I

: I

.: .\.

~ "

.:.E • \

, /j.

:

. ,,:

'"1"''-;'''

"

:

-0.7

-o..B

; ~

'\ ;Fully.; load~.: d tanke r : \

-0.5

:

i'"

,. ,.;-, ~"f' : "':":': r':':''':':'~'~

•

"

:.-

".

I

•

i";'

/, (deg) Figur e 3.10: Latera l wind force coefficient Cyw as a functio n of relative wind angle of attack 7R (OCIM F 1977).

3.3 Wind

83

. , '.

-~,

015

.

.

.

.

tan~er

Fully lbadecl ..

0.1

..~

"

'r"-'-'~""T'

,: ~'

:

0-,05

,:., o

BaUasted tanker ·0.05

o

20

40

~

80

100

lW

140

lW

180

" (deg)

Figure 3.11: Wind yaw moment coefficient of attack "YR (OCIMF 1977).

GN",

as a function of relative wind angle

The following example adopted from OCIMF (1977) illustrates how the wind forces and moment on a 280 000 (~wt) tanker with cylindrical bow configuration can be computed. ' Example 3.2 (Wind Load Calculations for a 280000 dwt Tanker) Consider a tanker in fully loaded condition with:

£=325 (m) From Figures 3.9-3.11 for a wind angle of 30 (deg) we obtain:

CXw = -0.73;

Cy ", = -0.31;

Assume that the wind speed at a 20 (m) elevation is V",(20) = 66 (knots). Hence, we can compute (see Equation (3.84)):

V",(10) = Vw (20)(10/20)1/7 = 60 (knots) which according to Formulas (3.95)-(3.97) with VR = V",(lO) results in.

../Ywind

Ywind Nwind

--

-

-073· (1224/76) .60 2 .1130'" -478 (kN) -031· (1224/7.6) , 60 2 .3160'" -568 (kN) 0032· (1223/7.6) .60 2 .3160 . 325", 19054 (kNm)

Environmental Disturbances

84

Here the density of air is taken to be Pw = 1.224 (kg/m 3 ) corresponding to 20° (e), see Appendix F. If the ship's heading is changed, the above calculations must be repeated for the new incident wind angle "(R.

o Another usefUl reference discussing wind resistance on large tankers in the 100 000 to 500 000 (dwt) class is Van Berlekom, Triigardh and Dellhag (1974). For medium sized ships of the order 600 to 50 000 (dwt) it is advised to consult Wagner (1967). Finally, an excellent reference for moored ships is De Kat and Wichers (1991).

3.4

Ocean Currents

Currents in the upper layers of the ocean are mainly generated by the atmospheric wind system over the sea surface; see pp. 3S-44 of 10th ISSC (19SS). Besides wind-genemted currents, the heat exchange at the sea surface together with the salinity changes develop an additional sea current component, usually referred to as thermohaline currents. This process also explains why varying water types are observed in different climatic regions. The oceans are conveniently divided into two water spheres, the cold and warm water sphere, which again are separated by the So C isotherm.. Since the earth is rotating, the Coriolis force will try to turn the major currents to the right in the northern hemisphere and opposite in the southern hemisphere. Finally, the major ocean circulations will also have a tidal I component arising from planetary interactions like gravity. In coastal regions and fjords tidal components can obtain very high speeds, in fact speeds of 2 to 3 m/s or more can be measured. A world map showing most major ocean surface currents is found in Defant (1961). 3.4.1

Current Velocity

The 10th ISSC (19SS) proposed that one write the surface current velocity Vc as a sum of the following velocity components: (3.101) where

Vi Viw If,

Vm V3et-up

Vd

= tidal component

component generated by local wind.. = component ge!1erated by nonlinear waves (Stokes drift). = component from major ocean circulation (e.g. Gulf Stream) = component due to set-up phenomena and storm surges . = local density driven current components governed by strong density jumps in the upper ocean. =

3.4 Ocean Currents

85

Tidal Component Let the vertical component z (m) be measured positive downwards.. Hence, the velocity profile of the tidal component can be written: ' V,(O) V, () z = { V,(O) 10glO 1 + d-io

(9' )

for 0::; z ::; d - 10 for d - 10 < z < d

(3.102)

Here V,(O) (m/s) is the surface speed of the tidal and d > 10 (m) is the water depth. Component Generated by Nonlinear Waves (Stokes Drift) As mentioned in Section 32, 2nd-order wave disturbances or so-called wave drift forces can be treated as an additional current component. The contribution to the surface drift (Stokes theory) resulting from the irrotational properties of the waves is written: N

l;.(z)

= I: ki Wi A~

N

exp(-2 ki z)

= I:

4 7f2 A~

i=l

i=l

)" '

1;,

exp( -4 7f z/ A;)

(3.103)

1

The derivation of this expression is found in Sarpkaya (1981), Component Generated by Local Wind The component generated by the local wind is written: v,w(z) = { V,w(O)

o

d~~z for

0::; z ::; do for do < z

(3,104)

Here do is the reference depth for the wind-generated current usually taken to be 50 (m), Collar (1986) has shown that V,w(O) can be approximated as: V,w(O) = 002 VlO

(3.105)

where VlO (m/s) is the wind velocity measured 10 (m) above sea level. 3.4.2

Current-Induced Forces and Moments

This section shows that the current-induced forces and moments can be included in the dynamic equations of motion by two methods, Both methods are based on the assumption that the equations of motion can be represented in terms of the relative velocity: (3,106) where Vc ities

----

=

tuc,

Vc> Wc>

0, 0,

Or is a vector of irrotational body-fixed current veloc-

Environmental Disturbances

86

Method 1: Section 2.L I has already shown that the earth-fixed linear velocity could be transformed to body-fixed linear velocities by applying the principal rotation matrices. Let the ear th-fixed current velocity vector be denoted by [u~, v~, w~]. Hence, we can compute the body-fixed components as:

(3.107)

where J 1 ( r/J, e, 1{;) =

c1{;cll s1{;clI [ -sll

-s1{;cq, + c1psllsq, C1PCq, + Sq,SIlS1P clIs,p

s1{;sq, + c1{;cq,sll ] -C1pSq, + SIIS1PCq, cllcq,

(3.108)

Let us assume that body-fixed CUIrent velocity is constant or at least slowlyvarying such that the following holds: (3109)

Hence, the nonlinear relative equations of motion (3.11) take the form:

IM

v + C(v

V

T)

r + D(v,)

V

T

+ 9('1]) =

T

I

\iJ = J(1]) vi

(3.110)

(3.111)

Notice that this model representation is based on the state variables (v, vc, 1]) with v,. = v - Vc Method 2: An alternative representation of the nonlinear equations of motion is obtained by defining (v r , v~, 1]) as the state variables. Moreover, from (3.11) we have that:

IM v + C(v,) T

V

T

+ D(v r) v, + 9(1]) =

T

I

(3.112)

Furthermore, we can write:

i]

= J(1]) v = J(1]) (v r + vc)

(3.113)

We recall that: (3.114)

where v~ = [u~,v~,w~,O,O,ojT Hence,

3.4 Ocean Currents

87

(3.. 115) Next, the kinematic equations can be modified to include the new state variable 1/, and a vector I/~ describing the earth-fixed current velocity, that is: (3.116)

Three-Dimensional Current Model (Submerged Body) If the vertical velocity profile Vz (z) is lmown, the average current velocity the draft of the vehicle can be computed as:

Vc = -1 foT Vz(z) dz

Vc

over

(3.117)

T.o where T is the hull draft. The earth-fixed fluid velocity components (u~, v!;, w~) can be related to Vc by defining two angles Cl< (angle of attack) and f3 (sideslip angle) describing the orientation of Vc about the y- and z-axis, respectively (see Figure 3.12).

Figure 3.12: Orientation of average current velocity with earth-fixed X, Y, and Z axes,

Using the results from Section 2.1.1, we can write:

[

:~ ] = Cy,a Cz,-~ [~ ]

(3 . 118)

where Ci,j is the transformation matrix defined in (2.8) and Vc is the average current velocity in the earth-fixed reference frame. Expanding this expression yields: E Uc E

Vc

wE c

-

Vc cos Cl< cos f3 Vc sin f3 Vc sin Cl< cos f3

(3 . 119) (3.. 120) (3.121)

88

Environmental Disturbances

Two-Dimensional Current Model (Sur'face Vessel) For the 2-D case, the earth-fixed cunent components can be described by two parameters only, that is average current speed Vc and direction of cunent 13. Consequently, the above 3-D expressions reduce to:

-

Vc cos 13 Vc sinfJ

(3.122) (3123)

Since we ale considering the horizontal motion of the vehicle, we can assume that both if; and e are zero which implies that (u c, vc) can be computed flom (3.107) as:

[

~: ] = [ - ~~~ ~ ~~~~ ] [ ~~ ]

Substituting the expressions for

u;; and v;; into (3.124) finally yields Vc

Uc

Vc

-

Vc

cos(fJ - 1ft) sin(fJ - 1ft)

(3.124) 3

,

(3.125) (3.126)

x

'----------~y

Figure 3.13: Definition of average velocity Vc and direction f3 of the CUlTent for a surface vesseL 'Here we have used the trigonometric formulas: (1) costa - b) =cosacosb+ sinasinb and (2) sin(a - b) = sinacosb - cosasin b.

3.4 Ocean Currents

89

Generation of Ocean Currents

For computer simulations the average CUIIent velocity can be generated by using a 1st-order Gauss-Markov Process For instance Ve(t) can be described by the following differential equation: V;,(t)

+ Ito v;,(t) = wet)

(3.127)

where wet) is a zero mean Gaussian white noise sequence and Ito ~ 0 is a constant. In many cases it is sufficient to choose Ito = 0 which simply corresponds to a random walk, that is time integration of white noise. For details on Gaussian processes see Gelb, Kasper, Jr., Nash, Jr, Price and Sutherland, JL (1988). This process must be limited such that Vmin :s; Ve(t) :s; Vma:< in order to simulate realistic ocean currents. The following algorithm utilizing Euler integration and a simple limiter can be used for this purpose: Algorithm 3.2 (Current Generator) 1. Initial value: vetO) = 0 5 (Vma:<

+ Vmin).

2. Euler Integration with sampling time h (see Appendix B.2):

3. Limiter: if (Ve(k

+ 1) > Vma.,)

aT

(Ve(k

+ 1) < Vmin)

then

4. k = k + 1, return to step 2

o Similar algorithms for aCt) and (3(t) can be used to simulate time-varying directions. We will now show how current disturbances and 1st-order wave disturbances can be included in the ship steering equations of motion. Example 3.3 (Augmented Model for Ship Steering) The linear sway-yaw dynamics of a ship with single screw propeller can be written in terms of (see Section 5.3.1): a,2

a22

1

0] [

0 0

VL -

rL

1/JL

v, ]

+[

, b ]

b2

lJ

(3 . 128)

0

where VL is the sway velocity, rL is the yaw rate, lPL is the heading angle, lJ is the rudder angle and V e is a parameter representing slowly-varying currents (see Equation 3.126). The subscript L is used to denote the low-frequency motion components. The high-frequency oscillatoric motion 7/JH of the waves can then be

~------------------------------_I

Environmental Disturbances

90

added to the model by simply augmenting (3.128) to the linear wave model (356) and (3.57), which yields: VL h

,pL

V,

=

,pH . {H

all a" 0 0 0 0

a12

a22 1 0 0 0

b, b,

+

0 0 0 0

6+

-all sin(,IJ -,pL -,pH) -a21 sin(,IJ -,pr -,pH) 0 -1'0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0

Kw 0

0 0 0 0

-w; 1

0 0 0 0 -2(w, 0

VL TL ,pL V, 1j; f! {f!

Ywind Nwind

[

:~

]+

0 0 0 0

(3 129)

where Wl and W2 are zero mean Gaussian white noise pr'ocesses and Y wind and N wind are two additional terms used to describe the wind force and moment in sway and yaw. For' this system the compass measurement equation is written: (3.130)

, I

where v is zero mean Gaussian white noise process This particular way of modeling the ship-wave interactions is attractive for' control systems design and state estimation.

o

3.5

Conclusions

In this chapter, we have discussed simple models for wave, wind and currentinduced forces and moments in terms of the 6 DOF marine vehicle equations of motion. This is done by applying spectral formulations of wind and waves whereas currents are modelled as random walks. The models discussed in this chapter are mainly intended for simulation and design of model-based control systems . Consequently, accurate prediction of marine vehicles in the presence of wind, waves and currents require more advanced modeling techniques.. A more detailed discussion on regular and irregular Airy (linear) wave theory can be found in Newman (1977) and Faltinsen (1990) while a detailed discussion on 2nd- and 5th-order Stokes theory is given by Sarpkaya (1981). In addition to this, it is recommended that one consults Ochi and Bales (1977) for a comparison and discussion of different wave spectra. Some standard references on wave, wind and current models with application to ship control are Zuidweg (1970), Kallstrom (1979), Jenssen (1980) and Blanke (1981); De Kat and Wichers (1991) and references therein is an excellent reference for moored ships.

I'

3.6 Exercises

3.6

91

Exercises

3.1 Plot the wave frequency w as a function of depth d for a wave with wave length A = 100 m. Show that if the ratio d/ >. is large enough, w win approach -Ikii, where 9 is the acceleration of gravity. This result should also be verified theoretically.

3.2 Show that: (3.131)

is the modal frequency for the Pierson-Moskowitz spectrum Find an analytical expression for the peak frequency of the JONSWAP and the modified Pierson-Moskowitz (MPM) spectrum. 3.3 Plot the spectral density function for the Harris spectrum together with a linear approximation of the same spectrum in a dB-Iog10(w) diagram for Vw (10) E {ID (knots), 20 (knots), 50 (knots) }. Is the linear approximation valid for the whole frequency range?

3.4 Compare the spectral formulations for the Harris, Davenport and, Ochi and Shin spectra by plotting Sew) versus w in a dB-Iog lO (w) diagram. Comment on the results 3.5 Plot the frequency of encounter as a function of negative and positive speeds U with fJ = 0 and wo = 06 (rad/s). Repeat the computation with U = 1 (m/s) and fJ E [00 3600). Comment on the results.

Chapter 4 Stability and Control of Underwater Vehicles

Conventional autopilot design based on linear theory starts with the assumption that the 6 DOF underwater vehicle equations of motion can be described as a linear model linearized around a point of equilibrium. This may be a rough approximation for many control applications Indeed, underwater vehicles performing coupled maneuvers at some speed are known to be highly nonlinear in their dynamics and kinematics" In such cases autopilots based on linear control theory can yield poor performance_ It is a common assumption that linear control design is much simpler than its nonlinear counterpart" However, exploiting the structure of the nonlinear equations of motion often yields a relatively simple and intuitive nonlinear autopilot design" This will clearly be shown in this chapter which emphasizes the following topics: e

Remotely oper ated vehicle (ROV) equations of motion

e

Stability of underwater vehicles

e

Conventional and nonlinear autopilot design of PID-type

e

Linear quadratic optimal autopilot design

.. Decoupled autopilot design .. Sliding mode control e

Feedback linearization

e

Nonlinear tracking

o

Adaptive autopilot design

This involves the design of automatic speed control systems, systems for dynamic positioning and tracking, as well as autopilot systems for automatic steering and depth controL ,

i

;

-I 1

----------------------

1

94

4.1

Stability and Control of Underwater Vehicles

ROV Equations of Motion

This section discusses different representations of the linear and nonlinear ROV equations of motion. 4.1.1

Thruster Model

Bilinear Thruster Model In the general case, the thruster force and moment vector will be a complicated function depending on the vehicle's velocity vector v E JR6 and the control variable n E JRP (p;::: 6). This relationship can be expressed as: (4.1) where b(·) is a nonlinear vector function. A 1st-order approximation of the developed thrust T and torque Q for a single-screw propeller can be derived from lift force calculations (see Blanke 1981). Let n denote the propeller revolution, D the propeller diameter, p the water density and Va the advance speed at the propeller (speed of the water going into the propeller). Hence the following expressions for the propeller thrust can be established: (4.2) where J o = Va/(nD) is the advance number and KT is the thrust coefficient; see Section 62.1 for details. In the general case KT will be a four quadrant nonlinear function as shown in Figure 4.1. KT n>O Va > 0

n>O Va <0

J. .,

n
n 0

"

v. KT

Figure 4.1: Four quadrant positive KT curve as a function of Jo For positive values of Jo experiments verify that KT is approximately linear in Jo while the results for negative Jo-values often show a nonlinear behavior; see Van Lammem et aL (1969) and Fossen (1991), pp 45-47.

.,

4.1 ROV Equations of Motion

95

In the general case the forward and bac.kward thrusts will be non-symmetricaL However, many ROV thruster systems are designed to give symmetrical thrust. Furthermore, KT usually shows linear behavior in J o such that the following approximation holds:

KT = 0<1 +0<2;;

(43)

where 0<1 and 0<2 are two constants given by the curves shown in Figure 4.L This implies that the thruster force (4.. 2) can be written:

T(n, Va) = 71nln

Inln + TlnlV.lnl Va

(4.4)

where 71nln > 0 and 71nlV. < 0 By using a similar approach it can be shown that the thruster torque can be written:

(45) Here Qlnln > 0 and QlnlV. < 0 are design parameters depending on propeller diameter, shape of the duct, water density etc

Figure 4.2: Schematic drawing of propeiler. The above coefficients will also depend on nand Va since (44) and (4,5) are only first-order approximations to a more general expression. However, experiments have shown that this dependency can be neglected for most practical conditions of operation, The advance speed Va is related through the speed of the vehicle V according to (see Figure 4,2): Va = (1- w) V

(4.6)

where w is the wake fraction number (typically: 0.1-04), Using the result, (44) implies that the propeller force developed by a single propeller can be described by the nonlinear function:

(4.7)

Stabi lity and Cont rol of Unde rwate r Vehic les

96

v
,

v>o

v=O

j l

./

--------~,.~.~~;;+'"=~-"":...--

(",/ !

v
n

,

v=o

v>o

n and speed of the Figur e 4.3: Propel ler force r as a functio n of propel ler revolu tion vehicle v.

where bl = 1]nln > 0, b2 = -1]nl V.(l- w) > multiv ariabl e case could be to write:

°

and v = V. An extension to the

Ir=B l u-B 2 (u)v!

(4.8)

and u E lR P is where B l and B 2 (u) are two matrices of appro priate dimensions a new contro l variable defined as:

(4,,9)

,i

Section 6 2.1. A simila r discussion on propeller forces can be made for ships, see

1 ,

Affine Thrus ter Mode l

I

Ui

=

Inil ni

(i = L.p)

ed by an In many practi cal applications the bilinear mcidel can be appro ximat ce, we can affine model , that is a system which is linear in its input, For instan appro ximat e (4.8) as: (410)

Ir=B ul by letting B

I i I

= B I and: (4.11)

Moreover, this Notice that for zero velocity, that is v = 0, this will always be true, j-th prope ller implies that the propeller force in the i-th DOF developed by the can be descri bed by: ri

= B ij Uj;

B ij

= Tjnjn

(4.12)

the theory In the rest of this chapter, we will mainl y consider affine systems since area active an of non-affine control systems is quite limite d. In fact, this is still I

--"

4.1 ROV Equations of Motion

97

of research However, regulation of non-affine systems will be discussed briefly in Section 4,3.3, Actuator Dynamics Most thruster systems are driven by small DC motors designed for underwater operating conditions The dynamic model of a speed-controlled DC motor can be written:

L adia dt

dn

2rrJm dt

-

R - 'a i a

-

KM i a

-

2rr

K

NI

n + Ua

(4.13) (4.14)

Q(n, Va )

where Da is the armature inductance, Ra is the armature resistance, U a is the armature voltage, KM is the motor torque constant, Jm is the moment of inertia of motor and thruster, n is the velocity of the motor in revolutions per second and Q(n, Va) is the load from the propeller defined in (4,5). i

Q (Nm)

,j ,~.

~prDpeJrer-

,~

!, 1

I

I

1

I

I

~-L~..,

: I I I I

1

- - -- - - --

Ua

n (rps)

(voltage):

> I

I I I 1

I 1

Figure 4.4: Propeller transfer functions Due to physical limitations of the DC motor, hard nonlinearities like actuator saturation, Coulomb friction, dead-zones and hysteresis should also be included in the complete model. Neglecting these effects, implies that we can apply Laplace's transformation to (4.13) and (4.14). Moreover,

(4.15)

n(s) = h".(s) ua(s) - hQ(s) Q(s) where s is the Laplace variable and: Kr h".(s) = (1 + T rs)(l + T2 s);

K 2 (1

+ T 3 s)

hQ(s) = (l + T rs)(l

+ T2 s)

(4.16)

Here K i (i = 1,2) are two gain constants and Ti (i = 1,2,3) are three time constants depending on the parameters in (4.13) and (4.14)

Stability and Control of Underwater Vehicles

98

Optimal Distribution of Propulsion and Control Forces For underwater vehicles where the control matrix B is non-square and P 2: n, that is there are equal or more control inputs than controllable DOF, it is possible to find an "optimal" distribution of control energy, for each DOF (Fossen and Sagatun 1991a), Consider the quadratic energy cost function: J=

~uTWu

(4,17)

2

which can be minimized subject to:

I II

(418)

r - Bu = 0

Here W is a positive definite matrix, usually diagonal, weighting the control energy, For underwater vehicles which have both control surfaces and thrusters, the elements in W should be selected such that using the control surfaces is much more inexpensive than using the thrusters, that is providing a means of saving battery energy" Define the Lagrangian: L(u, A) =

~UTWU + AT(r -

Bu)

I

(419)

i

where A denotes the Lagrange multipliers Hence, differentiating the Lagrangian L with respect to u yields:

i

8L 8u From this expression we obtain:

\

T

-=Wu-B A=O

(4,20)

\I

(4,21) By using the fact that: r = Bu = BW- 1BTA (4,22) and assuming that BW- 1B T is non-singular, we find the following optimal solution for the Lagrange multiplicators: ,

A = (BW-1BT)-lr

;.

(4,23)

1

1

,~

H,: ~,

(4,24)

which suggests that u can be computed as: u = Bt" r

,',.

,~,

Substituting this result into (4,21) yields the genemlized inverse:

IBt" = W- B T(BW- BT)-l I

"

(4,,25)

In the case when all inputs are equality weighted, that is W = I, (4,24) simplifies to:

:

4.1 ROV Equations of Motion

99

(426) This simplified result is known as the Moore-Penrose pseudo inverse" Notice that for the square case (p = n), Bt is simply equal to B- 1 4.1.2

Nonlinear ROV Equations of Motion

The nonlinear ROY equations of motion can be represented both in the hody-fixed and the earth-fixed reference frames. This has already been shown in Section 25" The body-fixed and earth-fixed vector representations are as follows: Body-Fixed Vector Representation

IM v+ C(v) v + D(v) v +g(ry) = rl

(4.27) (4" 28)

Earth-Fixed Vector Representation

IM~(1)) ij + C~(v, 1)) r, + D~(v, 1)) r, + g~(ry) = r

T

(1)) r

I

(4,"29)

The state vectors are v = tu, v, w, p, q, r]T and 1) = [x, y, z, if, 0, 1p]T The different matrices and their properties are discussed more closely in Section 2 5" 4.1.3

Linear ROV Equations of Motion

The linear equations of motion are obtained by linearization of the general expressions (4.27) and (4.28) about a time-varying reference trajectory or an equilibrium point, for instance:

vo(t) 1)o(t) -

[uo(t), vo(t), wo(t),Po(t), qo(t), To(tW [xo(t), yo(t), zo(t), ifo(t), Oo(t), 1/Jo(tW

(4,.30) (4",31)

6 DOF Perturbed Equations of Motion Let the perturbations from the reference trajectory lIo(t) and ryo(t) be described by the differentials:

.6.lI(t) = vet) - vo(t);

.6.ry(t) = 1)(t) -1)o(t);

M(t) = r(t) - TO(t) (4.32)

Introducing the following vector notation:

f Av) = C(v) v;

(4"33)

100

Stability and Control of Underwater Vehicles

implies that (4,27) can be linearized according to: M b.v

+

I

I

I

(4,34)

ryo + b.ry = J(TJo + b.TJ) [vo + b.v)

(4 . 35)

8tc(v) b.v + 8td(V) b.v + 8g(TJ) b.TJ = b.T 8v v o 8v v 0 8TJ TJ 0

Perturbating (4.28) yields:

Substituting ryo = J (TJo) Vo into this expression implies that:

b.ry = J(TJo + b.TJ) b.v + [J(TJo + b.7/) - J(TJo)] Vo

(4 . 36)

Linear theory implies that 2nd-order terms (b.1/i b.Vj "'" 0) can be neglected, Hence,

Ib.ry = J(TJo) b.v + J*(vo, TJo) b.TJ I

(4.37)

Here we have rearranged the last term in (4.36) according to:

[J(TJo + b.TJ) - J(TJo)] Vo £ J*(vo, TJo) b.TJ The following two special cases of (4 . 37) are particularly useful: (1) (2)

Vo = 0 Vo = TJo = 0

--> -->

(4.38)

b.ry = J(TJo) b.v b.ry = b.v

Linear Time-Varying ROV Equations of Motion

Defining

Xl

= b.v and

M

Xl

X2

= b.TJ, yields the following linear time-varying model:

+ C(t)

X2 = J(t)

Xl

Xl

+ D(t) Xl + G(t) X2

=

(4,39)

T

+ J*(t) X2

(440)

where

C(t) = 8 t c (v)1 8v Vo(t)

G(t) = 8g(TJ) I 8TJ TJo(t)

D(t) = 8td(V) I 8v Vo(t)

r(t) = J*(vo(t), TJo(t))

Defining

X

J(t) = J(TJo(t)) (441)

= [xi, x[]T and u = T, we obtain the following state-space model: -M-IG(t) ] [ J*(t)

Xl] X2

+

[

M-I] 0

u

(442)

.I

-'"

4.1 ROV Equations of Motion

101

which can be written in abbreviated form as: :i;

= A(t) x + B(t)

u

(4.43)

Linear Time-Invariant ROV Equations of Motion In many ROY applications it is reasonable to assume that the ROY is moving in the longitudinal plane with non-zero velocity components Uo and Wo in the x- and z-directions, respectively. Furthermore, let us assume that the steady-state linear and angular velocity components: Vo = Po = qo = TO = 0 and that the equilibrium point is defined by the zero roll and pitch angles, that is: ,po = eo = O. Hence, the time-varying matrices in (4.42) simplify to the following constant matrices:

M-

[

-xv

-Xp

m=G - Xrj

m-YiJ

-Yw

-m=c - Yp

-Yq

-Xw -Xp mZa - X q

-Yw

m-Zw

mYG - Zp

-mYG -Xf.

Xv Yv Zv Kv Mv Nv

Yu [ X.

~:

D = -

Mu

Nu

mYG - Zp

-Yoj m::z::c -Yf

-mXa - Zq

-1;:;v - 1<.;

-Zf

-1"= - K,

Xw Xp Yw Yp Zw Zp Kw Kp Mw Mp Nw Np

X, Yq Z, K, M, N,

-Xwuo

o

0

0

0

0

-(xoW - xBB)

(W

[~1 ~];

mXG

-Y"

-z;

-I;::r;: - Kf

-Iv;:. - M f

I;:. - NI'

C 12 C 22

] As!>vff\;l\j

tho.t

X"':lG-';C&~ 0

+ Xwwo + ZqWO -(XpuQ + Zpwo) XqUo

"«

]

"'1

0

f

B)

~

]

=

[

(4.44)

0 0

,po = eo = 0, J1

[ -~r2

Zw)wo

(zOW-zBB) -(VoW - YBB)

If we assume that ~'o = constant and matrix J takes the form:

=

+ (m -

-(m - Xli )uo

0

J

C=

0

[~ ~ ~ (z:~~-z:~) o

Z, K, M, N,

-(X"uo

X,:uo +Z,:wo [ -(XqUo + Z
1= -

Zri -1;:;y - K4 J y - A-fq -Iy ; - Mf

-mYG -Xi-

-m:tG -

Y, ] X,

+ ZfWO) o XpuQ + Zpwo

0

=

K.

-mza-Yp

0 XwuQ - (m - ZtiT)wo [ -X{,uQ - Ywwo

G

-xw

-xv

m-X.

the kinematic transformation

co'>!'o

-sin>!,o

sinwo

cos1/Jo

o

0

0] 0 1

(4.45)

whereas J' = O. Consequently, the linear time-invariant model can be written as: :i: -

Ax+Bu

(4.46)

~ (4.47)

(2. \02)

Stability and Control of Underwater Vehicles

102

where A and B are constant matrices. Notice that C will be zero if we require that Uo = Wo = 0 in addition to Vo = o.

4.2

Stability of Underwater Vehicles

Stability of an underwater vehicle can be defined as the ability of returning to an equilibrium state of motion after a disturbance without any corrective action, such as use of thruster power or control surfaces. Hence, maneuverability can be defined as the capability of the vehicle to carry out specific maneuvers. Excessive stability implies that the control effort will be excessive while a marginally stable vehicle is easy to control. Thus, a compromise between stability and maneuverability must be made. Furthermore, it makes sense to distinguish between controls-fixed and controls-free stability. The essential difference between these terms is that: o

Controls-fixed stability implies investigating the vehicle's stability when the control surfaces are fixed and when the thrust horn all the thrusters is constant.

o

Controls-free stability refers to the case when both the control surfaces and the thruster power are allowed to vary. This implies that the dynamics of the control system must also be considered in the stability analysis.

These terms will be described more closely in the next sections. 4.2.1

Open-Loop Stability :.t;

Open-loop (controls-fixed) stability analysis of marine vehicles concerns the problem of finding static stability criteria based on the hydrodynamic derivatives For linear models this is quite simple thanks to the well known techniques of Routh and Hurwitz. This section shows that an alternative approach based on Lyapunov's direct method can be applied in the nonlinear case. For marine vehicles the Lyapunov function V can be chosen to represent the system's total mechanical energy. Consider the Lyapunov function candidate

V(1],il) =

~ ilT

Mry(1]) il +

f g~(z)

dz

(4.48)

where Mry and gry are defined in Chapter 2. Here V can be interpreted as the sum of the kinetic and potential energy of the vehicle. Hence, zero energy corresponds to the equilibrium point 1] = 0 and il = O. Instability corresponds to a growth in mechanical energy while asymptotic stability ensures the convergence of mechanical energy to zero. Differentiating V with respect to time (assuming Mry = M~ > 0) yields:

11 =

il T [Mry(1])ij + g,,(1])J

+ ~ il T

Mry il

(4.49)

, •

4.2 Stability of Underwater Vehicles

103

Hence the expression for V can be rewritten as:

Applying the skew-symmetric property: iJT (!VI~

11 =

iJT [M~ij

-

2C~) iJ = 0 V iJ, yields

+ C~(v, rJ)iJ + g,/(1J)]

(451)

In controls-fixed stability analyses, the dynamics of the control inputs is neglected Hence, we simply consider the system:

+ C~(ll, 1J)iJ + D~(v, 1J)iJ + g~(1J) = 0 to the expression for 11, finally yields:

M,/(1J)ij Applying this equation

(4.52)

Theorem 4.1 (Controls-Fixed Stability) According to Lyapunov stab,i[ity theory, Appendix Cl, sufficient conditions for controls-fixed stability are:

(1) V > 0 for all iJ, 1J E JRn whereas iJ

#0

and 1J

# 0.

Hence: (4.54)

if and only if the inertia matrix: M>O

(4.55)

Notice that J-l(1J) is defined for all1J E JRn while J(1J) is undefined for 0= ±90°

(ii) V < 0 for all v E IRn if and only if the damping matrix: D(v) > 0 V v E IR.n

(iii) V

->

co as

1I1J1I

->

co and

lIiJlI

->

(4.56)

co. This is satisfied for (4.48).

o The first condition simply states that the inertia including hydrodynamic added mass must be strictly positive. For underwater vehicles we can assume constant added mass (independent of the wave frequency) which implies that !VI = 0 and M = M T > O. The second condition simply states that the system must be

Stability and Control of Underwater Vehicles

104

dissipative which is also true for an uncontrolled undisturbed ROV Moreover, energy should not be generated by the system itself. The uncontrolled system above is said to be autonomous since it does not explicitly depend on time t. Hence, we can apply Lyapunov's direct method to prove stability (see Appendix C.1.1), If tracking of a time-varying reference trajectory is of interest, the new dynamics associated with the tracking error will be non-autonomous. By non-autonomous we mean a system with state equation: :i: = j(x, t)

(4.57)

where the nonlinear function j(x, t) explicitly depends on time. In order to prove convergence or stability of this system non-autonomous theory must be applied. 4.2.2

Closed-Loop Tracking Control

In this section, it will be shown how Barbiilat's lemma can be used to derive a non-autonomous tracking control law" The design methodology is best illustrated by considering a simple example. Example 4.1 (Velocity Tracking Control) Assume that we want to contml the vehicle's linear and angular velocities. Let the error dynamics be denoted by v(t) = v(t) - Vd(t) where Vd(t) is the desired state vector. For' marine vehicles as well as mechanical systems in general, we can define a Lyapunov function candidate,;

V(V,t)=~VTMV

(4,58)

which can be interpreted as the ''pseudo-kinetic'' energy of the vehicle. Differentiating V with respect to time (assuming M = M T and M = 0) yields:

V= vT MD expression for V yields:

(4.59)

vT[r - MVd - C(V)Vd - D(V)Vd - g(1])] - vTD(v)v

(4.60)

Substituting (4.27) into the

V=

Here we have used the skew-symmetric p1'Operty: vTC(v)v suggests that the control law could be be selected as! ,:

=0

'if v E]Rn

This

(4.61)

where K d is a positive r'egulator gain matTix of appmpriate dimension, Hence:' (4.62) type of control action is usually referred to as the Slotine and Li algorithm in robotics (Slotine and Li 1987). However, in this case the special structure of the underwater vehicle dynamics is exploited in the design.

4.3 Conventional Autopilot Design

105

Notice that V ::; 0 implies that Vet) ::; V(O) V t ::::: 0, and therefore that v i.s bounded. This in-turn implies that if is bounded. Hence, If must be uniformly continuous Finally, application of Barbalat's lemma (see Appendix Cl) shows that V ..... 0 which implies that v ..... 0 as t ..... (Xl

o

4.3

Conventional Autopilot Design

This section starts with a brief review of PID-control design before we discuss extensions to nonlinear control theory. 4.3.1

Joy-Stick Control Systems Design

It is common to classify control systems into two fUndamental types, open-loop and closed-loop (feedback) control systems. Figure 45 shows an open-loop ROY autopilot system where the commanded feedforward force and moment vector re is generated by the RaY piloL The output from the joy-stick system u is computed by applying the generalized inverse of the input matrix. Openloop systems work satisfactory if the environmental disturbances are not too large and if the numerical expression for B is known with sufficient accuracy. Improved robustness and performance in the presence of environmental disturbances can be obtained by applying a closed-loop control system of PID-type (proportional, derivative and integral) instead, see Figure 45 In this case, the pilot joy-stick is used to generate the commanded position and attitude 'TIc (or alternatively linear and angular velocity). Closed-loop control requires that sen, sor/navigation data are available for feedback. In Figure 4.5 a reference pre-filter is included to smooth out the commanded input. This is done to avoicl saturation in the actuator as a result of large tracking errors caused by steps in the commanded input. For a second-order system, the reference pre-filter is usually chosen as:

Btv

(4.63)

where TJd E JR6 is the desired output from the pre-filter, 'TIc E JR6 is the commanded input, A = diag{ (I> .,., (6} is the desired damping ratios and n = diag{wr, ... , W6} is the desired natural frequencies. The design of the PID-control law for tracking of the desired state 'TId is the topic for the next section. 4.3.2

Multivariable PID-Control Design for Nonlinear Systems

Most existing RaY-systems use a series of single-input single-output (SISO) controllers of PID-type where each controller is designed for the control of one DOF. This implies that the control matrices K p , K d and K i in the PID-controllaw:

-

Stability and Control of Underwater Vehicles

106

i'····',···,··,',··,,·,········

environmental disturbances

from

pilot

r ~

, ; feedforward lSQ~lr9.1.§1~l~W

I--JI

.. g(.) , ;: ROY d . ji 1.~~~~~:

t;

"

•

!;

;l~2~.~]!!.~!!1.~~i5:~.

..

1,

AOpen·loop (feedforward) force/moment control

environmental disturbances

sensor noise

from pilot

BClosed·loop (feedback) position/attitude control

Figure 4.5: Open-loop and closed-loop ROV control systems design.

r.

.'

" T PID

= K p e(t) + K d e(t) + K i

t

e(-r) d-r

(4.64)

should be chosen positive and diagonaL Here e = TJd - TJ is the tracking error.. However, most ROV systems for offshore applications use only simple P- and PI-controllers for automatic heading and depth control since it is difficult to measure (estimate) the velocity vector v. A standard PID-control design can be improved by using the vehicle kinematics together with gravity compensation.. Moreover, we will show that perfect set-point regulation can be achieved in terms of Lyapunov stability theory if TpID is transformed according to: (4.65)

In addition, we will assume that the control input vector is related to the thruster

.'

'~

4.3 Conventional Autopilot Design

107

forces and moments according to (4.l0)c Hence, the inverse mapping: u

= Btv

(4.66)

T

where Btv is the generalized inverse (see Section 4.Ll), can be used to calculate the desired controls u In the next section we will also show that excellent performance can be obtained for the whole flight envelope by including the vehicle kinematics and restoring forces in the PID-control designc Moreover, it is not necessary to perform a gain-scheduling technique to counteract the time-varying behavior of the dynamics and kinematicsc However, precautions against saturation and integral wind-up should be made. This is illustrated in Figure 4c6 where the PID-controllaw of the EAVE-EAST vehicle at the University of New-Hampshire is shown. This design is performed under the assumptions (without loss of generality) 'TJd = constant, J('TJ) = I and g('TJ) = o.

Vehicle Dynamics

"::'

11

Reset function

Figure 4.6: The EAVE-EAST Proportional Integral Derivative Controller (Venkatachalam et al. 1985)

4.3.3

PID Set-Point Regulation in Terms of Lyapunov Stability

In this section, we will investigate the closed-loop dynamics of the control law (464) and (4.65) under the assumption that the desired state vector:

'TJd = constant

(4.67)

This control problem will be referred to as regulation as the opposite of tracking control which involves the design of a feedback and feedforward controller for tracking of a time-varying smooth reference trajectory 'TJAt). Consider an affine ROV-model:

Stability and Control of Underwater Vehicles

108

M

v + C(v) v + D(v) v + 9(17) = B

(4.68)

'U

where 17 E JR n, V E JRn and 'U E JRT Let us assume that the input matrix Band gravitational forces 9(17) are known whereas M', C and D are unknown. Hence, the following considerations may be done: PD-Control of Nonlinear Square System (r = n) Assume that B is invertible and let the control law be chosen as a PD-controllaw where the term g( 17) is included to compensate for gravity and buoyancy, that is:

'u

= B-

1

[J 7 (17) K

p

e - K

dv + 9(17)] I

(4 . 69)

Notice that (4 . 64) and (4.65) are equivalent to (4.69) if K d = JTKdJ > 0 and K i = O. This control law is motivated from time differentiation of a Lyapunov function candidate:

(4.70) which yields

11 = v 7 [Mv - J7(17) K p e] Here we have used the fact that this expression for 11, yields:

11 = v 7

[B u -

e=

(4 . 71)

-7] = -J(17) v. Substituting (4.68) into

D(v) v - 9(17) - JI(17)

K e] p

(472)

Notice that vI C(v) v = 0 for all v E JRn . From this it is seen that the proposed PD-controllaw with appropriate cboices of K p = K~ > 0 and K d > 0 ensures that:

.1

(4.73) This means that the power is dissipated passively by the damping matrix D and actively by the virtual damping matrix K d We now only have to cbeck that the system cannot get "stuck" at 11 equal to zero, whenever e # 0 . From (4.. 73) we see that 11 = 0 implies that v = O. Hence, (4.68) with (4.69) yields:

'/

(4.74)

Consequently v will be non~zero if e # 0 and 11 = 0 only if e = O. Therefore the system cannot get "stuck" and the system state vector 17 will always converge to 17d in view of V --t O. This result was first proven by Tagegaki and Arimoto (1981) who applied the result to robot manipulator controL However, nonlinear control of underwater

, 1,

4.3 Conventional Autopilot Design

109

vehicles in terms of Euler angle feedback was first discussed by Fossen and Sagatun (1991b)0 Later this work has been extended to quaternion feedback regulation in terms of vector quatemion, Euler rotation and Rodrigues parameter feedback by Fjellstad and Fossen (1994b)0 PID-Control of Nonlinear Square System (r = n) Arimoto and Miyazaki (1984) have shown that the results above can be generalized to include integral action, Let:

(475) denote the generalized momentum of the vehicle, Hence, it can be shown by time differentiation of a Lyapunov function candidate: M-I

1

V(x)=-x

[

2

a

~

T

aI

K

0

K

la] K p

i

i

aK

(4,76)

X

i

where a is small positive constant and

rp,

x =

17,

l

(4,77)

err) drf

that V .:; 0 and that 17 converges to 17d = constant. This is based on the assumption that the PID-control law is taken to be: (4,78) where K

p,

K

i

and K

d

are matrices satisfying:

K

d

K

i

Kp

(4079)

> M~ > 0

(4.80)

2 > Kd+-K i

(4.81)

a

where a is a small positive constant chosen so small that:

1 a -(1- a)K d - aM ~ + 2

2

6

81\!I~

2::;(T)i - T)id)-i=1 8T)i

>0

(40.82)

It should be noted that this solution only guarantees local stability in a limited region about the origin, For details on the proof see Arimoto and Miyazaki (1984).

Stability and Control of Underwater Vehicles

110

Overdetermined System (r > n) If we have more control inputs than states to be controlled, we showed that B- 1 could be replaced by the generalized inverse Btv defined in Section 4 1.1. Hence, it is straightforward to show that the above results are valid for the non-square case T > n,

Non-Affine Systems (r > n) For the regulation problem it is straightforward to extend the above results to non-affine systems where the control input is given by (see (4.8)):

r=B l u-B 2 (u)v

(4,83)

In fact the nonlinear control law:

tu = Bt [J (7]) K T

applied to (4.68) implies that

p

e -

v

K'd + 9(7])]

I

(4.84)

V can be written: (4.. 85)

where K'd > 0 must be chosen such that this e..xpression becomes negative. Notice that the additional coupling term B 2 ( u)v only contributes to the system damping if B 2 (u) > O. If B 2 (u) < 0 we must choose K'd > -B 2 (u) to ensure stability Perfect Collocation (r = m) In some cases, we can design an output feedback control law that overcomes the problem that all states must be measured. This design is based on the assumption that the number of inputs u E lH" are equal to the number of measured outputs yE lR m . To do this, we will apply passivity theory; see Appendix C.3 for details. This suggest that the plant and control system can be described by two blod
u

-

ROV

y

(passive)

feedback control law ---.,. (strictly passive)

,

{.

Figur'e 4.7: Passive and strictly passive block

"

"

For simplicity, we will assume that the ROV can represented by a linear model which is quite realistic if only positioning is of concern The model is:

!

'."

:~

4.3 Conventional Autopilot Design

111

Mil+Nv+GT/=Bu

(486)

T/=V

(4.87)

Here the inertia matrix M, Coriolis, centripetal and damping matrix N = G +D, gravitational matrix G and input matrix B are assumed to be constant (the earthfixed coordinate system is orientated such that J(T/) = I whenever the ROV is perfectly positioned). Hence, we can write the above system in state-space form as:

x= A x +B

u

(4.. 88)

where x = [v T , T/TJT and A and B are given in (4.47). Consider the Lyapunov function candidate: 1

V = 2" x T P x

(4 . 89)

v = ~xT(ATp+ PA)x + xTpBu

(4.90)

with P = pT > O. Hence

Let us assume that the sensors and actuators can be located such that:

y=Gx

(4.91)

where G is a constant lmown matrix defined by: (4.92) and P satisfies the Lyapunov equation: (4.93) with Q = QT 2:

o.

Hence: .

1

V=yTu--xTQX (4 . 94) 2 This is referred to as perfect collocation between the sensors and actuators This result is also known as the Kalman- Yakubovich lemma (see Appendix C.3) which is used to check if a system is positive real. For linear causal systems positive realness is equivalent with passivity. We now turn ou~ attention to the last block representing the output feedback control law. According to Appendix 0.3 a system is strictly passive if and only if there exists a scalar Cl< > 0 and some constant f3 such that:

(4.95)

Stability and Control of Underwater Vehicles

112

For a linear output feedback control law: (4.96) where H(s) = diag{hi(sn (i hi ( s) must satisfy:

= L.T) to be

Re{hi(s -

un

strictly passive the transfer functions

~ 0 'r/

w~ 0

(4.97)

for some u > 0 and: Lhi(jW) < 90° 'r/ w ~ 0

(4.98)

This is satisfied, for instance, if hi(s) is chosen as a PID-controllaw with limited derivative and integral action, that is:

f3 I+Ti s I+Td s (4.99) 1 + f3 Tis 1 + Cl< Tds < 1 and f3 > 1. Finally, Definition G8 ensures

h(s)-K ,

-

p

Here K p > 0, T i > Td, Cl< that y E LT. It should be noted that it is straightforward to generalize these results to a nonlinear ROV model by using the general framework of passivity. A related work on collocation is found in S!ilrensen (1993) who has applied this design methodology to control high-speed surface effect ships (see Section 7.1). 4.3.4

Linear Quadratic Optimal Control

Linear quadratic (LQ) optimal control design is based on minimization of a linear quadratic performance index representing the control objective. Consider the linear state-space model:

:i:: -

Ax+Bu+Ew

y

ex

-

(4.100) (4101)

where x is the state vector, u is the input vector, w is the disturbance vector and y is used to describe the control objective. Let J be a performance index weighting the tracking error vector against the control power, that is: min J=

~

(T(fJTQfJ

2 la

+ uTpu)

dr

(4102)

Here P > 0 and Q 2: 0 are the weighting matrices and fJ = y - Y d is the tracking error vector. The commanded input vector is denoted Yd' An approximate optimal solution to the tracking problem (4.102) for 0 « T <
(4.103) ,,

z

4.3 Conventional Autopilot Design

113

w

E

x

y

x

c

B

Figure 4.8: Linear Quadratic Optimal Autopilot A block diagram of the control system is shown in Figure 48 Under the assumption that Yd = constant and w = constant, the following steady-state solution is obtained; see Appendix D for details:

G1

_

_p-l B T R

G2

-

G3

-

-P-1BT(A+BG1)-TCTQ P-lBT(A+BGltTRooE

oo

(4.104) (4.105) (4.106)

Here R oo is the steady-state solution of the matrix algebraic Riccati equation:

Optimal state estimation (Kalman filtering) can be used to realize the autopilot in the case when not all states are measured. For instance, the LQG/LTR (Loop Transfer Recovery) design methodology have been applied to underwater vehicles by Milliken (1984) and Triantafyllou and Grosenbaugh (1991). Loop shaping techniques like the LQG/LTR design methodology allow the designer to deal with robustness issues in a systematic manner. Moreover, robust stability (RS) can be guaranteed if bounds on the uncertainties are known. On the contrary, robust performance (RP) is still an unsolved problem. A linear controller design can be checked for RP by performing a structured singular value analysis. This technique is often referred to as the M-analysis technique in the technical literature (see e.g. Maciejowski 1990). Nevertheless, the design of a so-called M-optimal controller is still an active area for resear ch.

114

4.4

Stability and Control of Underwater Vehicles

Decoupled Control Design

Healey and Marco (1992) suggest that the 6 DOF linear eq@tions of motion can be divided into three non-interacting (or lightly interacting) subsystems for speed control, steering and diving Each systems consists of the state variables: 1) Speed system state: u(t) . 2) Steering system states: v(t), r(t) and 1f;(t). 3) Diving system states: w(t), q(t), e(t) and z(t). The rolling mode, that is p(t) and q,(t) is left passive in this approach. This decomposition is motivated by the slender form of the Naval Postgraduate School (NPS) AUV (see Figure 4.9).

Figure 4.9: Schematic drawing of the NPS AUV II (Healey and Lienhard 1993) The mathematical model and specifications of the vehicle are given in Appendix E. 2.

Healey and Lienard (1993) have applied the theory of sliding regimes to control the NPS AUV IL This control system has been successfully implemented and tested at the NPS in Monterey A related work discussing the problems on adaptive sliding mode control in the dive plane is found in Cristi, Papoulias and Healey (1990) We will discuss this design methodology in a later section. The above configuration suggests that the three subsystems can be controlled by means of two single-screw propellers with revolution n(t), a rudder with deflection 8R (t) and a stern plane with deflection 8s (t). This particular choice of actuators is inspired by those used in flight and submarine control. Of course other combinations of control surfaces, thrusters etc. can be used to control the above subsystems. Nevertheless, we will use this simple actuator configuration to illustrate how a decoupled control design can be performed in terms of: (1) proportional, derivative and integral control and (2) sliding mode controL The AUV examples in this section are based on the NDRE-AUV (see Figure 4.10). This is a test vehicle designed by the Norwegian Defence Research Establishment where the main purpose has been to test a propulsion system using sea water batteries (Jalving and St0rkersen 1994)

'

,

---

4.4 Decoupled Control Design

115

Figure 4.10: Schematic drawing of the NDRE-AUV (Jalving and Stilrkersen 1994). Specifications: length of hull = 4.3 m, maximum hull diameter = 0.7 m, propeller diameter = 0.6 rn, cruise speed = 2.0 m/s and hull contour displacement = LO m3 . 4.4.1

Forward Speed Control

Neglecting the interactions from sway, heave, roll, pitch and yaw suggests that the speed equation can be written as: (4.108) X u) iL = X 1u1u lulu + (1 - t) T + X ext Here we have assumed that quadratic damping is the dominating dissipative effect. Furthermore n represents the propeller revolution, u is the surge velocity, X ext is external disturbances due to waves and currents and t is the thrust deduction number t. Recall that the thruster force T and moment Q can be written:

(m -

T

= 71nln Inln + TlnlVo InlVa;

Q

= Qlnln Inln + QlnlVo InlVa

For simplicity, we will assume that 71nlVo = 0 (affine system) notation Xlnlnl = (1 - t)71nln, finally yields:

I(m -

X u) iL = X 1u1u

lulu + X 1n1n Inln + X exl I

(4.109)

Introducing, the

(4.110)

We will now demonstrate how a speed control sY2tem can be designed for this modeL Inner Loop PI Control System The propeller revolution n( s) can be measured with a pulse counter Or a tachogenerator. Hence, an inner loop feedback control system can be designed by applying a PI-controller (see Figure 4.. 11). h

( )

ua(s) Kp (1 + Tis) = - nd(s) - n(s) Tis

'''1' S -

,

(4.111)

Here nd( s) is the desired propeller revolution and Ua (s) is the armature voltage (see Section 4.11). This implies that (4.15) can be written: (4112)

Stability and Control of Underwater Vehicles

116 ...... ............. " ...• " "

"

"

..

' "."""."

..

~

".... ..

.,

,.,

"

",

"

"

".

Outer loop velocity control system

Va

. .

"

. \

" • • •" "

~" • • . ,,

, • .. • • • • • • • • .. • ..• • • • • • • • • .. • • • .. • . . "

• ••••• H

•••• ""

"

" ••

Figure 4.11: Speed Control System where

(HI3) PIOper tuning of the PI-contIOl parameters will ensure that n( s) tracks the desired propeller revolution nd(s). The main advantages with an inner servo loop is that sensitivity to varying load conditions Q(s), nonlinear actuator dynamics and hysteresis are reduced. The desired forward speed Ud corresponding to the propeller revolution nd can be solved from (4110) under the assumption that X ext = 0 and that all model parameters are known. Unfortunately, it is quite obvious that this result will be uncertain.. If accurate speed is important, an outer loop control system must be design in addition to the inner servo loop. Outer Loop Velocity Control System If u is measured or at least estimated, an outer velocity servo loop can be designed for proper tracking of .the desired velocity Ud. This is illustrated in Figure 4.11, where the nonlinear term Inln is generated by a PI-control law under the assumption that n tracks nd perfectly, that is:

(4114)

,;

;~

i

1

4.4 Decoupled Control Design

117

Here ii = Ud - u is~_he tracking error and Cl > 0 and C2 > 0 are the regulator proportional and integral gains, respectively, Example 4.2 (The NDRE-AUV Speed Control System) The performance of the inner-loop speed controller, Equation (4111), has been demonstrated by NDRE 2 who has designed a long-range AUV for testing of new battery technology and advanced control theory, Figure 4,10. The NDRE-AUV has a low drag hull and a battery capable of delivering energy jor a long mission In May 1993 the NDRE-A UV was successfully tested in the open sea between Norway and Denmark It then traveled a distance of 109 nautical miles (rpm)

, 34r---,.-----,---,----,-----,-----,--~--_,_--,

'22 '1200! ;----'OCO'-.S---;-----".-',.,-S--"'2----,2",S,,----*---,;,,""',S---!4 (hour.,)

Figure 4.12: Desired and actual propeller revolution versus time for the NDRE-AUV. The propeller revolution n for a typical mission is shown in Figure 412 where the desired propeller revolution is held constant at nd = 127 rpm,

o 4.4.2

Automatic Steering

Automatic steering or heading control can be done by means of a rudder or a pair of thrusters In the next section we will show how a rudder control system can be designed for course-changing maneuvers and then illustrate the performance of the control system by considering a small example Steering Equations of Motion According to (4.44), the linear steering equations of motion can be expressed in a compact form as:

[

m - Yu mxe - Nu

o

mxe - Y; I, - N r 0

2The Norwegian Defence Research Establishment, Kjeller, Norway

Stability and Control of Underwater Vehicles

118

where v. is the sway velocity, T is the angular velocity in yaw, 'ljJ is the heading angle and oRis the rudder deflection. ReaITanging this expression into state-space forID, yields: = Y =

:i;

Ax+boR eT x

(4.116)

(4.117)

where x = [v, T, ,p]T and y =,p. Moreover, (4.118)

where the choices of aij and bi should be quite obvious. Consequently, the transfer function between ,p and oRis obtained as:

:!L(s) = eT(sJ _ A)-lb = (a21 bl - allb2) + b2s OR S[S2 - (all + a22)S + alla22 - a12a21]

(4.119)

or equivalently:

:!L(s = K(l + Tas) OR) s(1+T1 s)(1+T2 s) where K is a gain constant and T j (j

=

(4.120)

1,2,3) are three time constants.

Autopilot Design The heading contIOI system can be designed by applying a PID-controllaw:

OR(S) = K~ 1 + Tis 1 + Tds [,pd(S) _ ,p(s)] Tis 1 + Tfs

(4.121)

where K p (controller gain), Ti (integral time constant), Td (derivative time constant) and Tf ~ 0,1 Td (low-pass filter time constant). The loop transfer function is:

I(s) = KKp (1 + 1is)(l + Tds)(l + Tas) TiS 2 (1 + Tf s)(l + T 1 s)(1 + T2 s)

(4.,122)

Hence, the closed-loop dynamics is described by:

I(s) ,p(s) = 1 + I(s) 7j;d(S)

(4.123)

It is seen from the final value theorem that the yaw angle will converge to the desired value for a step response: ,pd(S) = ,po/s, that is: lim ,p(t) = limslj;(s) = lim

t-lOO,_o

5-0

1(~~),pO =,po s

1+

(4,124)

404 Decoupled Control Design

119

The yaw angle can be measured by a compass while rate measurements usually are obtained by a rate gyro or a rate sensor, If the compass measurements are of good quality, rate estimates can be obtained from numerical differentiation or state estimation (Kalman filtering), Example 4.3 (The NDRE-AUV Heading Control System) A simpl~fied version of (4·121) without integral action has shown to perform satisfactory for the NDRE-A UV This control law is simply taken to be:

(4..125) where K p and K d are the proportional and derivative gain, respectively, Hence, steady-state errors due to environmental disturbances and neglected dynamics cannot be compensated for. The main reason for omitting integral action is that the rudder servo has an on-off or relay nonlinearity which would cause a limit-cycle (chattering) if integral action is added. However, the magnitude of the steady-state errors, 1-2 degrees, was in the same order as the accuracy of the flux-gate compass. Consequently, there is not much gained by including integral action. For the PD-controllaw the loop-transfer junction (4·122) reduces to: l(s)

+ K d s)(1 + T 3 s) s(1 + Trs)(1 + T2 s)

__ K(Kp

(4.126)

Hence, tuning ojthe yaw controller (4.125) in terms of the controller parameters K p and K d can easily be done by plotting l(jw) in a Bode-diagram. The performance of the autopilot is shown in Figure 4,13 for a long time mission where the course is changed in steps of 10 degrees each 10 minutes. A typical step response is shown in Figure 4- J4 where the heading angle is changed from .1 9S to 185 degrees.

o 4.4.3

Combined Pitch and Depth Control

In this section we will design a control system for decoupled pitch and depth controL This design will also be based on a linearized model of the vehicle.. Diving Equations of Motion The diving equations of motion should include the heave velocity w, the angular velocity in pitch q, the pitch angle 0, the depth z and the stern plane deflection Os Assume that the forward speed is constant and that the sway and yaw modes can be neglected. Hence, the pitch and heave kinematics can be perturbed according to:

(Ba + t:.B) = cos( 0 + t:.
+ t:.q)

(4.127)

Stability and Control of Underwater Vehicles

120

300

-;;;, 250 ~

:3 ;;;

""

200 150 0

150

100

50

200

time (min)

2

,• .1

,"

-4

o

I,

.I,L! 1"1 '

.1. ,"",

r 'f"')

,

r:: I

50

,,;.

'"

..

I

:1 '! 11'

"

-I

Jl

uN

~,;

200

150

100

I,"

A.

' rt'IT

time (min)

Figure 4.13: Long time mission [or the NDRE-AUV where the heading is changed in steps of 10 degrees each 10 minutes.

200,-------------, 195

ou,',

-a

~.O.5

ID

:s!. 190

ID

~

"2.

~1

•

185

-1 5

-2!.----f::.-----:;:'::;---:-::::;----,:::l

1 800.':---;:5::::0----,1"'0'::0,--1"'5'::0,----,2:::l00 time (s)

o

50

" : w1l

100

150

200

time (s)

4,--,----'----,--,---.---..,-----.---,------,---,-----, ,

,

, "

:'" "

, ,

':'

: ,

, ' ,

" " "

,.'. ..

•..

J.I1iv~ . . .

.

-2 .':---==---:':,---:6.;;0,---:!=--:-O-:,---:-:':-=----:1.':4::-0----.=-----,:-:---::::l 20 40 80 100 120 160 180 200 0 time (s)

Figure 4.14: Full-scale experiment showing one typical step response for the NDREAUV. The heading angle is changed 10 degrees fr'om 195 to 185 degrees.

4.4 Decoupled Control Design In steady-state we have that 00

121

= constant,

qo

= 0 and

rjJo

= 0, Hence:

If::"O = cos f::"rjJ f::"q "" f::"q I

(4.128)

for small f::"rjJ. Similarly, the perturbed heave dynamics is:

(io+f::"z) = - sin(Oo+M)( uo+f::"u)+cos(Oo+f::"O) cos( rjJo+f::"rjJ)(wo+f::"w) (4.129)

Using the trigonometric formulas:

sin( 00 + f::,,0) cos( 00 + f::,,0) -

sin 00 cos f::,,0 + cos 00 sin f::,,0 "" sin 00 + cos 00 f::,,0 cos 00 cos f::,,0 - sin 00 sin f::,,0 "" cos 00 - sin 00 f::,,0

(4,130) (4.131)

together with rjJo = 0 yields:

(io+f::"z) = -(sin Oo+cos 00 f::,,0)( Uo +f::"u) +(cos 00 -sin OoM)(Wo + f::"w ) (4.132)

Applying the steady-state condition: Zo

= - sin 00 Uo + cos 00 Wo

(4.133)

to tills expression together with the assumption that 2nd-order terms in f::" can be neglected, finally yields: If::"z = -sinOo f::"u - cos 00 uof::"O+cosOo f::"w -sinOo Wo f::" 0

I

(4.134)

Tills suggests the following linear model (dropping the f::,,-notation for notational simplicity): U

(>11

0'12

0'13

W

(>21

(>22

(>23

q iJ z

=

(>31

(>32

(>33

0

0

1 0

-sOo cOo

0

U

/31

0'24

W

/32

0 0 0 (>34 0 0 -(sOowo + cOouo) 0

q

+

0

z

/33

Os

(4135)

0 0

where aij is found from the general expression for the 6 DOF linear equations of motion and Pi should be determined for the actual stern plane. Simplified Diving Equations of Motion (Zero Pitch Angle)

A further reduction could, be to assume zero pitch (0 0 = 0) constant forward speed (uo = constant) Hence, the above state-space model reduces to:

Stability and Control of Underwater Vehicles

122

[~ ] = [~~: ~~~ ~~: ~] [;]e + [~~ ] iJ Z

0 1

1 0

0 0 -uo 0

0 0

Z

Os

(4.136)

with obvious definitions of aij and bi . Alternatively, we can write this model in terms of the hydrodynamic derivatives a.s: m mXG -

[

Zw Mw

o o

+[

00 0] 0

Zq I. - M q 0 0

mXG -

-Zw -Mw

1 0 0 1

muo - Zq Mq

mXGUO -

0 -1

-1 0

o BG,W

o

l][f].[T]" (' m:

Here BG z = Za - ZB is used to denote the vertical distance between the center of buoyancy and center of gravity. Pitch-Depth Control Design

The pitch-depth controller can be designed with background in the linear model (4.136), by simply choosing y = z. Hence: x

y -

Ax + bos eT x

(4.138) (4.139)

-,

where: (4.140)

eT=[OOOl]

Applying the Laplace transformation to this model yields the tr ansfer function: ':"(s) =

b1s2 + (b 2a 12 - b1a 22 - b2Uo)s + (b 2u oall - b1a 21 u O - b1a 23

OS

S[S3 - (all

+ a22)s2 + (alla22 -

a23 - a21a12)s

+ (alla23 -

+ b2a 13) a21a13)]

(4.141) For simplicity, we will a.ssume that the heave velocity during diving is small and that Xa = O. This is quite realistic since most small underwater vehicles move slowly in the vertical direction. This a.ssumption implies that the linear model (4.137) reduces to: , (=a-ZB) IV Ill-Mrj

o

-uo

0] o [q] e + [~] 0 o

Z

0

Os

(4.142)

.,

'~

:'

~.

--;11

4.4 Decoupled Control Design

123

Consequently, the transfer functions 8/6s and z/ 6s are obtained as follows:

~(s) = -'- UQ 8(s) 6s s 6s (s)

8 Ko -6(s)=" 2( ,; s s- + 0 Wo s + 0.1 0

(4,143)

where the gain constant is Ko = Mo/(Iy - M q) The natural frequency Wo and relative damping ratio (0 for the pure pitching motion are defined as: ,. ,0

-Mq

= 2

VBGzW(Iy - M q)

(4144)

Consequently, the natural period in pitch is: (4145)

From this expression it is seen that a reduction in the moment of inertia (Iy - Mq) or an increase of the vertical distance between the center of gravity and the center of buoyancy BC., and the vehicle's weight W, will reduce the natural pitch period, The terms Mqq and BC z W8 are often referred to as passive damping and restoring forces, respectively, since modifications of these parameters require the vehicle to be redesigned, Similar effects can be obtained by designing an active feedback control system of PID-type for combined pitch and depth control, For instance, the control law:

allows the designer to modify the damping and the restoring forces through the derivative and proportional action in the controller, that is, adjusting the controller gains C; (i = 1..,5), Notice that feedback from w is omitted since this state is not usually measured, In the implementation of the controller, the depth z can be measured by a pressure meter, the pitch angle 8 can be measured by an inclinometer while the pitch rate q requires a rate gyro or a rate sensor. If heave velocity w is measured in addition, the more complex model (4,.137) can be used in the control design instead. One way to obtain velocity measurement in the vertical plane is by simply combining a pressure meter with an accelerometer to form a velocity state estimator. Kalman filter algorithms are well suited for this purpose. A more expensive solution is using a Doppler log for directly obtaining velocity measurements.

-

Stability and Control of Underwater' Vehicles

124

Example 4.4 (The NDRE-AUV Depth Control System) The NDRE-A UV depth contr'olleT was designed without integml action accoTding to:

85 = C) We also notice that into (4.143) yields.:

()d

(Zd -

z) - C 3 ()

-

C4 q

(4.147)

= 0 in this implementation. Substituting this expTession

(4.148)

Z

() (deg)

(m)

60

5 0

40

-5 -10

20

-15 00

2

4

8 (deg)

6

8 (minutes)

-20 0

2

468 (minutes)

q (degjs)

4

2

2

1

.. .~

-1 -4 0

2

4

6

8 (minutes)

-2 0

2

4

6 8 (minutes)

"

Figure 4.15: Full-scale depth changing maneuver for the NDRE-AUV. The bias in the pitch rate time series is due to a small off-set in the rate sensor.

Hence, we can choose C 3 and C 4 such that the closed-loop pitch dynamics is stable. Next, we can use the Telationship.:

Z(s) = _ uo ()(s) s to tune C). MaTeoveT, the cubic characteristic equation:

(4 . 149)

4.5 Advanced Autopilot Design for RaVs

125

must have all its roots strictly in the left half-plane to ensure that (z = Zd = constant) in steady-state. However, environmental disturbances can cause steadystate erron for this approach since integral action is omitted. A full-scale depth changing maneuverfor the NDRE-AUV is shown in Figure 4·15.

o

4.5 4.5.1

Advanced Autopilot Design for ROVs Sliding Mode Control

Sliding control has been applied successfully in the control of underwater vehicles by Yoerger and Slotine (1984, 1985) who propose to use a series of single-input single-output (SISO) continuous time controllers. Recent work by Yoerger and Slotine (1991) discusses how adaptive sliding control can be applied to underwater vehicles. Cristi et al. (1990) have applied an adaptive sliding mode controller to control an AUV in the dive plane. Sliding mode controllers have been successfully implemented for the JASON vehicle, Woods Hole Oceanographic Institution by Yoerger, Newman and Slotine (1986) and the MUST vehicle at Martin Marietta, Baltimore by Dougherty, Sherman, Woolweaver and Lovell (1988) and Dougherty and Woolweaver (1990). Besides this, successful implementations have been reported for the NPS AUV II at the Naval Postgraduate School, Monterey by Marco and Healey (1992) and Healey and Lienard (1993) All these experiments show that sliding mode controllers have significant advantages to traditional linear control theory. Single-Input Single-Output (SISO) Affine Systems Yoerger and co-authors propose to use a simplified ROV model: MiiVi

+ ni(vi) = Ti

and

iJi

= Vi

(i

= 1..6)

(4.151)

where all kinematic and dynamic cross-coupling terms terms are neglected. Here, ri is the input, M ii is the diagonal element of the inertia matrix M and ni corresponds to the quadratic damping term in the nonlinear vector n, that is:

.M=

["1"

0 m-Y,j 0 0 0 0

0 0 m-Ztit

0 0 0

0 0 0

Ix -

J(p

0 0

0 0 0 0 Iy0

M,

(4.. 152)

"IN"

Uncer tainties in the model are compensated for in the control design.. For notational simplicity, let us write the ROV model according to:

-

._--

-------

Stability and Control of Underwater Vehicles

126

mx + d Ixlx = T where m> 0; d> 0 = Ti, m = NIii and d Ixlx = ni(vi). We also assume

Here x = 7)i, T and x are measured.

(4.154) that both

x

,>0

s=o ,<0

slope ).

Figure 4.16: Graphical interpretation of the sliding surface. Define a scalar measure of tracking: (40155) where x = x - Xd is the tracking error and ,\ > 0 is the control bandwidth. For s = 0 this expression describes a sliding surface with exponential dynamjcs:

!i(t) = exp( -'\(t - to)) !i(ta)

(4..156)

which ensures that the tracking error x(t) converges to zero in finite time when s = 0 (sliding mode). In fact, the error trajectory will reach the time-varying sliding surface in finite time for any initial condition x(t a) and then slide along the surface towards x(t) = 0 exponentially Hence, the control objective is reduced to finding a nonlineaI control law which ensures that:

.. f,

lim set) - 0

1_00

(4.157)

~

A graphical interpretation of the sliding surface is given in Figure 4.16. In the design of the sliding control law, it is convenient to define a virtual reference x, satisfying: (4..158) Hence, the following expression for m

s is obtained: "

4.5 Advanced Autopilot Design for ROVs

m!;

= =

127

x- m x = ( T - d Ixlx ) - mx, -d IXIs + (r - mx, - d Ixlx,)

m

T

(4.159)

Consider the scalar Lyapunov-like function candidate:

V(8, t)

= ~m 8 2

Differentiating V with respect to time (assuming

v = m!; 8= -d Ixl8

Ih-

2

T

-

m=

0) yields

+ 8(T - mx, - d Ixlx,)

Xr

:i:

IImx""+dIT xX=TII

(4 . 161)

S

x

..

(4.160)

m> 0

,

-

Xd - Ai:

xr' = Id - .Ai s = X - xr

:i:,

x.L,.-

m dl:i:1 Kd I K sgn() I

Figure 4.17: 8180 sliding control applied to underwater vehicles Taking the control law to be:

IT = mx, +dlxlx, - K d s - K where

m and d are the estimates of m and sgn(s) = {

~

-1

sgn(s)!

(4.162)

d, respectively, and: if s> 0 if s = 0

(4.163)

otherwise

yields:

v = -( K

+ d Ixl ) 82 + ( fix, + d lXIX, ) s - K Isl d = d- d Conditions on the switching gain K d

Here fi = m- m and by requiring that V :::; O. The particular choice:

(4.164)

are found

128

Stability and Control of Underwater Vehicles

(4165) with

1]

> 0 implies that: (4.166)

This is due to the fact that ( K d + d lxl ) ;, 0 V x. Notice that, 11 So 0 implies that V(t) So V(O), and therefore that s is bounded. This in turn implies that V is bounded. Hence 11 must be uniformly continuous, Finally, application of Barbiilat's lemma then shows that s -;. 0 and thus x -;. 0 as t -;. 00 Chattering It is well known that the switching term K sgn(s) can lead to chattering. Chattering must be eliminated for the controller to perform properly Slotine and Li (1991) suggest smoothing out the control law discontinuity inside a boundary layer by replacing the sgn(.) function in the control law with: sat(sjcp)

if Is/if;1 > 1 otherwise

= { sgn(s) s/if;

(4 . 167)

where rP should be interpreted as the boundary layer thickness. This substitution will in fact assign a low-pass filter structure to the dynamics of the sliding surface s inside the boundary layer (see below). Moreover, replacing the Ksgn(s) term in (4.162) with If. sat(sjcp) yields the following expressions for the s-dynamics and V. • Inside the boundary layer: m

K

s+ ( K d + d Ixl + ~ ) s =

-

mX r + d lxix,

v ~ -( Kd + d Ixl + ~ )

8

2

(4.168)

(4.169)

• Outside the boundary layer: m s+ (Kd

+ d Ixl ) s + K sgn(s) = mX r + dlxix,

(4170)

(4.171)

The boundary layer thickness can also be made time-varying to exploit the ma..'{imum control bandwidth available. See Slotine and Li (1991) for a closer description on time-varying boundary layers.

4.5 Advanced Autopilot Design for RaVs

129

Example 4.5 (Sliding Mode Control Applied to RaVs) Consider the simplified model of an underwater vehicle in surge: m

with m as"

= 200 kg

and d

x + d Ixlx =

= 50 kg/m.

The SISO sliding controller can be complLied

Ir = mx, +d lxix, - K where K d

~

(4..172)

T

d

s - Ksat(slifJ)

I

(4.173)

O. The following two cases were studied

(1) PD-Controller:

m=O d=O

Kd = 500 K=O

Notice that this simply corresponds to the PD control law: T

= -Kd S = -Kd j; -

.\ K d

x

(4.174)

(2) Sliding Mode Controller:

m =0 . 6m

m

d

d Kd

K

= 1.5 d = (m x, +

I

d Ixl x,)

1+0.1

~

0.5m 0.5 d = 200 ~

In the simulation study the closed-loop bandwidth was chosen as .\ = 1 for both controllers. The boundary layer thickness was chosen as