Ship Kinetics

4 Ship Kinetics by T. Perez and T.I. Fossen Chapter 3 reviewed the kinematics of ship motion, i.e. the geometrical aspects of motion: variables, reference frames and transformations. This chapter reviews the kinetics of ship motion, i.e. the study of the forces acting on the ship and the motion they produce. By combining the material presented Chapter 3 with that of this chapter, one can obtain ship dynamic models for control design and testing (cf. Section 1.3). 4.1 An Overview of Ship Modeling for Control As discussed in [174], the models for ship motion control system design use superposition of either motion or force. Therefore in each case, magnitudes can be conceptually decomposed as X = Xw + X̄, • X̄ = Xs−vd + Xc . (4.1) First-order wave-induced force/motion (w). This force/motion is oscillatory. This is commonly modelled as a time series disturbance obtained by combining the wave spectrum with the vessel Force Response Amplitude Operators (Force RAO) or the Motion Response Amplitude Operators (Motion RAO), which are transfer functions that map the wave elevation or wave slope into force and motion. • Slowly-varying disturbance force/motion (s-vd). This force/motion is produced by second-order wave effects (wave mean-drift and slowlyvarying forces), current and wind. • Control-induced force/motion (c). This is the force/motion induced by the control system, which is usually designed to counteract only the the effect of the slowly-varying disturbances. 60 4 Ship Kinetics The motion superposition model is the most commonly adopted model for control system design—see Figure 4.1. In this model, a manoeuvring model is used to describe the relationship between the control action and the motion induced by this control action whereas seakeeping models are used to describe the motion due to the waves. Se a k e e ping m ode l (m ot ion in w a ve s) Wave Spectrum Motion RAO Motion spectrum Σ t M a noe uvring m ode l (c a lm w a t e r m ode l) Control forces and moments Linear equations of motion based on hydrodynamics derivatives (ω → 0) Low-frequency motion First-order linear wave motion time series (wave-frequency motion) Reference frame transformation Total motion Slowly-varying disturbance motion Nonlinear terms (Damping, Coriollis) Fig. 4.1. Motion superposition model of a marine vessel [174]. The type of model shown in Figure 4.1 presents two shortcomings. The first one is that the model may not be used for multibody system interactions. This requires energy exchange, i.e., it requires elements with common forces and speeds, and this is not captured by a models that uses wave-induced motion as a disturbance. This issue, however, is more relevant to marine operations like up-loading and deploying heavy equipment, pipe-lying, positioning of deep suspended loads, than to the problems considered in this book. The second shortcoming of this type of model is that the manoeuvring part does not incorporate fluid memory effects associated with the wave-frequency induced motion. Indeed, the radiation forces due to the frequency dependent mass and damping of the ship are only considered in the seakeeping model (this is embedded in the motion RAO). This results in miss-modelled dynamics, which are of interest for ship motion control in a seaway. An alternative approach consists of using a model with force superposition rather than motion superposition, as illustrated in Figure 4.2. In this model, the Force RAO are combined with the sea spectrum to give the wave-excitation forces. Also, a time-domain representation is used for the fluid memory effects 4.1 An Overview of Ship Modeling for Control 61 associated with the wave-frequency motion of the ship. This time domain representation uses a convolution integral to compute the radiation damping forces. Therefore, these forces are computed using both the wave- and controlinduced motion instead of only the wave-induced motion as in the model shown in Figure 4.1. This is key issue for obtaining good models for wavemotion damping control problems like roll stabilisation. The model shown in Figure 4.2 is well known in marine technology and it is part of the state of art time-domain ship motion simulators. However, its use for control system design has not yet been widely adopted. Recent reported results [127] proposed a state-space representation for the convolution integral which results in model amenable to control system design tools. Forc e supe rposit ion m ode l (U nifie d se a k e e ping-m a noe uve ring m ode l) Wave spectrum Force RAO First-order wave exciting loads Control Forces and Moments Wave exc. spectrum Σ t Reference frame transformation Linear equations of motion with memory effects Motion (ω ≥ 0) Slowly-varying exciting loads Nonlinear terms (damping, coriollis) Fig. 4.2. Force superposition model of a marine vessel [174]. In the rest of the chapter, we will present and discuss the elements shown in Figures 4.1 and 4.2. We will first describe the seakeeping and manoeuvring models so as to obtain the classical motion superposition model shown in Figure 4.1. Then, we will discuss the elements of the force superposition model shown in Figure 4.2. 62 4 Ship Kinetics 4.2 Seakeeping Theory Models As already mentioned in Section 3.3.3, seakeeping theory studies the motion of surface vessels in waves. The motion of a ship in a seaway is the response to sea loads, which are forces and moments that arise due to changes in pressure over the surface of the hull. If the sea state is not extreme, the motion responses can be considered within a linear framework. This allows one to apply the principle of superposition to study ship motion: the motion due to an irregular sea can be described as the summation of the responses to many regular waves. This concept was introduced to the field of ship motion study by St. Denis and Pierson in the 1950s [213], as a method to predict ship responses in a realistic seaway. The linear seakeeping theory of ship motion makes three essential assumptions: • The sea surface elevation is assumed to be a realisation of an ergodic Gaussian stochastic process with zero mean. Thus, the process is entirely described by its power spectral density—the sea spectrum. • The wave-excitation loads and ship motion response are assumed to be linear. • The ship keeps a steady course and moves at a constant average speed (which includes the case of zero speed). The first assumption was already discussed in Chapter 2. The second assumption allows one to obtain the power spectral density of the waveexcitation loads and the motion components by multiplying the response amplitude operators (force and motion RAO) by the wave spectrum. Since the excitation is Gaussian and the motion response linear, the motion is also Gaussian; then, once the motion spectrum is obtained, all the statistics of motion can be calculated and data for time series generation be obtained. The linear assumption is a limitation since it neglects viscous effects (which need to be incorporated separately) and characteristics of a real free surface. Despite this, linear theory has proven to be a tool that yields reasonable predictions for analysis at preliminary stages of both ship and ship motion control system design. For roll motion, in particular, viscous effects play a very important role because of bilge keels and appendages, and these effects cannot be neglected for they provide most of the total roll damping. This leads to nonlinearities, which can be reduced to linear equivalent terms via stochastic or equivalent linearisation [135, 182, 192]; and hence, the problem can be still considered within a linear framework. 4.2 Seakeeping Theory Models 63 4.2.1 Equations of Motion and Hydrodynamic Forces in the h-frame A distinctive characteristic of seakeeping theory is the use the hydrodynamic frame, h-frame (xh , yh , zh ), to describe the motion of the fluid and the ship— see Section 3.1. Since, the h-frame is inertial, the vector equation of motion in this frame is MhRB ξ̈ = τ hhyd . (4.2) The coordinates ξ are the generalised perturbation coordinates or seakeeping coordinates given in the h-frame defined in Section 3.3.3. The components of the vector τ hhyd are the generalised hydrodynamic forces (forces and moments) for the six degrees of freedom expressed in the h-frame. The matrix MhRB is the generalised rigid-body mass matrix (mass and inertia) with respect to the h-frame:   mI3×3 −mS(rhg ) MhRB  Ih mS(rhg ) ⎤ ⎡ m 0 0 0 mz hg −my hg ⎢ 0 0 mxhg ⎥ m 0 −mz hg ⎥ ⎢ (4.3) h h ⎢ 0 0 m my g −mxg 0 ⎥ ⎥ ⎢ . =⎢ h h Ixh −I hxy −I hxz ⎥ ⎥ ⎢ 0 −mz g my g h h h h h ⎣ mz 0 −mx −I I −I ⎦ g −my hg mxg g 0 yx y −I hzx −I hzy yz Izh The mass of the ship m can be calculated as the product of the water density and the displaced volume ∇, i.e. m = ρ∇ (where ρ is the water density). The coordinates of the centre of gravity with respect to the h-frame are given by the vector rhg = [xhg , ygh , zgh ]t , and S(·) is the skew-symmetric matrix defined in (3.3). The inertia tensor with respect to the h-frame Ih is ⎤ h Ixx −I hxy −I hxz h Ih  ⎣−I hyx Iyy −I hyz ⎦ , h h h −I zx −I zy Izy ⎡ Ih = (Ih )t > 0, (4.4) where the moments and products of inertia can be calculated using the displaced volume as follows:   Ixx = [(y h )2 + (xh )2 ]ρdV, Ixy = Iyx = (y h xh )ρdV, ∇ ∇   h 2 2 Iyy = [(z ) + (x ]ρdV, Ixz = Izx = (xh z h )ρdV, (4.5) ∇ ∇   [(y h )2 + (z h )2 ]ρdV, Iyz = Izy = (y h z h )ρdV. Izz = ∇ ∇ 64 4 Ship Kinetics The total hydrodynamic force vector τ H is assumed to be the linear combinations of the following components [63]: τ hyd = τ 1w + τ 2w + τ r + τ v + τ hs . (4.6) • First-order wave excitation forces (τ 1w ). These forces are the zeromean oscillatory forces caused by the waves, which are separated into two components. The first component gives the, so-called, Froude-Kriloff forces, which are forces due to the incident waves under the assumption that the hull is restrained from moving and that the presence of the hull does not disturb the flow field. The second effect is a correction to account for the modification of the flow field due to the hull; otherwise there would be water mass transfer through the hull. This second component gives the so-called diffraction forces. • Second-order wave excitation forces (τ 2w ). These forces include mean wave-drift loads, slowly varying (difference frequencies) and rapidly varying (sum frequencies) wave loads. Depending on the ship motion control problem considered these forces represent a significant contribution of the wave excitation loads. • Radiation forces (τ r ). These forces appear as a consequence of the change in the momentum of the fluid and the waves generated due to the motion of the hull. These forces are proportional to the accelerations of and velocities the ship. Due to this, the radiation forces are separated into the so-called added-mass forces (forces proportional to accelerations) and potential-damping forces (forces proportional to velocities). • Viscous forces (τ v ). These are nonlinear damping forces that appear due to nonlinear non-conservative phenomena by which kinetic energy of the hull is transferred to the fluid due to viscous effects (skin friction, flow separation and eddy making). These forces are depend nonlinearly on the relative velocities between the hull and the fluid. • Hydrostatic forces (τ hs ). These are restoring forces due to gravity and buoyancy that tend to preserve the equilibrium of the ship. Except for the viscous forces, the rest can be studied within a linear framework. Indeed, the wave excitation forces can be obtained by considering linear potential flow theory. As already discussed in Section 2.1.4, this theory seeks a potential function Φ(x, y, z, t), that satisfies the Laplace equation in the fluid and boundary conditions on the surface of the hull, on the free surface of the water and at the seafloor. Due to linearity, this potential is separated into a radiation potential and a wave-excitation potential: 4.2 Seakeeping Theory Models Φ = Φr + Φw . 65 (4.7) The pressure can be calculated by substituting the potential into the linearised Bernoulli Equation (see Section 2.1.4), and the hydrodynamic forces and moments are obtained by integrating the pressure over the wetted surface Sw of the hull. Thus, the components, for each degree of freedom, of the radiation forces can obtained as follows:   ! r − Sw ∂Φ + gz (n)i ds i = 1, 2, 3. h ∂t !  τri = (4.8) ∂Φr − Sw ∂t + gz (r × n)i ds i = 4, 5, 6. Where n = [n1 , n2 , n3 ]t is a unit vector normal to the wetted surface of the hull and positive into the fluid. The notation (n)i and (r × n)i refers to the i-th component of n and r × n respectively. The first-order wave excitation forces and moments are obtained in a similar way using the wave-excitation potential:   ! − Sw ∂Φ∂t1w + gz (n)i ds i = 1, 2, 3. h !  (4.9) τ1wi = − Sw ∂Φ∂t1w + gz (r × n)i ds i = 4, 5, 6. The first-order wave excitation potential is further separated into the incident potential and the diffraction potential Φ1w = ΦI + ΦD . The incident potential corresponds to the undisturbed incident waves, cf. (2.26). The integration of the pressure due to this potential gives the so-called Froude-Kriloff forces, whereas the integration of the pressure due to diffraction potential gives the wave diffraction forces—see [195, 63] for details. Within the first order potential theory, these forces depend linearly on the magnitude of the waves. As we shall see in the next section, expressions (4.8) and (4.9) can be evaluated numerically for sinusoidal waves using hydrodynamic programs. Because of the linearity assumption made, these forces are proportional to the amplitude of the sinusoidal waves, and the data obtained is tabulated as transfer-function model (force RAO)—see Section 4.2.2. The second-order wave excitation forces are obtained by integration the pressure due to second order potentials, but this goes beyond the scope of this book. These forces, however, have a little impact for the ship motion control problems considered in this book. The reader is referred to [195, 63] for further details of second-order wave excitation forces. Finally, the hydrostatic forces and moments are proportional to the displacements ξ. Thus, the linearised forces and moments can be expressed as τ hhs = Gh ξ. The only non-zero linear restoring coefficients are [63] (4.10) 66 4 Ship Kinetics Gh33 = ρgAW P h C53 Gh35 = Gh44 Gh55 = ρg∇GM t = −ρg  xh ds AW P (4.11) = ρg∇GM l, where AW P is the water-plane area, ∇ is the displaced volume, and GM t and GM l are the transverse and longitudinal metacentric heights. Figure 4.3 illush trates the components involved in the roll restoring moment (τ4hs = Gh44 ξ4 ). In this case, the restoring moment is the product of the righting arm GZ ≈ GM t ξ4 times the force due to the displaced volume ρg∇. GM t M t-Trans. Metacentre ξ4 = φ CG Water level ρg∇ GZ CB Fig. 4.3. Transverse metacentric height GM t and transverse righting arm GZ. 4.2.2 Wave Force Response Amplitude Operator (Force RAO) As already mentioned, the wave-excitation forces can be studied using superposition due to the linearity assumption made in seakeeping theory. Therefore, it is common to study the forces and the motion only for sinusoidal excitations and use the frequency domain approach to calculate force transfer functions, which can then be combined with the wave spectrum to obtain the loads for a particular sea state and sailing condition. If the sea elevation at the origin of the h-frame is given by ζ(t) = ζ̄ cos(ωe t + ǫ), (4.12) using complex notation we can express the harmonic wave, the wave excitah tions and the components of motion in terms of the complex variables ζ̃, τ̃1wi ˜ and ξi , such that 4.2 Seakeeping Theory Models # "   ζ(t) = ℜ ζ̃ejωe t = ℜ ζ̄ejǫ ejωe t " #  h jωe t  h h h τ1wi = ℜ |τ̃1wi (t) = ℜ τ̃1wi e | ej(arg τ̃1wi +ǫ) ejωe t , # # " " ξi (t) = ℜ ξ˜i ejωe t = ℜ |ξ˜i | ej(arg ξ̃i +ǫ) ejωe t , 67 (4.13) where ℜ{·} denotes the real part of the argument. Using these, we can define the wave-force response amplitude operators (Force RAO) for each degree of freedom as   h  τ̃ (ωe , χ)  j arg τ̃ h (ω ,χ) e e 1wi Fi (ωe , χ) =  1wi ,  ζ̃ or, expressed in the wave frequency domain,   h  τ̃ (ω, U, χ)  j arg τ̃ h (ω,U,χ) e 1wi . Fi (ω, U, χ) =  1wi  ζ̃ (4.14) (4.15) h Note that the complex wave excitation forces τ̃1wi depend on the encounter frequency ωe and also the encounter angle χ. These force transfer functions can be obtained from standard seakeeping programs, which calculate them based on the geometry of the hull and loading condition. Details of this are beyond the scope of this book; here we will only use of the results of these programs, but the interested reader can see for example [64, 115, 144] and references therein. As an example, Figure 4.4 shows the magnitude of the force RAO in sway for different headings corresponding to the benchmark example given in Appendix B at 15kt. With these force RAO, one can obtain the power spectral density of the wave loads given the sea spectrum Sζζ (ωe , χ): Sτ1W τ1W i (ωe ) = |Fi (ωe , χ)|2 Sζζ (ωe , χ), or simulate first-order wave excitation loads due to irregular seas as  h τ1wi (t) = |Fi (ωej , χ)| ζj cos(ωej + ǫ + arg[Fi (ωej , χ)]), (4.16) (4.17) j for i = 1, 2, . . . , 6. 4.2.3 Motion Response Amplitude Operator (Motion RAO) Apart from evaluating the force RAO, as discussed in the previous section, one can also use the frequency domain approach to evaluate the radiation forces, and then the equation of motion (4.2). This leads to the so-called motion response amplitude operators (motion RAO) 68 4 Ship Kinetics FRAO2 (Sway force amp/wave amp) [N/m] for U = 15kt x 10 5 9 8 7 6 5 4 3 2 1 0 3 2.5 150 2 1.5 100 1 50 0.5 Frequencies [rad/s] 0 0 Enconuter angle [deg] Fig. 4.4. Sway force RAO of the naval vessel benchmark at 15kts for different headings. The data is given as a function of the wave frequency. These force RAO were obtained using by ShipX-VERES [64]. The radiation forces are proportional to the accelerations and velocities of the ship motion. For sinusoidal motion, the vector of radiation forces (with components given in Expression (4.8)) can be expressed as follows [195, 63]: τ hr = −Ah (ωe )ξ̈ − Bh (ωe )ξ̇, (4.18) where—if the ship is symmetric with respect to the xh -zh plane—the matrix Ah (ωe ) is of the form [195]: ⎤ ⎡ h A11 (ωe ) 0 Ah13 (ωe ) 0 Ah15 (ωe ) 0 ⎢ 0 Ah22 (ωe ) 0 Ah24 (ωe ) 0 Ah26 (ωe )⎥ ⎥ ⎢ h h h ⎢ A31 (ωe ) 0 A33 (ωe ) 0 A35 (ωe ) 0 ⎥ ⎥ , (4.19) Ah (ωe ) = ⎢ ⎢ 0 Ah42 (ωe ) 0 Ah44 (ωe ) 0 Ah46 (ωe )⎥ ⎥ ⎢ h ⎣A51 (ωe ) 0 Ah53 (ωe ) 0 Ah55 (ωe ) 0 ⎦ 0 Ah62 (ωe ) 0 Ah64 (ωe ) 0 Ah66 (ωe ) 4.2 Seakeeping Theory Models 69 and the damping matrix Bh (ωe ) is similar to (4.19) modulo substitution Ahik h by Bik . Due to the relationship (4.18), the matrix A(ωe )h is called added-mass matrix, and the first term in (4.18) gives the so-called added mass forces. The added mass forces reflect the change of momentum in the fluid due to the motion of the hull—it does not mean that the water moves along with the ship. The matrix B(ωe )h is called potential-damping matrix . Since potential theory does not account for viscous effects (skin friction and flow separation effects), the second term in (4.18) represents the transfer of kinetic energy of the hull into the generated waves that appear due to the motion of the hull. In the case of roll and surge motion, the potential damping is very small and corrections need to be made to account for viscous effects. These corrections are based on empirical approaches, and can be incorporated into the equations of motion via an equivalent linear viscous damping coefficient h h and B44visc such that the energy dissipated through viscous effects is B11visc the same as that dissipated by the linear viscous term [135, 182]. By substituting the wave-excitation, radiation and hydrostatic forces into (4.2), and using the complex notation (4.13), we obtain the frequency-domain equations of motion: −ωe2 [MhRB + Ah (ωe )]ξ̃ + jωe Bh (ωe )ξ̃ + Gh ξ̃ = τ̃ h1w (ωe , χ), (4.20) from which the the responses can be evaluated: ξ̃ = [−ωe2 (MhRB + Ah (ωe )) + jωe Bh (ωe ) + Gh ]−1 τ̃ h1w (ωe , χ). (4.21) It should be emphasised that by multiplying both sides of Expression (4.20) by ejωe t and taking the real part, one would not obtain the true equations of motion. Indeed, this would result in time-domain equations which describe the motion of the ship only if the wave excitation forces are sinusoidal, and provided the coefficients assume a proper value according to the frequency of the excitation. The true time-domain equations of motion (valid for any type excitation) have constant added mass and the damping is given by a convolution integral. This will be further discussed in Section 4.4. In practice, Equation (4.21) is evaluated numerically only for a discrete set of frequencies (usually between 20 and 40). Details on how this is calculated is beyond the scope of the book, and interested reader should consult, for example, the hydrodynamic programs described in [64, 115, 144]. Due to the linearity of (4.21) and also the linearity of the wave excitation forces, the results give the amplitude and phase of each component of motion per unit of wave height or wave slope. These give the frequency response of the ship, which is represented in terms of the so-called Motion Response Amplitude Operators motion RAO: 70 4 Ship Kinetics    ξ˜ (ω )   i e  j arg ξ̃i (ωe ) Hi (ωe , χ) =  e  ζ̃  (for i = 1, 2, . . . , 6) (4.22) Note that it is also common practice to give the motion RAO of rotational motions (roll, pitch and yaw) normalised in terms of the wave slope (ke ζ̄). In this case, the motion RAO is defined as    ξ˜ (ω )   i e  j arg ξ̃i (ωe ) Hi (ωe , χ) =  (for i = 4, 5, 6), (4.23) e  ke ζ̃  where ke = ωe2 /g is the regular wave number—see Section 2.2. The motion RAO can be expressed either in the wave frequency domain or in the encounter frequency domain: Hi (ωe , χ) ≡ Hi (ω, U, χ). Then, for example, if the sea elevation, at the origin of the h-frame is described by ζ(t) = ζ̄ cos(ωe t + ǫ), (4.24) the motion components are, then, obtained as ξi (t) = ζ̄|Hi (ωe , χ)| cos(ωe t + arg[Hi (ωe , χ)] + ǫ) for i = 1, 2, . . . , 6 (4.25) If the motion RAO of the rotational components are normalised by the wave slope, (4.25) can be replaced by ξi (t) = ζ̄|Hi (ωe , χ)|ke sin(ωe t + arg[Hi (ωe , χ)] + ǫ) for i = 4, 5, 6. (4.26) In the formulation presented thus far, the motion RAO obtained are given at the origin of the h-frame. To obtain the motion RAO at any other location point of interest x on the ship given by the vector rhx = [xhx , yxh , zxh ], relative to the h-frame, the following relationship can be used: [H1x , H2x , H3x ]t = [H1 , H2 , H3 ]t + [H4 , H5 , H6 ]t × rhx , (4.27) where H1x , H2x , and H3x are the motion RAO at the desired location. Figure 4.5 shows the roll, sway and yaw motion RAO, at the center of gravity, for one speed and different encounter angles of the benchmark example vessel described in Appendix B. When the phases are negative, it means that the motion lags the wave elevation (or slope) at the origin of the h-frame. As described in Chapter 2, 4.2 Seakeeping Theory Models 71 when sailing in beam seas, the encounter and wave frequency are the same, and we can see from Figure 4.5 that the roll natural frequency is close to 1 rad/s. When sailing in bow seas, wave frequencies are mapped in to higher frequencies. This is why in Figure 4.5 the wave frequencies that are close to the resonance in roll are lower than those in the case of beam seas. For the sailing conditions shown, no wave frequencies are close to the roll resonance in quartering seas. Figures 4.6, 4.7 and 4.8 show the motion RAO for encounter angles with increments of 15 deg. It should be noted from these examples that for the sway and yaw motion RAO in quartering seas with encounter angles of 30 deg and 45 deg, the encounter frequency approaches zero for the wave frequencies approaching 1.5 rad/s and 2.5 rad/s respectively, and then the response in sway and yaw increases rapidly. This can be seen from (4.21), by noting that there are singularities when ωe → 0 because there are no restoring coefficients G2k or G6k , and. For simulations in irregular seas, this may not be a problem if the wave spectrum has low energy close to those frequencies. In practice, the helmsman or the autopilot will always limit these large motions [135]. 4.2.4 Ship Motion Spectra and Statistics of Ship Motion The motion RAO can be combined with the sea elevation or sea slope spectrum to obtain the motion spectra: Sξξi (ωe ) = |Hi (ωe , χ)|2 Sζζ (ωe ) for i = 1, 2, . . . , 6. (4.28) If the motion RAO of the rotational components are normalised by the wave slope, a similar expression holds for the wave slope spectrum (cf. (2.33)): Sξξi (ωe ) = |Hi (ωe , χ)|2 S′ζζ (ωe ) for i = 4, 5, 6. (4.29) Since the power of any magnitude is invariant with respect to the reference frame from which it is observed, it follows that the spectral moments of order n of the motion components can be computed either in the encounter or in the wave frequency domain:  ∞ n ωen |Hi (ωe , χ)|2 Sζζ (ωe ) dωe mξi = (4.30) 0 ∞ ω n |Hi (ω, U, χ)|2 Sζζ (ω) dω. = 0 It is implicitly assumed in the above that either the wave spectrum or the wave slope spectrum should be used according to the definition of the MRAO of each motion component, and we will not make any explicit difference henceforth. The factor |Hi (ω, U, χ)|2 Sζζ (ω) in the second line of (4.30) is called a pseudo-spectrum and will be denoted by 72 4 Ship Kinetics Roll Sway Yaw 6 60deg 90deg 135deg |H2 (ω, U, χ)| |H4 (ω, U, χ)| 5 4 3 2 0 1 2 4 0.6 3 0.4 2 0 3 0.2 0 1 ω 0 −100 2 3 0 −200 2 3 2 3 100 −100 ω 1 200 100 1 0 ω arg H6 (ω, U, χ) 100 0 0 3 200 arg H2 (ω, U, χ) arg H4 (ω, U, χ) 2 ω 200 −200 60deg 90deg 135deg 0.8 1 1 0 60deg 90deg 135deg 5 |H6 (ω, U, χ)| 6 1 0 −100 0 1 2 3 −200 ω 0 1 ω Fig. 4.5. Roll, sway and yaw motion RAO of the naval vessel benchmark at 15kts for different encounter angles. The motion RAO are given as a function of the wave frequency. S⋆ηηi (ω) = |Hi (ω, U, χ)|2 Sζζ (ω). (4.31) The pseudo-spectrum can be used to determine the moments of the motion component as readily seen in (4.30), yet it is not the actual psd of the motion component observed from the ship. To obtain the latter, the pseudo-spectrum must be converted to the encounter frequency, using (2.62). Figure 4.9 depicts a particular example for roll motion in quartering seas. The top and middle plots on the left hand side show the ITTC spectrum and the corresponding wave slope spectrum for a particular significant wave height and average wave period. The bottom left plot shows the square of the magnitude of the roll motion RAO for the particular sailing condition. The top plot on the right-hand side shows the roll pseudo-spectrum according to (4.31). The middle plot on the right hand side shows the transformation between wave and encounter frequency for the adopted sailing condition according to |H4 (ω, U, χ)| 4.2 Seakeeping Theory Models 73 30deg 45deg 60deg 75deg 90deg 105deg 120deg 135deg 150deg 165deg 8 6 4 2 3 0 200 2.5 2 150 1.5 100 χ 1 50 ω 0.5 0 0 Fig. 4.6. motion RAO Roll of the naval vessel benchmark at 15 kt for different encounter angles. The motion RAO are given as a function of the wave frequency. (2.29), and the bottom plot shows the roll power spectral density, which is the transformation of the pseudo-spectrum according to (2.62). Because of the frequency content of the pseudo-spectrum and the transformation to the encounter frequency for the particular sailing condition, most of the frequency components falls close to ωe =4.5 rad/s—see Figure 2.3 for details on how to calculate this value. Due to the singularity that can appear in the case of quartering seas, as shown in Figure 4.9, the encounter spectrum is seldom used in computations; all statistics are calculated using the pseudo-spectrum. We will see in the next section that the encounter spectrum is not necessary to simulate the time series or motion. Figure 4.10 shows a similar example but for the bow sea case. As we can see from these two examples, the roll power spectral density can vary significantly depending on sea state, speed and encounter angle. This can result in problems for rudder-based roll stabilisers as we shall see in the second part of the book. 4.2.5 Time-series of Ship Motion using Seakeeping Models The ship motion in a seaway can be simulated by time series. The method is the same as that presented in Chapter 2, Section 2.10, to simulate sea surface elevation. 74 4 Ship Kinetics 30 45 60 75 90 105 120 135 150 165 |H2 (ω, U, χ)| 2 1.5 1 0.5 3 0 200 2.5 2 150 1.5 100 1 χ 50 ω 0.5 0 0 Fig. 4.7. motion RAO sway of the naval vessel benchmark at 15 kt for different encounter angles. The motion RAO are given as a function of the wave frequency. 30 45 60 75 90 105 120 135 150 165 |H6 (ω, U, χ)| 1.5 1 0.5 3 0 200 2.5 2 150 1.5 100 χ 1 50 ω 0.5 0 0 Fig. 4.8. motion RAO yaw of the naval vessel benchmark at 15 kt for different encounter angles. The motion RAO are given as a function of the wave frequency. 4.2 Seakeeping Theory Models 6 S⋆ξξ4 (ω) Sζζ (ω) [m2 s] 1 0.5 x 10 75 −3 RMS roll:2.8103 [deg] 4 2 ITTC Hs [m]:2.5, T [s]:7.5 0.8 1 1.2 ω[rad/s] 0 1.4 0 1 0 1 ω[rad/s] 2 3 2 3 1.5 ωe [rad/s] S′ζζ (ω) 0 0.4 0.6 −3 x 10 3 2 1 1 0.5 ITTC Hs [m]:2.5, T [s]:7.5 0 0.4 0.6 0.8 1 1.2 ω[rad/s] 0 1.4 0.1 Sξξ4 (ωe ) |H4 (ω)|2 2 ω[rad/s] Speed [kt]: 15, Enc. Angle [deg]:45 1.5 1 RMS roll:2.515[deg] 0.05 0.5 0 0 1 ω[rad/s] 2 3 0 0 0.5 1 We [rad/sec] 1.5 Fig. 4.9. Roll motion encounter psd in quartering seas. By formulating the problem in the wave encounter frequency domain, the time series for the different components can be generated as indicated in the following: ξi (t) = M N   n=1 m=1 η̄inm (ωn ) cos    ω2 U ωn − n cos(χm ) t + ϑinm (ω) + ǫn , g (4.32) for i = 1, 2, . . . 6, with  η̄inm (ωn∗ ) = 2|Hi (ωn∗ , U, χm )|2 Sζζ (ωn∗ , χ∗m )ΔχΔω. ϑinm (ω) = arg Hi (ωn∗ , U, χ∗m ), and ωn∗ chosen randomly in the interval   Δω Δω , ωn + ωn − . 2 2 (4.33) (4.34) 4 Ship Kinetics 0.08 S⋆ξξ4 (ω) 1 0.5 S′ζζ (ω) 0 0.4 0.6 −3 x 10 3 0.8 1 1.2 0 1 ω[rad/s] 0.6 0.8 1 1.2 2 3 4 2 0 1.4 0 1 ω[rad/s] 0.04 Sξξ4 (ωe ) 40 |H4 (ω)|2 3 6 ω[rad/s] 30 Speed [kt]: 15, Enc. Angle [deg]:135 20 10 0 2 ω[rad/s] 8 ITTC Hs [m]:2.5, T [sec]:7.5 0 0.4 RMS roll:5.4802 0.04 0 1.4 2 1 0.06 0.02 ITTC Hs [m]:2.5, T [s]:7.5 ωe [rad/s] Sζζ (ω)[m2 s] 76 0.03 RMS roll:5.4802 0.02 0.01 0 1 2 ω[rad/s] 3 0 0 1 2 3 We [rad/sec] Fig. 4.10. Roll motion encounter psd in bow seas. The phases ǫn are independent and identically distributed with uniform pdf in the interval [0, 2π], and once chosen, these are the same for all the different motion components. If the motion RAO of the rotational components are normalised by the wave slope, one can use the following expression for i = 4, 5, 6: M N      ωn2 U cos(χm ) t + ϑinm (ω) + ǫn . η̄inm (ωn )kn sin ωn − ξi (t) = g n=1 m=1 (4.35) The expressions given above are valid for the general case of short-crested seas, provided the directions are discretised and the RAO are calculated for each of the discrete directions χm . For the case of a long-crested sea M = 1, and Δχ = 1. To chose the number of components N and Δω, the rules given in Section 2.10 can be used. As an example of time series, using the MRAO consider the results shown in Figures 4.11 and 4.12 for long-crested seas. These figures show the time series for roll, sway and yaw for beam seas, and roll time series for quartering, beam and bow seas. For the simulation of the wave excitation force times 4.2 Seakeeping Theory Models 77 series, a similar procedure can be followed by combining the sea spectrum with FRAO. ITTC Hs [m]:2.5, T [s]:7.5, Speed [kt]: 15, Enc. Angle [deg]:90 20 ξ4 (t)[deg] Sξξ4 (ωe ) 0.08 0.06 RMS roll:8.3597 0.04 0.02 0 0 1 ωe [rad/s] 2 3 10 0 −10 −20 ξ2 (t)[m] 1 RMS Sway:0.654 0.5 0 Sξξ6 (ωe ) 4 0 −4 x 10 1 40 0 20 40 0 20 40 t[s] 60 80 100 ωe [rad/s] 2 60 80 100 60 80 100 0 −1 −2 3 RMS yaw:0.46027 0 1 t[s] 2 2 0 20 2 ξ6 (t)[deg] Sξξ2 (ωe ) 1.5 0 1 ωe [rad/s] 2 3 1 0 −1 −2 t[s] Fig. 4.11. Roll sway and yaw motion power spectral densities and time series for beam seas at 15kts. The wave spectrum used is ITTC with H1/3 = 2.5 m and T = 7.5 s. It should be noted that the time series can also be implemented using linear shaping filters driven by white noise as described in Section 2.6. Depending on the vessel and sailing conditions the, psd of the motion components can be approximated by simple second-order shaping filters. For these cases, the algorithm given in Section 2.6 can be used to tune such filters according to the given sea state and sailing conditions. For cases in which the second-order approximation does not yield good results, higher order models should be estimated—see, for example, [37]. With this, we have defined all the elements of the seakeeping model indicated in Figure 4.1. The following summarises the properties of the seakeeping models [174]: • The motion is described from an equilibrium frame traveling with the average speed of the ship—that is to say, fixed speed and heading. 78 4 Ship Kinetics χ = 135 deg χ = 90 deg 30 20 20 20 10 0 −10 10 0 −10 −20 −30 ξ4 (t)[deg] 30 ξ4 (t)[deg] ξ4 (t)[deg] χ = 45 deg 30 50 t[s] 100 −30 0 −10 −20 0 10 −20 0 50 t[s] 100 −30 0 50 t[s] 100 Fig. 4.12. Roll motion time series under different encounter angles for the benchmark example at 15 kt. The wave spectrum used is ITTC with H1/3 = 2.5 m and T = 7.5 s. • Linearity is assumed between the motion responses and the wave amplitude; the problem is analysed in the frequency domain, i.e. the frequency-domain equations of motion are analysed in steady state and for sinusoidal waves excitations. • To calculate the hydrodynamic radiation-induced forces, potential flow theory is used, and corrections for viscous effects are made based on empirical procedures. • The mass and damping coefficients of the frequency-domain equations of motion are frequency dependent, and the equations are solved for a discrete set of frequencies. Since for a fixed frequency the equations are linear, the results give the amplitude and phase of the motion components per unit of wave amplitude as a function of the frequency—the motion RAO. • By combining the motion RAO with the sea spectrum, all the statistics of motion can be calculated, and time series can be obtained for time domain simulations, but only describe steady-state motion. • Depending on the hydrodynamic program used, the models may not accurate at low frequency. This is characteristic of strip theory codes (2D potential Theory) [195]. Despite all the apparent limitations of the method, it should be emphasised that the potential flow formulation solved by strip theory gives very good results compared to other, more complex formulations and more importantly, compared to experimental results. See [195, 160] for strip theory formulations, and [18] for details 4.3 Manoeuvring Theory Models 79 about other computational methods for ship motion in seaway. • The linearity assumption limits the validity of the model to waves of small steepness and small motion amplitudes, so superposition can be applied to the potential. 4.3 Manoeuvring Theory Models In manoeuvring theory, as well as in guidance and navigation [67], the equations of motion are described using the variables η defined in the north-eastdown frame (n-frame) and the variables ν defined in the body-fixed frame (b-frame)—see Section 3.3. Using the notation of Fossen [66, 67], the vector equation to describe the motion of a marine vehicle using these variables can be expressed as [MbRB + MbA ]ν̇ + Cb (ν)ν + Db (ν)ν + gb (η) = τ b η̇ = Jnb (Θnb )ν, (4.36) where • MbRB is the rigid-body generalised mass matrix (mass and inertia) with respect to the origin of the b-frame, • MbA is the generalised added-mass matrix, • Cb (ν) is the total (rigid body and added mass) Coriolis and centripetal acceleration matrix, • Db (ν) is the damping matrix, • gb (η) is the restoring function, • Jnb (Θnb ) is the velocity transformation matrix given in (3.24), • τ b is the vector of forces and moments acting on the hull originated by the control devices, the propulsion system, and hydrodynamic effects. Because the motion is described in a non-inertial frame (b-frame), the equations of motion include fictitious accelerations. We will next describe the different elements of equation (4.36). 4.3.1 Rigid Body Dynamics in the b-frame The rigid-body equations of motion expressed in the body-fixed reference frame or b-frame are: MbRB ν̇ + CbRB (ν)ν = τ b , (4.37) where MRB is the generalised mass matrix   mI3×3 −mS(rbg ) MbRB  , mS(rbg ) Ib (4.38) 80 4 Ship Kinetics which is of the form of (4.3). The Coriollis and centripetal acceleration matrix can be expressed in different ways; one representation is   mS(ν 2 ) −mS(ν 2 )S(rbg ) CbRB (ν)  , (4.39) mS(rbg )S(ν 2 ) −S(Ib ν 2 ) where ν 2  [p, q, r]t —see [67] for alternative representations of (4.39). The terms in CRB (ν)ν (4.37) are fictitious forces and moments arising from expressing the equations of motion in the non-inertial b-frame. The components of (4.37) are [66]   m u̇ − vr + wq − xbg (q 2 + r2 ) + ygb (pq − ṙ) + zgb (pr + q̇)   m v̇ − wp + ur − ygb (r2 + p2 ) + zgb (qr − ṗ) + xbg (qp + ṙ)   m ẇ − uq + vp − zgb (p2 + q 2 ) + xbg (rp − q̇) + ygb (rq + ṗ) b b b Ixb ṗ + (Izb − Iyb )qr − (ṙ +pq)Ixz + (r2 − q 2 )Iyz + (pr − q̇)Ixy  b b +m yg (ẇ − uq + vp) − zg (v̇ − wp + ur) b b b Iyb q̇ + (Ixb − Izb )rp − (ṗ + qr)Ixy + (p2 − r2 )Izx + (qp − ṙ)Iyz  +m zgb (u̇ − vr + wq) − xbg (ẇ − uq + vp) b b b Izb ṙ + (Iyb − Ixb )pq − (q̇ +rp)Iyz + (q 2 − p2 )Ixy + (rq − ṗ)Izx  b b +m xg (v̇ − wp + ur) − yg (u̇ − vr + wq) = τ1b = τ2b = τ3b = τ4b = τ5b = τ6b (4.40) For the motion control problem addressed in this book (course keeping and roll stabilisation), it is a common practice to neglect the pitch and heave motion components. This yields a model in four degrees of freedom (4DOF): surge, sway, roll and yaw. Under this assumption, m[u̇ − vr − xbg r2 − ygb ṙ + zgb pr] = τ1b m[v̇ + ur − ygb (r2 + p2 ) − zgb ṗ + xbg ṙ] = τ2b b b b Ixb ṗ − ṙIxz + r2 Iyz + prIxy + m[ygb vp − zgb (v̇ + ur)] = τ4b (4.41) b b b Izb ṙ − rpIyz − p2 Ixy − ṗIzx + m[xbg (v̇ + ur) − ygb (u̇ − vr)] = τ6b . The position of the origin ob of the b-frame can be chosen so as to simplify the equations of motion. For instance, if ob coincides with CG, and the axes xb , yb , zb coincide with the principal axes of inertia; then, the simplest form of the equations of motion (4.41) is obtained. Some disadvantages of this choice, however, are that the axes xb , yb , zb may differ from the symmetry axes of the ship and also that the location of the CG may vary with the loading condition. These effects has to be compensated for by the control system if such equations of motion are used for control applications [66]. As a consequence, it is often more convenient to chose the origin ob such that 4.3 Manoeuvring Theory Models 81 the inertia products are negligible and the axes xb , yb , zb correspond to the longitudinal, lateral and normal direction of the vehicle. As shown by Fossen, [66], this can be achieved by choosing ob such that the coordinates of CG satisfy the following relationships: g 2 g g xg = −Ixy Ixz mIyz g 2 g g mIxz yg = −Ixy Iyz g 2 g g mIyz xg = −Ixz Iyz where the superscript g denotes that the moments of inertia are taken with the body frame fixed at CG. The relationship between moments of inertia with respect to the b-frame and those located at CG are related via b g = Ixx + m[(ygb )2 + (zgb )2 ] Ixx b g Iyy = Iyy + m[(zgb )2 + (zgb )2 ] b g = Izz + m[(xbg )2 + (ygb )2 ]. Izz The following equations of motion are valid when body-fixed axes correspond to the longitudinal, lateral, and normal directions [66]:   m u̇ − vr + wq − xbg (q 2 + r2 ) + ygb (pq − ṙ) + zgb (pr + q̇) = τ1b   m v̇ − wp + ur − ygb (r2 + p2 ) + zgb (qr − ṗ) + xbg (qp + ṙ) = τ2b   m ẇ − uq + vp − zgb (p2 + q 2 ) + xbg (rp − q̇) + ygb (rq + ṗ) = τ3b   b b Ixx ṗ + (Izz − Iyy )qr + m ygb (ẇ − uq + vp) − zgb (v̇ − wp + ur) = τ4b   b b b Iyy q̇ + (Ixx − Izz )rp + m zgb (u̇ − vr + wq) − xbg (ẇ − uq + vp) = τ5b   b b b Izz ṙ + (Iyy − Iyx )pq + m xbg (v̇ − wp + ur) − ygb (u̇ − vr + wq) = τ6b (4.42) and in 4DOF m[u̇ − ygb ṙ − vr − xbg r2 + zgb pr] = τ1b m[v̇ − zgb ṗ + xbg ṙ + ur − ygb (r2 + p2 )] = τ2b b Ixx ṗ − mzgb v̇ + m[ygb vp − zgb ur] = τ4 (4.43) b Izz ṙ + mxbg v̇ − mygb u̇ + m[xbg ur + ygb vr] = τ6b Expression (4.43) is the model that will be adopted for the rigid-body dynamics in the sequel. The total vector of forces and moments τ —appearing in the right-hand side of expression (4.43)—is generated by different phenomena and can be 82 4 Ship Kinetics separated into components according to their originating effects as, and the total effect be studied using superposition [133]: τ b = τ bhyd + τ bhs + τ bc + τ bp , where subscripts stand for • • • • Hydrodynamic forces and moments, Hydrostatic forces and moments, Control devices forces and moments, Propulsion forces and moments. In this book, it will be assumed that τ bp compensates for the hydrodynamic resistance of the hull, and also that the dynamics associated with the surge component of motion are much slower than the dynamics of the other motion components. This assumption allows us to decouple the surge component and to treat the variable u as a constant equal to the ship service speed ū, i.e. in the sequel it will be assumed that u̇ ≈ 0 and u ≈ U . 4.3.2 Manoeuvring Hydrodynamics Due to the complex phenomena associated with the hydrodynamics forces during manoeuvring, these forces motions are usually determined from experimental scaled-model tests. In this approach, The forces are modelled as a general nonlinear functions: τ bhyd = fhyd (ν̇, ν, η). (4.44) These are the radiation and viscous forces at low frequency, and the hydrostatic or restoring forces. The only restoring force relevant to manoeuvring is the roll restoring moment is the roll moment: b τ4hs = GZ(φ)ρg∇, (4.45) where GZ(φ) is the so-called roll righting arm shown in Figure 4.3. If data for GZ(φ) is not available, the linear approximation given in (4.11) is sufficient for surface vessels. The first term in (4.44) is often expanded in a series. There are two approaches to express such a series. The first approach, proposed by Abkowitz, [1], consists in using a truncated Taylor series with only odd terms of third order. The second approach, uses the, so-called second-order modulus terms. This method was proposed by Fedyaevsky and Sobolev [65], and later by Norrbin [162]. The manoeuvring model of the benchmark example (see Appendix B) is given in this latter form [33]: 4.3 Manoeuvring Theory Models 83 Sway terms b τ2hyd = Yv̇ v̇ + Yṙ ṙ + Yṗ ṗ +Y|u|v |U | v + Yur U r + Yv|v| v |v| + Yv|r| v |r| + Yr|v| r |v| +Yφ|uv| φ |U v| + Yφ|ur| φ |U r| + Yφuu φ U 2 . (4.46) Roll terms b τ4hyd = Kv̇ v̇ + Kṗ ṗ +K|u|v |U | v + Kur U r + Kv|v| v |v| + Kv|r| v |r| + Kr|v| r |v| (4.47) +Kφ|uv| φ |U v| + Kφ|ur| φ |U r| + Kφuu φ U 2 + K|u|p |U | p 3 +Kp|p| p |p| + Kp p + Kφφφ φ − ρg∇GZ(φ). Yaw terms b τ6hyd = Nv̇ v̇ + Nṙ ṙ +N|u|v |U | v + N|u|r |U | r + Nr|r| r |r| + Nr|v| r |v| +Nφ|uv| φ |U v| + Nφu|r| φ U |r| + Np p + N|p|p |p|p + N|u|p |U |p +Nφu|u| φ U |U | . (4.48) The linear coefficients in expressions (4.46) to (4.48) are called the hydrodynamic derivatives since, for example, Yṗ = b ∂τ2hyd , ∂ ṗ and Kp = b ∂τ4hyd ∂p are the force in sway due to the roll rate derivative (added mass term), and the roll moment due to the roll rate (damping term). The hydrodynamic derivatives proportional to the accelerations are the added masses and moments of inertia similarly to the seakeeping theory. The difference, however, is that in manoeuvring, these are measured by making the scale model oscillate at low frequencies. The coefficients of the nonlinear terms are simply obtained from curve fitting and therefore should not be called hydrodynamic derivatives. The use of second-order modulus terms have proven to represent the cross-flow drag at large angles of attack [50]. As mentioned by Clarke [50] it has been found that for simulation purposes, it is much more straightforward to store the data and use look-up tables to interpolate rather than to fit explicit curves. For a recent review on the foundations of the manoeuvring equations, see [50]. 4.3.3 Nonlinear Manoeuvring State-space Models In order to perform time domain numerical simulations and to design model based control strategies, it is convenient to use state-space methods. From the models provided in the previous section, the components of vector equation (4.36) in sway, roll and yaw can be expressed as 84 4 Ship Kinetics − (m − Yv̇ )v̇ − (mzG + Yṗ )ṗ + (mxG − Yṙ )v̇ = τ2hyd − mur + τ2c − −(mzG + Kv̇ )v̇ + (Ix − Kṗ )ṗ − Kṙ ṗ = τ4hyd + mzG ur + τ4c − (mxG − Nv̇ )v̇ − Nṗ ṗ + (Iz − Nṙ )ṙ = τ6hyd − mxG ur + τ6c φ̇ = p ψ̇ = r cos(φ) (4.49) − Here, the terms τihyd , i = 2, 4, 6 correspond to the nonlinear hydrodynamic terms, for example, given in Expressions (4.46) to (4.48) without the terms that are proportional to the accelerations—which have been included on the left-hand side of (4.49). The forces and moments produced by the control surfaces, reviewed in Chapter 5, are represented by the terms  b b b t τ4c τ6c . (4.50) τ bc = τ2c Equations (4.49) are already in a form similar to a state-space form. Indeed, by defining the state vector vector  t x vprφψ , then Expression (4.49) can be written as ⎡ −1 ⎤  −1  M M 0 ẋ = f (x) + ⎣ 01×3 ⎦ τ bc , 0 I2×2 01×3 with and ⎤ (m − Yv̇ ) −(mzgb + Yṗ ) (mxbg − Yṙ ) b ⎦ − Kṗ ) −Kṙ M  ⎣−(mzgb + Kv̇ ) (Ixx b b − Nṙ ) −Nṗ (Izz (mxg − Nv̇ ) ⎡ f (x) = fhyd (x) + fc (x). (4.51) (4.52) (4.53) with the obvious definitions for the functions fhyd (x) and fc (x): ⎡ − ⎤ τ2hyd ⎢τ − ⎥ ⎢ 4hyd ⎥ − ⎥ fhyd (x) = ⎢ ⎢τ6hyd ⎥ ⎣ 0 ⎦ 0 ⎡ ⎤ −mur ⎢ mzgb ur ⎥ ⎢ ⎥ b ⎥ fc (x) = ⎢ ⎢−mxg ur⎥ ⎣ ⎦ 0 0 (4.54) 4.3 Manoeuvring Theory Models 85 4.3.4 Linear Manoeuvring State-space Models The model presented in the previous subsection is a comprehensive model that can be used as a calibration model to test different control strategies. In order to design a control system and to draw conclusions about intrinsic limitations, it is convenient to work with a linear model: ẋ = A x + B τ c . (4.55) By taking the first term of the Taylor expansion around the equilibrium point x = 0, u = u, the matrices of this linearised model are defined as: ⎡ −1 ⎤   −1  M ∂f (x)  M 0 ⎣ 01×3 ⎦ , (4.56) A B  0 I2×2 ∂x x=0 01×3 where, for example, by using (4.46)–(4.48)  ∂f (x)  = ∂x x=0 ⎡ Y|u|v |u| 0 (Yur − m)u Yφuu u2 ⎢K|u|v |u| Kp + K|u|p |u| (Kur + mzG )u Kφuu u2 − ρg∇GM t ⎢ ⎢N|u|v |u| Np + N|u|p |u| N|u|r |u| − mxG u Nφu|u| u|u| ⎢ ⎣ 0 1 0 0 0 0 1 0 ⎤ 0 0⎥ ⎥ 0⎥ ⎥. 0⎦ 0 (4.57) Model (4.55) incorporates all the relevant couplings of interest between roll, sway and yaw to analyze the problem of rudder-based stabiliser control system design [33]. This is an important characteristic for studying design performance limitations associated with the dynamics of the ship. This topic will be the subject of the second part of this book. With this, we have completely defined all the elements of the manoeuvring model indicated in Figure 4.1. The following summarises the properties of the manoeuvring models [174]: • The equations of motion are formulated in a reference frame fixed to the ship (b-frame), and not in an equilibrium frame like in seakeeping (h-frame). • The equations can be linear or non-linear depending on the application. For autopilots (course-keeping), and autopilots with rudder rudder roll stabilisation, linear models are usually sufficient, since large deviations requiring non-linear terms can be regarded as poor course-keeping [50]. 86 4 Ship Kinetics • The coefficients of the equations are estimated from captive scale-model tests, by measuring forces while the model is subjected to low frequency forced oscillations in 3DOF (surge, sway and yaw) or 4DOF (with the addition of roll). These coefficients can also be obtained using hydrodynamic programs combined with look-up tables, or by scalling data from known models [208]. • The model is valid in calm water conditions because the coefficients are estimated from data collected from low frequency experiments. Therefore, if there is a significant motion induced by the waves, the memory effects associated with the radiation forces are not accounted for. We finish this section by showing some simulations results obtained using the manoeuvring models presented. Figure 4.13, shows the value of the state variables for the linear and the non-linear model of our benchmark example (see Appendix B) under a zig-zag test [133]. All the coefficients are given in Appendix B. 4.4 A Force-superposition Model for Slow Manoeuvring in a Seaway As commented in Section 4.1, the combination of seakeeping model as an output disturbance for the manoeuvring model is the commonly used approached to model ships for control system design, but this approach suffers from two drawbacks. The first one is that the model may not be used for multibody system interactions, and the second is that there is some uncertainty in the response to the control actions because of miss-modelled dynamics related to changes in added mass and damping with the frequency of the wave excitation. A more realistic description it that shown in Figure 4.2, in which the wave excitation loads are used as input disturbance and a unified model that accounts for memory effects is used. This type of model can therefore be used for slow manoeuvring in a seaway. The slow manoeuvring restriction is because the hydrodynamic frame h is used to obtain the wave-excitation and radiation forces. Thus, it is implicitly assumed that the speed of the manoeuvre is much slower than the wave-induced motion so the h-frame can be considered inertial. This is a reasonable assumption at least for large vessels. In this section, we further describe the elements of such model. The material that follows has been motivated by the work of Bailey et al.[11] and Kristansen and Egeland [127]. 4.4.1 Time Domain Seakeeping Models in the h-frame We have seen that for sinusoidal wave excitation forces, the linear equations of motion in the h-frame can be expressed as [195, 63], 4.4 A Force-superposition Model for Slow Manoeuvring in a Seaway p[deg/s] v[m/s] 0.5 0 −0.5 −1 −1.5 200 250 3 3 2 2 1 0 −1 0 −1 −2 −3 −3 200 250 −4 300 10 200 250 300 200 250 300 20 5 0 −5 ψ[deg] 5 φ[deg] rudder [deg] 1 −2 −4 300 r[deg/s] 1 87 0 −5 −10 −15 200 250 300 −10 10 0 −10 200 250 300 −20 Fig. 4.13. Simulation results corresponding to the naval vessel under an IMO 10 deg zig-zag test. Nonlinear model (solid). Linearised model (dotted). [MhRB + Ah (ωe )]ξ̈ + Bh (ωe )ξ̇ + Gh ξ = τ hw . (4.58) However, as already mentioned, these are not truly equations of motion, because they are valid only if the wave excitation forces are sinusoidal, and provided the coefficients assume an appropriate value according to the frequency of the excitation. Further, these equations describe motion only in steady state. Cummins, in 1962 [58] considered the behaviour of the fluid and the ship in the time domain. He made the assumption of linearity, and considered impulses in the components of motion. This resulted in a boundary value problem in which the potential was separated into two parts; one valid during the duration of the impulses and the other valid after the impulses extinguished. By expressing the pressure as a function of these potentials and integrating it over the wetted surface of the vessel, he obtained a vector integro-differential equation, which is known as the Cummins Equation: 88 4 Ship Kinetics [MhRB h + Ā ]ξ̈ + t  −∞ K̄h (t − τ ) ξ̇(τ ) dτ + Gh ξ = τ hw . (4.59) The matrix Ah is constant and depends only on the ship geometry. The entries of the matrix K̄h (t − τ ), in the convolution integral, are retardation functions of time, which depend on the forward speed and geometry of the vessel; this functions are the impulse response functions of the velocities. The Cummings equation reveal the structure of the true equations of motion of a ship and are valid for any bounded excitation τ hw . Since equation (4.59) is valid for arbitrary excitations, this includes sinusoids as a particular case. Ogilvie, [168], took the Fourier transform of equation (4.59) for sinusoidal excitations, and found the following relationships  ∞ 1 K̄h (τ ) sin(ωe τ )dτ (4.60) Ah (ωe ) = Āh − ωe 0  ∞ K̄h (τ ) cos(ωe τ )dτ (4.61) Bh (ωe )= 0 Since the first equation must be valid for all ωe , it follows that Āh = lim Ah (ωe ) ≡ Ah (∞) ωe →∞ (4.62) The second equation is rewritten using the inverse Fourier transform giving:  2 ∞ h B (ωe ) cos(ωe τ )dωe (4.63) K̄h (t) = π 0 This expression is recognized as a matrix of retardation functions. Thus, the Cummins Equation (4.59) can be expressed as  t K̄h (t − τ )ξ̇(τ )dτ + Gh ξ = τ hw . (4.64) [MhRB + Ah (∞)]ξ̈ + −∞ An alternative representation is [MhRB h h + A (∞)]ξ̈ + B (∞)ξ̇ +  t −∞ where Kh (t) = 2 π  Kh (t − τ )ξ̇(τ )dτ + gh (ξ) = τ hw . (4.65) ∞ [Bh (ωe ) − Bh (∞)] cos(ωτ )dωe , (4.66) 0 The term gh (ξ) indicates that nonlinear restoring forces may be used, instead of the linear approximation Gh ξ. The above form of the equations of motion results in better numerical properties for the convolution integral because [Bh (ωe ) − Bh (∞)] goes to zero as ωe goes to infinity. 4.4 A Force-superposition Model for Slow Manoeuvring in a Seaway 89 4.4.2 Seakeeping Model in the b-frame Using the transformations for the velocities and accelerations (3.60), the equations of motion (4.58) expressed in the b-frame become (Jhb )t [MhRB +A h ˙ (ωe )] [Jhb δν h + U Lδν] + B (ωe )  Jhb δν U ˙ − 2 Lδν ωe  + (Jhb )t gh (ξ) = (Jhb )t τ hw . (4.67) The transformation of the generalised mass matrix from the h-frame to the b-frame is given by Then, ˙ + [Cb + Cb (ωe )]δν + Bb (ωe )δν [MbRB + MbA (ωe )] δν RB A + gb (η) = τ bw , (4.68) where MbRB = (Jhb )t MhRB Jhb , Bb (ωe ) = (Jhb )t Bh (ωe )Jhb MbA (ωe ) = (Jhb )t Ah (ωe )Jhb − U h t h (J ) B (ωe )Jhb , ωe2 b CbRB = U MbRB , (4.69) CbA = U (Jhb )t Ah (ωe )Jhb gb (η) = (Jhb )t gh (ξ), τ bw = (Jhb )t τ hw . Notice that this transformation of (4.58) to the b-frame generates two new matrices CbRB and CbA (ωe ), which are recognized as the Coriolis and centripetal matrices due to rigid-body and frequency dependent added mass. There matrices appear as a consequence of expressing the equations of motion in a non-inertial frame. Let us for convenience define the term Nb (ωe )  CbA (ωe ) + Bb (ωe ), = (Jhb )t [Bh (ωe ) + U Ah (ωe )L]Jhb (4.70) and express (4.68) as ˙ + Nb (ωe )δν + Db (ωe )δν + gb (η) = τ b , [MbRB + MbA (ωe )] δν w (4.71) Therefore, a time-domain representation for the above model valid for any type of excitation is the Cummins equation in the b-frame: 90 4 Ship Kinetics ˙ + Cb δν + Nb (∞)δν + [MbRB + MbA (∞)]δν RB  t Kb (t − τ )δν(τ ) dτ 0 + gb (η) = τ bw , (4.72) where MbRB = (Jhb )t MhRB Jhb , MbA (∞) b = N (∞) = Kb (t) = (Jhb )t Ah (∞)Jhb , (Jhb )t [Bh (∞) + U Ah (∞)L]Jhb ,  2 ∞ b b π (N (ω) − N (∞)) cos(ωt)dω, (4.73) (4.74) (4.75) (4.76) 0 The above is a unified model that can be used for slow manoeuvring in a seaway. Indeed, the classical linear manoeuvring model—linear version of (4.36)—can be obtained from (4.71), by taking the limit ωe → 0. By doing this, we find the relationship between the added masses and potential damping obtained from strip theory and the hydrodynamic derivatives: U MA = (Jhb )t Ah (0)Jhb − lim 2 (Jhb )t Bh (ωe )Jhb ωe →0 ωe ⎤ ⎡ Xu̇ Xv̇ Xẇ Xṗ Xq̇ Xṙ ⎢ Yu̇ Yv̇ Yẇ Yṗ Yq̇ Yṙ ⎥ ⎥ ⎢ ⎢ Zu̇ Zv̇ Zẇ Zṗ Zq̇ Zṙ ⎥ ⎥ ⎢ = −⎢ ⎥ ⎢ Ku̇ Kv̇ Kẇ Kṗ Kq̇ Kṙ ⎥ ⎣ Mu̇ Mv̇ Mẇ Mṗ Mq̇ Mṙ ⎦ Nu̇ Nv̇ Nẇ Nṗ Nq̇ Nṙ Dblin = (Jhb )t [Bh (0) + U Ah (0)L]Jhb ⎤ ⎡ Xu Xv Xw Xp Xq Xr ⎢ Yu Yv Yw Yp Yq Yr ⎥ ⎥ ⎢ ⎢ Zu Zv Zw Zp Zq Zr ⎥ ⎥ , = −⎢ ⎢ Ku Kv Kw Kp Kq Kr ⎥ ⎥ ⎢ ⎣ Mu Mv Mw Mp Mq Mr ⎦ Nu Nvp Nw Np Nq Nr where Dblin ν is the linear part of Db (ν)ν in (4.36). It should be noticed, however, that if coefficients in Ah (ωe ) and Bh (ωe ) calculated using strip-theory codes (2D potential theory), they are not accurate as ωe → 0, because this is not contemplated in the assumptions used to derive strip theory—see, for example, [195] for details. 4.4 A Force-superposition Model for Slow Manoeuvring in a Seaway 91 4.4.3 A Unified Nonlinear State-pace Model Kristiansen and Egeland, [127], have proposed a state-space formulation for the convolution term in (4.72). If we consider: γ(t) =  t Kb (t − τ )δν(τ )dτ causal =  t Kb (t − τ )δν(τ )dτ ; (4.77) 0 −∞ Then, for causal systems it follows that: Kb (t − τ ) = 0 for t < 0 (4.78) If δν as a unit impulse, then (4.77) will be an impulse response function. Consequently, γ(t) can be represented by a linear state-space model: µ̇ = Abr µ + Bbr δν, γ= Cbr µ + Dbr δν µ(0) = 0 (4.79) (4.80) where (Abr , Bbr , Cbr , Dbr ) are constant matrices of appropriate dimensions— see [127] for details on how matrices of the state-space can be obtained. As stated in [68], a good approximation to the convolution term is obtained using typically 5 states per entry of the matrix Kb . For example, if we would like to design an autopilot with this model, we could consider the sway velocity, b b the yaw angle and the entries K66 , and K22 . Thus the control design model will have 13 states. Using these results, we can obtain a state-space representation in the bframe that is valid for slow manoeuvring (slow turning) in a seaway with forward speed : [MbRB + MbA (∞)]ν̇ = −(Nb (∞) + Dbr + CbRB )(ν − ν̄) − Cbr µ − gb (η) + τ bW + τ bP , µ̇ = Abr µ + Bbr (ν − ν̄), η̇ = Jnb (Θhb )ν, (4.81) The term τ bP represents the propulsion forces that take the ship to the speed ν̄. In the first equation in (4.81), we can now add the control action forces τ bC and slowly-varying environmental forces τ bE (due to wind, current and second order wave excitation loads) and additional viscous damping terms DbV (ν)ν. The second equation in (4.81) provides additional dynamics associated with the radiation potential, i.e. fluid memory effects. The following summarises the characteristics of the model given above: 92 4 Ship Kinetics • Compared to the motion superposition model shown in Figure 4.1, the model (4.81) provides a more accurate calculation of the radiation forces when manoeuvering in a seaway. Indeed, in the motion superposition model, the radiation forces are computed using only the motion induced by the waves (this is embedded in the motion RAO), whereas in the force superposition model (4.81), the radiation forces are computed using the motion induced by the combination of the wave and the control forces. This is well known in marine technology, but the state-space representation proposed in [127], and the formulation of the model in the body-fixed frame presented in this chapter, allows one to use the model for control system design. To the best of the authors’ knowledge, these high-order models have never before been used for model-based control system design . • The model uses data readily available from standard seakeeping programs. This is an advantage to the designer of motion control systems because a preliminary control design and evaluation can be performed before doing system identification. The only data required by these programs are the drawing of the hull and the loading condition. • Different corrections to incorporate viscous effects can be included either in the time domain or in the frequency domain (by modifying Bh (ωe )) before calculating retardation functions [11, 68, 69]. Also, if the potential theory data are complemented by data from a PMM test, the combination of these can contribute to improve the quality of the model, especially at low frequency. Details about this can be seen in [11, 68]. • A limitation of the model is that it is valid for slow manoeuvring (slow turning). The reason for this is that the wave-excitation forces and the retardation functions are obtained using programs based on seakeeping theory. Therefore, the model is limited to slow turning so that the h-frame can be considered approximately inertial. • Other limitations of the model are heavily dependent on the seakeeping software used. These are related to the type of ship, the theory, and the code implementation. Therefore, one should be aware of these limitaions. For example, if a strip theory software is used (2D potential theory), there may be limitations in the accuracy of the data at low encounter frequency, and the maximum Froude number should be limited to 0.3, which means that the speed of the vessel should be restricted to $ (4.82) U < 0.3 gLpp.