Simulation of Submarine Manoeuvring Sebastian Thun´e A thesis presented for the degree of Master of Science

(1)

Simulation of Submarine Manoeuvring

Sebastian Thun´

e

A thesis presented for the degree of

Master of Science

(2)

(3)

Abstract

The best way to save money in a submarine design project is to, in an as early stage as possible, be able to find errors in a design or construction. A powerful method to find such errors is to put the design through rigorous testing in simulation software. This project has created a simulation program for a submarines cruising fully submerged. The software is using well-known and established equations of motion (EOM’s) that has been developed through the years. These equations describe the forces and moments acting on the submarine in six degrees of freedom and can be described with Newton’s second law of motion. The equations accelerations are integrated using the scipy.integrate.odeint module in Python to obtain the translational and rotational velocities for a time series. By integrating the EOM’s with different initial values and conditions a library of standard manoeuvring tests have been created together with graphically tailored representation for each test result. The software has been verified by comparing a few chosen test results with simulation software currently used by FOI. Through this comparison it is seen that the software correlates well with the current one in every aspect except one. There is an uncorrelated phenomenon when stern-plane diving is

simulated. The error is identified to the stern-plane drag term in the surge equation, Xδsδs. Evidence shows

(4)

Sammanfattning

Simulering av ub˚

atsman¨

ovrering

Det bästa sättet att spara pengar inom ett ub˚atdesignprojekt är att, i ett s˚a tidigt skede som möjligt, hitta fel

i designen. En kraftfull metod att hitta s˚adana fel är att föra designen genom rigorösa tester i ett simulering

program. Detta projekt har skapat ett s˚adant program för en ub˚at som rör sig i undervattensläge.

Pro-grammet använder sig av välkända och etablerade rörelseekvationer som har utvecklats genom ˚aren. Dessa

ekvationer beskriver krafter och moment som ¨ar verksamma p˚a en ub˚at i sex frihetsgrader som kan beskrivas

med Newton’s andra r¨orelselag. Ekvationernas accelerationer integreras med hj¨alp av scipy.integrate.odeint

modulen i Python f¨or att erh˚alla ub˚atens translation- och rotations-hastigheter f¨or en tidsserie. Genom att

integrera r¨orelseekvationerna med olika begynnelsev¨arden och villkor s˚a har ett bibliotek av man¨

ovrering-stester blivit framtaget tillsammans med anpassad grafisk representation av testresultaten. Programmet har

verifierats genom jämförelse av n˚agra valda test gentemot programmet som i dagsläget används av FOI.

Denna j¨amf¨orelse har visat att programmen korrelerar bra i alla fall utom ett. Vid dykning genom att

anv¨anda aktra djuproder s˚a korrelerar inte programmen gentemot varandra. Felet har blivit identifierat

till dragmotst˚andet som uppst˚ar p˚a rodren i x-ekvationen, Xδsδs. Det ¨ar inte otroligt att felet kan ligga

i programmet som används i dagsläget och en dialog med skaparna har initierats. Tills detta problem är

utrett bör man vara medveten om att resultat fr˚an simuleringar där aktra djuproder används kan ha lägre

(5)

Acknowledgement

The author would like to offer his deepest gratitude to FOI for the opportunity to perform this thesis at their institution, and allowing me to test and improve my engineering skills.

A special thanks goes out to my supervisor and Senior Scientist Lennart Boss´er for his endless support,

discussions and guidance, which made my time at FOI very pleasant and fun. I would also like to thank him for believing that I was capable to fulfil the task at hand. He has been a great mentor during my time at FOI.

A thanks goes out PhD Mats Nordin and to my friend and classmate MsD Johan Fridstr¨om for their inputs

and support during the thesis. Johan should have praise for not going crazy while sharing an office with me. I would also like to thank him for suggesting me for the task to FOI.

Naturally I would like to thank my examiner PhD Ivan Stenius for the discussions, encouragement and enthusiastic support.

Finally I would like to thank the entire underwater department at FOI for all the intriguing small talks at coffee and lunch breaks. The knowledge within the facility is truly amazing and it was sincerely interesting to hear all the discussions that took place.

Stockholm June 2015

(6)

Nomenclature

¯

vFW Lateral velocity due to lift on the bridge fairwater [m/s]

β Drift angle [rad]

δ1 First stern hydroplane deflection [rad]

δ2 Second stern hydroplane deflection [rad]

δ3 Third stern hydroplane deflection [rad]

δ4 Fourth stern hydroplane deflection [rad]

δb Deflection of bow-plane [rad]

δr Deflection of rudder [rad]

δs Deflection of stern-plane [rad]

η The ratio uc

u [ - ]

φ Roll angle [rad]

ψ Yaw angle [rad]

ρ Density of water [kg/m3_]

τ (x) Time interval for vortex [s]

θ Pitch angle [rad]

ξ Bottom clearing b_c [ - ]

B Buoyancy force [N]

b Beam of the vessel [m]

C Revolution scaling ratio [ - ]

c Clearing between baseline and bottom [m]

CD Drag coefficient [ - ]

CL Lift coefficient [ - ]

D Propeller diameter [m]

g Gravitational acceleration [kg/s2_]

(9)

Iyy Moment of inertia about y-axis [kgm2]

Iyz Product of inertia with respect to y- and z-axis [kgm2]

Izz Moment of inertia about z-axis [kgm2]

J Advance ratio [ - ]

K Moments about the x-axis [Nm]

KQ Torque coefficient [ - ]

KT Thrust coefficient [ - ]

l Length of the vessel [m]

M Moments about the y-axis [Nm]

m The vessel mass [kg]

N Moments about the z-axis [Nm]

n Propeller revolution velocity [rpm]

p Angular velocity about x-axis [rad/s]

Q Propeller torque [Nm]

q Angular velocity about y-axis [rad/s]

r Angular velocity about z-axis [rad/s]

T Propeller thrust [N]

tdf Thrust deduction factor [ - ]

u Linear velocity in x-direction [m/s]

uc Commanded velocity [m/s]

v Linear velocity in y-direction [m/s]

v(x) Velocity component v at any x coordinate [m/s]

W Weight of vessel [N]

w Linear velocity in z-direction [m/s]

w(x) Velocity component w at any x coordinate [m/s]

wf Wake fraction [ - ]

X Force in x-direction [N]

x1 Vortex starting point [m]

x2 Vortex detachment point [m]

xB x coordinate for centre of buoyancy [m]

xG x coordinate for centre of gravity [m]

Y Force in y-direction [N]

(10)

yG y coordinate for centre of gravity [m]

Z Force in z-direction [N]

zB z coordinate for centre of buoyancy [m]

zG z coordinate for centre of gravity [m]

Abbreviations

CB Centre of Buoyancy

CG Centre of Gravity

DOF Degrees of Freedom

EOM Equations of Motion

FOI Totalf¨orsvarets forskningsinstitut (Swedish defence research agency)

HS Hydrostatics

HD Hydrodynamics

(11)

Chapter 1

Introduction

A quick Google search of the sentence ”definition of submarine” gave the short but rather good answer: ”A nautical vessel capable of operating submerged”. In other words, a submarine is a manned atmospheric multi-purpose vessel that removed the constraints of just moving in the 2D horizontal plane of the vast ocean to the 3D ocean space in a controlled way. The idea of submarines can be traced all the way back to the Renaissance, where sketches of an underwater vessel drawn by nobody else then the universal-genius Leonardo Da Vinci has been found. The first practical submarine designed and manufactured was, however, done 1620 by the Dutch inventor Cornelis van Drebbel. This vessel was propelled by paddles. Today, the submarines are considered one of the most complex machines that exist, but seldom seen as they manoeuvre the dark oceans beneath.

The Swedish history of submarine design begins with the development and building of Nordenfelt I in 1883. Based on an original idea from the English reverend George Garrett, the Swedish inventor and entrepreneur Thorsten Nordenfelt built the Swedish first machine driven submarine at AB Palmcrantz & Co at Karlsvik on Kungsholmen, Stockholm. The steam engine came from J. & C. G. Bolinders Mekaniska Verkstads AB

at lake Klara at Kungsholmen, Fleminggatan 4, in Stockholm where today the Tekniska n¨amndhuset is

situated. The submarine was finally assembled at Ekensbergs yard AB, Kungsholmen, Figure 1.1.

(12)

During 1884 the Nordenfelt I was further refined at Eriksbergs yard in G¨oteborg. The submarine was finally shown for a collection of high-level individuals and royalties from Europe, just north of Landskrona in 1885. In spite of a successful demonstration interest among the nobility did not result in any further business. The submarine was later sold to Greece after being inspected by the Royal Swedish Navy with the comments that the interior of the submarine was too hot.

The Swedish navy history of submarine design begins with HMUb Hajen, Figure 1.2, launched July 16, 1904. The vessel was designed based on blue-prints that the Swedish marine engineer, Carl Richson had got his hands on during his apprentice period at Mr. Hollands submarine yard in America 1900. After some redesign and royal acceptance by the Swedish King Oscar II, the submarine was built at the Naval yard in Stockholm 1902-04. After a one year test period performed by the KMF (later FMV), and performed by the former chief engineer of Nordenfelt I, HMUb Hajen became part of the Swedish navy.

Figure 1.2: HMUb Hajen on display in Karlskrona.

In 1945 the German submarine U3503 sank in Swedish territorial waters and was later salvaged by Sweden in 1946[13]. This was the starting point of a huge initiative for national development of submarines within the Swedish navy , Swedish industry and governmental agencies and technical university’s. A little bit more

than twenty years later, 1967, Sweden launched their new model, Sj¨oormen (project A11), again influenced

by American submarine development, with a new hull-form giving the submarine superior manoeuvrability compared to the older models. After this, Sweden continued to develop its own submarines and followed

Sjöormen with three more classes: Näcken (project A14), Västergötland (project A17) and Gotland (project

A19). Three out of the five submarines that are in use today in the Swedish navy are of the Gotland class.

The other two submarines are the modernisation of the V¨asterg¨otland class submarines that was updated

and re-launched as the S¨odermanland class[8]. This class will be replaced by a new type, project A26,

launching in 2019-2021.

(13)

(14)

Chapter 2

Thesis structure and goals

This thesis has two goals. The first is to create the core functions of a submarine manoeuvring simulation software that can simulate some of the common submarine standard manoeuvring test. The second is to verify the program simulation results against the currently used software.

(15)

Chapter 3

Limitations

The programs biggest limitation is that it is only applicable on submarines fully submerged. There is how-ever room to evolve the program into handling wave interactions, surface propulsion hydrodynamics and hydrostatics.

The program cannot handle the use and manipulation of the ballast tanks to use as trim water or compen-sation water or how the free liquid movements will affect the submarine.

Submarines today are dependent of a control-system. However, in the current version of this program, there is no autopilot available. This is a limitation for some of the tests that should benefit from the option of using an autopilot.

The program can only handle forward motion of the submarine. To simulate the submarine moving in a backward motion, the backward EOM’s have to be implemented to the program in similar fashion as de-scribed in this thesis.

(16)

Chapter 4

Equations of motion structure

Modelling a submarines motion involves a study of both statics and dynamics. Statics are considered when the body is in rest or moving at a constant velocity, the dynamics are concerned with the bodies accelerations. A submarine moving through the oceans has six DOF (degrees of freedom), giving six different components. These six components are conventionally defined as: surge, sway, heave, roll, pitch and yaw. The notations used for the vessel in this report are based on Fossen (1994)[6] and are shown in Table 4.1.

Table 4.1: Notations for the six DOF forces/moments, linear/angular velocities, position and Euler angles.

DOF forces &

moments

linear & angular vel.

position & Euler angles

1 motion in the x-direction (surge) X u x

2 motion in the y-direction (sway) Y v y

3 motion in the z-direction (heave) Z w z

4 rotation about the x-axis (roll) K p φ

5 rotation about the y-axis (pitch) M q θ

6 rotation about the z-axis (yaw) N r ψ

The first three coordinates and their corresponding time derivatives in Table 4.1 are the position and trans-lation motions along x, y, and z-axis. The last three corresponds to the orientation angles and rotational motions, respectively. When analysing the motions it is convenient to define two coordinate systems: The body-fixed and the earth-fixed reference frames. The body-fixed frame is fixed to the moving vessel in its

origin O, often chosen to coincide with the vessels center of gravity, CG. The body axes X0, Y0 and Z0 are

(17)

Figure 4.1: The body- and earth-fixed reference frames.

(18)

J =         c(ψ)c(θ) −s(ψ)c(φ) + c(ψ)s(θ)s(φ) s(ψ)s(φ) + c(ψ)c(φ)s(θ) 0 0 0 s(ψ)c(θ) c(ψ)c(φ) + s(φ)s(θ)s(ψ) −c(ψ)s(φ) + s(θ)s(ψ)c(φ) 0 0 0 −s(θ) c(θ)s(φ) c(θ)c(φ) 0 0 0 0 0 0 1 s(φ)t(θ) c(φ)t(θ) 0 0 0 0 c(φ) −s(φ) 0 0 0 0 s(φ)/c(θ) c(φ)/c(θ)         (4.1)

where s(·) = sin(·), c(·) = cos(·) and t(·) = tan(·). For the deriving of the transformation matrix in detail see Fossen (1994)[6]. The earth-fixed velocities will then be transformed according to Equation 4.2

˙

η = J ν, (4.2)

where ˙η is the time derivative of the positions and Euler angles. ν is the translation and rotational motions

according to Equation 4.3 and 4.4

η = x y z φ θ ψ T, (4.3)

ν = u v w p q r T

. (4.4)

This transformation is vital for extraction of the position and movements of a moving vessel as it is described mathematically with the six EOM’s. The following sections describes the structure of the derived EOM’s and the theory behind their terms such as hydrostatics, hydrodynamics and rigid-body dynamics. Through this chapter the dynamic equation of motion expression is used repeatedly and is, for this report, written as

M ˙ν + C(ν)ν + D(ν)ν + g(η) = τ , (4.5)

where M is the inertia matrix including the added mass, C(ν) is the Coriolis and centripetal terms including the added mass, D(ν) is the damping matrix, g(η) is the vector of gravitational and buoyant forces and moments, also know as hydrostatics, and τ is the vector of control inputs.

4.1 Rigid body dynamics

The rigid body dynamics are derived through application of Newton’s second law of motion that relates the

mass m, and the acceleration ˙ν, to the force fCas

fC= ˙νm. (4.6)

Euler’s first and second axioms, states that Newton’s second law of motion can be expressed in terms of

conservation of the linear and angular momentums, pCand hC respectively as

pC= mvC,

hC= ICω,

(4.7)

where vC, ICand ω is the linear velocity, inertia tensor and angular velocity about the body’s CG,

respec-tively. The time derivatives of Equation 4.7 yields the forces and moments at the body’s CG, fC and mC

(19)

f0= m( ˙v0+ ω × v0+ ˙ω × rG+ ω×(ω × rG))), (4.9)

where rGis the distance from the origin to the centre of gravity expressed in the body coordinates [xG, yG, zG].The

rotational motions derived from the second axiom is

m0= I0ω + ω×(I˙ 0ω) + mrG× ˙v0+ ω × v0, (4.10)

where I0is the body’s inertia tensor defined as:

I0=   Ixx −Ixy −Ixz −Iyx Iyy −Iyz −Izx −Izy Izz  . (4.11)

For complete deduction of the expressions in Equations 4.9 and 4.10 see Fossen[6]. The other notations are written in component form as:

f0= τ1= [X, Y, Z]T External Forces

m0= τ2= [K, M, N ]T Moments from external forces

v0= ν1= [u, v, w]T Linear velocities

ω = ν2= [p, q, r]T Angular velocities

Input of these components into 4.9 and 4.10 gives six equations that can be broken down into a vectorial representation in form of Equation 4.5 as

MRBν + C˙ RB(ν)ν = τRB, (4.12)

where τRB is a sum of all the external forces and moments acting on the submarine further explained in

section 4.7. MRBis the rigid-body inertia matrix:

MRB=         m 0 0 0 mzG −myG 0 m 0 −mzG 0 mxG 0 0 m myG −mxG 0 0 −mzG myG Ixx −Ixy −Ixz mzG 0 −mxG −Iyx Iyy −Iyz

−myG mxG 0 −Izx −Izy Izz

        , (4.13)

and where the rigid-body Coriolis and centripetal matrix CRB(ν) is

CRB(ν) =     0 0 0 m(yGq+zGr) −m(xGq−w) −m(xGr+v)

0 0 0 −m(yGp+w) m(zGr+xGp) −m(yGr−u)

0 0 0 −m(zGp−v) −m(zGq+u) m(xGp+yGq)

m(yGq+zGr) −m(xGq−w) −m(xGr+v) 0 −Iyxq−Ixzp+Izzr Iyzr+Ixyp−Iyyq

−m(yGp+w) m(zGr+xGp) −m(yGr−u) −Iyzq+Ixzp−Izzr 0 −Ixzr−Ixyq+Ixxp

−m(zGp−v) −m(zGq+u) m(xGp+yGq) −Iyzr−Ixyp+Iyyq Ixzr+Ixyq−Ixxp 0

    . (4.14)

4.2 Hydromechanics

The hydromechanics can, as for the rigid-body, be described with Equation 4.5 where the sum of the hydromechanic forces and moments can be written as

τH = −MAν − C˙ A(ν)ν − D(ν)ν − g(η), (4.15)

where MAand CA(ν) are the hydrodynamic added inertia and Coriolis and centripetal terms, respectively.

g(η) is the hydrostatic term presented in section 4.2.3. D(ν) is the combined damping matrix consisting

of the potential damping DP(ν), the skin friction DS(ν), the wave drift damping DW(ν) and the vortex

shedding damping DM(ν) written as

(20)

4.2.1 Added mass

Added mass should be understood as the the pressure-induced forces and moments due to a forced harmonic motion of the body, which are proportional to the acceleration of the body. In other terms, as the body is performing a motion in the fluid it is pushing a fluid mass as it is moving through it, generating forces and moments. If the vehicle is completely submerged, the added mass coefficients will be constant and

independent of wave frequency. To derive the expression of the added inertia matrix, MA, and the the

matrix of the hydrodynamic Coriolis and centripetal terms, CA(ν), the fluid kinetic energy approach in

terms of Kirchhoff’s equation is used. Fluid kinetic energy is the energy needed for the fluid to move aside and close behind the moving vessel. The expression of fluid kinetic energy may be expressed as[12]

TA=

1

2ν

T_M

Aν, (4.17)

where the added inertia matrix, MA, is defined as a 6x6 matrix

MA=         Xu˙ Xv˙ Xw˙ Xp˙ Xq˙ X˙r Yu˙ Yv˙ Yw˙ Yp˙ Yq˙ Y˙r Zu˙ Zv˙ Zw˙ Zp˙ Zq˙ Z˙r Ku˙ Kv˙ Kw˙ Kp˙ Kq˙ K˙r Mu˙ Mv˙ Mw˙ Mp˙ Mq˙ M˙r Nu˙ Nv˙ Nw˙ Np˙ Nq˙ N˙r         . (4.18)

The matrix expressions follow the SNAME (1950) notation[16]. As an example; the hydrodynamic added

mass force XAalong the x-axis due to the acceleration ˙u in the x-direction is written as XA= Xu˙u. Another˙

common way to express the notation is Aij = −{MA}ij. This gives the notation A11 = −Xu˙ for the same

example as above. However, bodies that are used for underwater vehicles are usually symmetric in multiple planes. With a dorsal fin at the maximum diameter of the body and four fins located at the rear end body, a finned prolate spheroid is a good approximation for a submarine body. This would then give the reduced expression[10][9] MA=         Xu˙ 0 0 0 0 0 0 Yv˙ 0 Yp˙ 0 Y˙r 0 0 Zw˙ 0 Zq˙ 0 0 Kv˙ 0 Kp˙ 0 K˙r 0 0 Mw˙ 0 Mq˙ 0 0 Nv˙ 0 Np˙ 0 N˙r         . (4.19)

The added force and moments can be then be derived through potential flow theory, which describes the velocity field as a gradient of a scalar function. This method is based on the assumptions that the fluid is inviscid, has no circulation and that the body is completely submerged in an unbounded fluid. However, these assumptions are not valid near submerged objects or at the surface. To solve this problem double-body

theory[6] is applied with the assumption that MA= MAT, which gives

2TA= − Xu˙u2− Yv˙v2− Zw˙w2− Kp˙p2− Mq˙q2− N˙rr2

− 2K˙rrp − 2p(Yp˙v) − 2q(Zq˙w) − 2r(Y˙rv). (4.20)

Through the kinetic energy of the fluid, the forces and moments are derived through the application of

Kirchhoff’s equation. This equation relates the fluids energy forces and moments acting on the vessel.

(21)

d dt ∂TA ∂u = r ∂TA ∂v − q ∂TA ∂w − XA, d dt ∂TA ∂v = p ∂TA ∂w − r ∂TA ∂u − YA, d dt ∂TA ∂w = q ∂TA ∂u − p ∂TA ∂v − ZA, d dt ∂TA ∂p = w ∂TA ∂v − v ∂TA ∂w + r ∂TA ∂q − q ∂TA ∂r − KA, d dt ∂TA ∂q = u ∂TA ∂w − w ∂TA ∂u + p ∂TA ∂r − r ∂TA ∂p − MA, d dt ∂TA ∂r = v ∂TA ∂u − u ∂TA ∂v + q ∂TA ∂p − p ∂TA ∂q − NA. (4.21)

Substituting Equation 4.20 into Equation 4.21 will gives the following added mass terms

XA= Xu˙u + Z˙ w˙wq + Zq˙q2− Yv˙vr − Yp˙pr − Y˙rr2,

YA= Yv˙˙v + Yp˙p + Y˙ ˙r˙r − Zw˙wp − Zq˙pq + Xu˙ur,

ZA= Zw˙w + Z˙ q˙q − X˙ u˙uq + Yv˙vp + Yp˙p2+ Y˙rrp,

KA= Yp˙˙v + Kp˙p + K˙ ˙r˙r − vw(Yv˙ − Zw˙) − wr(Y˙r− Zq˙) − Yp˙wp − vq(Y˙r− Zq˙) + K˙rpq − qr(Mq˙− N˙r),

MA= Zq˙( ˙w − uq) + Mq˙q − wu(Z˙ w˙ − Xu˙) + Yp˙vr − Y˙rvp − K˙r(p2− r2) + rp(Yp˙ − N˙r),

NA= Y˙r˙v + K˙rp + N˙ ˙r˙r − uv(Xu˙ − Yv˙) + Yp˙up + Y˙rur + Z˙rur − Yp˙vq − pq(Kp˙ − Mq˙) − K˙rqr. (4.22)

These six equations can be written in component form according to Equation 4.15, where MA is according

to 4.19 and the hydrodynamic Coriolis and centripetal matrix can thus be written as[6]

CA(ν) =         0 0 0 0 −a3 a2 0 0 0 a3 0 −a1 0 0 0 −a2 a1 0 0 −a3 a2 0 −b3 b2 a3 0 −a1 b3 0 −b1 −a2 a1 0 −b2 b1 0         , (4.23)

with the matrix elements

a1= Xu˙u, a2= Yv˙v + Yp˙p + Y˙rr, a3= Zw˙w + Zq˙q, b1= Yp˙v + Kp˙p + K˙rr, b2= Zq˙w + Mq˙q, b3= Y˙rv + K˙rp + N˙rr. (4.24)

(22)

Table 4.2: Cross-terms for added mass couplings derived from added mass terms in the CA(ν) matrix. X Y Z K M N Zw˙ = Xwq −Zw˙ = Ywp −Xu˙ = Zq −(Yv˙ − Zw˙) = Kvw −Zq˙ = Mq Y˙r= Nr Zq˙ = Xqq −Zq˙ = Ypq Yv˙ = Zvp −(Y˙r+ Zq˙) = Kwr −(Zw˙ − Xu˙) = Mw −(Xu˙ − Yv˙) = Nv −Yp˙ = Xpr Xu˙ = Yr Yp˙ = Zpp −Yp˙ = Kwp Yp˙ = Mvr Yp˙ = Np −Y˙r= Xrr Y˙r= Zrp (Y˙r+ Zq˙) = Kvq −Y˙r= Mvp Zq˙ = Nwp −Yv˙ = Xvr K˙r= Kpq −K˙r= Mpp −(Kp˙ − Mq˙) = Npq −(Mq˙ − N˙r) = Kqr K˙r= Mrr −K˙r= Nqr (Kp˙ − N˙r) = Mrp

The added mass terms can be calculated with both numerical and semi-empirical methods, which would be Strip Theory and Planar Motion Mechanism respectively. However, this thesis will not cover that part since all hydrodynamic coefficients and derivatives are given for the task. A description of both the numerical and the semi-empirical method is described in Appendix A.

4.2.2 Damping

In the general case, see Equation 4.16, four terms give contributions to the damping. However, for a fully submerged vessel, only two of the terms are effective:

• The quadratic and non-linear skin friction, DS(ν).

• The viscous damping due to vortex shedding, DM(ν).

This will reduce the damping term to D(ν) = DS(ν) + DM(ν), which are the viscous terms. Generally,

the motion of an underwater vessel moving in six DOF will be highly non-linear and coupled. Therefore, the damping is considered non-linear as well and is described as a third order function. For an example the damping in the x-direction for a velocity u is expressed as[11]

XD= Xuu + Xu|u|u|u| + Xuuuu3. (4.25)

In practice, see Fossen[6] and Ridley & Fontan[15], the third order term is often simplified down to the second order. The third order term could be used if motions are highly non-linear, which occurs due to the speed dependence for Reynolds number. Like the added mass terms, the damping terms can be calculated either semi-empirically or numerically with the strip theory through Morison’s equation[6]

F (U ) = −1

2ρCDAU |U |, (4.26)

where U is the velocity of the vehicle, A is the projected cross-sectional area and CD is the drag coefficient.

4.2.3 Hydrostatics

The hydrostatics term g(η) is the balance between the gravitational and buoyancy forces as explained by Archimedes, also known as the restoring forces. The gravitational force W is acting through the vessels

CG, with its coordinates rG = [xG, yG, zG]. In the same way the buoyancy force B acts in CB, with its

(23)

(a) xz-plane. (b) yz-plane.

Figure 4.2: Gravity and buoyancy forces acting on a submarine.

Resolving the forces onto the submarines body frame with the Euler angles gives the restoring forces and moments in the six DOF as[6]

g(η) =         (W − B) sin(θ) −(W − B)) cos(θ) sin(φ) −(W − B) cos(θ) cos(φ)

−(yGW − yBB) cos(θ) cos(φ) + (zGW − zBB) cos(θ) sin(φ)

(zGW − zBB) sin(θ) + (xGW − xBB) cos(θ) cos(φ)

−(xGW − xBB) cos(θ) sin(φ) − (yGW − yBB) sin(θ)

        . (4.27)

4.3 Non-linear cross flow drag

The manoeuvring characteristics of a submarine depend on, as shown in previous sections, the hydrodynamic forces and moments that are developed at the hull. In a turning manoeuvre the hydroplane lift force will create a moment around the submarines centre axis and thus develop a yawing angular velocity. As the submarine is turning it will also move with a drift angle β. This will create a lifting force on the hull, pushing the submarine towards the centre of the turn, Figure 4.3.

Figure 4.3: Projected turning radius of a submarine and its drift angle, β.

(24)

with a large drift angle the term become non-linear. The non-linear cross flow may, according to Feldman [5] and Bohlmann [2], be calculated using strip theory i.e. the hull is split into multiple sections where the drag force is a function of the local velocity. The total forces are then calculated by integrating over the length of the vessel

YCF= − ρ 2Cd Z l h(x)v(x)pw(x)2_{+ v(x)}2 _dx, ZCF= − ρ 2Cd Z l b(x)w(x)pw(x)2_{+ v(x)}2_dx, MCF= ρ 2Cd Z l xb(x)w(x)pw(x)2_{+ v(x)}2 _dx, NCF= − ρ 2Cd Z l xh(x)v(x)pw(x)2_{+ v(x)}2 _dx, _(4.28)

where x is the axial location on the hull, Cdis the drag coefficient, h(x) and b(x) are the height and beam

respectively for each strip. The local cross-flow velocities, v(x) and w(x), at each strip are calculated as

v(x) = v + xr,

w(x) = w − xq. (4.29)

This model is assuming that the drag coefficient is constant over the whole length (and direction) of the body. However, there are models that use separate drag coefficients for better cross-flow predictions. One of the models is[6]

YCF = − ρ 2 Z l Cdyh(x)v(x)2+ Cdzb(x)w(x)2 v(x) Ucf(x) dx, ZCF = − ρ 2 Z l Cdyh(x)v(x)2+ Cdzb(x)w(x)2 w(x) Ucf(x) dx, MCF = ρ 2 Z l Cdyh(x)v(x)2+ Cdzb(x)w(x)2 w(x) Ucf(x) x dx, NCF = − ρ 2 Z l Cdyh(x)v(x)2+ Cdzb(x)w(x)2 v(x) Ucf(x) x dx, (4.30)

where Ucf(x) =pv(x)2+ w(x)2 is the cross-flow velocity.

4.4 Horse-shoe vortex

On the submarine’s slender body it carries a couple of appendages. The shapes on these independent

(25)

Figure 4.4: The horse-shoe vortex at the sail-body junction[14].

As the submarine turns, lift is created on the bridge fairwater causing this vorticity of varying strength to shed and convect downstream along the hull as a trailing vortex. This vorticity causes circulation around the hull and therefore also lift on the hull. This lift at a location after the sail is considered a product of

the sail’s lift coefficient, Cl, the local lateral velocity component, w(x) or v(x), and the local lateral velocity,

¯

vFW(t − τ [x]), which was developed at the sail at the time the vortex was shed from the sail. With other

words, the lateral velocity created at the sails quarter chord travels downstream along the hull and thus, at a location on the hull, the vortex velocity will be velocity created at the sail a certain time before, Figure 4.5.

Figure 4.5: Illustration of vortex starting point x1with the lateral velocity at the time t and the detachment

point x2with lateral velocity at the time t − τ [x2].

The time at a location x from the vortex starting point x1, can be calculated according to Feldman[4] as

Z t−τ (x)

t

u(t)dt = x1− x. (4.31)

(26)

YV= − ρ 2Cl Z x1 x2 w(x)¯vFW(t − τ [x]) dx, ZV= − ρ 2Cl Z x1 x2 v(x)¯vFW(t − τ [x]) dx, KV= KiuvFW(t − τT) + ρ 2l 2 z1Cl Z x1 x2 w(x)¯vFW(t − τ [x]) dx, MV= − ρ 2lCl Z x1 x2 xw(x)¯vFW(t − τ [x]) dx, NV= − ρ 2lCl Z x1 x2 xw(x)¯vFW(t − τ [x]) dx. (4.32)

In the KV term the Ki contribution is the interaction between the trailing vortex and the stern control

surfaces. The lateral velocity components are calculated as mentioned in Equation 4.29 and the velocity of the hull-bound vortex is calculated as[4]

¯

vFW(t − τ [x]) = v(x) + x1r − zF Wp(x), (4.33)

where zFW is the z-coordinate of the 42% span of the bridge fairwater.

4.5 Propulsion

The most common way to propel a submarine is to use a propeller propulsion system. The propeller will with a rotating motion of its aerofoil shaped blades create a high pressure behind the vessel and at the same time a low pressure at the stern generating a forward motion. In this section, the equations and variables are based on Garme (2007)[7]. The propeller will propel the vessel with a thrust T, calculated as

T = KTρD4n2. (4.34)

With that thrust follows a torque Q as

Q = KQρD5n2. (4.35)

In Equation 4.34 and 4.35 D is the propeller diameter, n is the propeller revolution per second and KTand

KQ are the thrust and torque coefficients, respectively. The thrust and torque coefficients are taken from

(27)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0

0.5 1

Advanc e r at i o J

P rop eller co eﬀcients

K T & η0 KQ KT η0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.90 0.05 0.1 K Q

Figure 4.6: Thrust and torque coefficients as a function of J.

The advance ratio is a dimensionless ratio between the distance the propeller moves forward during one revolution and the diameter of the propeller:

J = u(1 − wf)

nD , (4.36)

where wf is the wake fraction, taking into consideration that the relative fluid velocity will decrease due

to the wake created from the vessel body. The low pressure created at the stern of the submarine by the

propeller gives an increases resistance. This effect is called the thrust deduction factor, tdf

tdf =

∆R

T , (4.37)

where ∆R is the added resistance. The relation between the drag resistance R and the propellers thrust

force can be described as T = R + T tdf. Assuming that the propeller is mounted in the body frame x-axis

and the force is parallel to x, the force in the x-direction is given by

XP= T (1 − tdf). (4.38)

Following the propellers mount assumption the torque will only give a pure contribution in the roll moment equation i.e.

KP= −Q. (4.39)

4.6 Control surfaces

(28)

(a) +-configuration. (b) ×-configuration.

Figure 4.7: The cruci- and cross-form stern hydroplane configuration, where α is the hydroplane mounting angle from the submarines vertical line.

The main difference between the two configurations is that for the ×-configuration, all four hydroplanes are used for every movement whether it is horizontal or vertical. Because of this, there is a built in redundancy, meaning that if one or even two hydroplanes breaks down, the vessel is still manoeuvrable in contrast to the +-configuration. The fact that all four hydroplanes are used for every movement indicates that the hydroplanes does not have to be as big to give the same amount of force to create the rotating moment, and thus creating less drag. The ×-configuration allows the vessel to come closer to the sea bed than a +-configuration without damaging the hydroplanes.

The hydroplanes work as a wing, creating two force components: Drag D, parallel to the incoming flow, and lift L, perpendicular to the incoming flow. The lift and drag can be calculated with the following equations[1]

D =1 2ρu 2_AC D, (4.40) L = 1 2ρu 2_AC L. (4.41)

In these equations A is considered the surface area or the hydroplane. The drag coefficient CD is calculated

as the sum of the parasitic drag CD0 and the lift induced drag CDi as

CD= CD0+ C2 L πaee | {z } CDi . (4.42)

Here, e is the Oswald efficiency factor, which depends on the lift distribution and is in the range of 0.9-0.95

for low aspect-ratio hydroplanes. ae is the aspect ratio of the hydroplane and CL is the lift coefficient. The

parasitic drag is often very small compared to the induced drag and can be neglected. The lift coefficient depend on the angle of attack α as

CL= CLαα, (4.43)

(29)

CLα=

1.8πae

1.8 + cos(Ω) a2e

cos4_(Ω)+ 4

1/2, (4.44)

where Ω is the hydroplanes the quarter chord sweep angle. With these equations the hydroplane forces and moments can be calculated by knowing their placements with respect to the rotation centre. As an example, the force in the x-direction due to the bow-plane can be deduced from Equation 4.40 to 4.43 as

Xδbδb= ρ 2 (CLδbδb) 2 πaee u2A. (4.45)

4.7 Equations of motion

Newton’s second law of motion is described in Equations 4.9 - 4.10, which are reintroduced here for reference:

f0= m( ˙v0+ ω × v0+ ˙ω × rG+ ω×(ω × rG))),

m0= I0ω + ω×(I˙ 0ω) + mrG× ˙v0+ ω × v0,

where f0 and m0 are the external forces and moments, respectively. The right hand side of the equation

is the rigid body dynamics explained in section 4.1. The external forces and moments are described as the sums X X = XHD+ XHS+ XP+ XC, X Y = YHD+ YHS+ YCF+ YV+ YC, X Z = ZHD+ ZHS+ ZCF+ ZV+ YC, X K = KHD+ KHS+ KV+ KP+ KC, X M = MHD+ MHS+ MCF+ MV+ MC, X N = NHD+ NHS+ NCF+ NV+ NC, (4.46)

where the index HD refers to the hydrodynamic forces/moments as described in section 4.2, HS refers to hydrostatics as described in section 4.2.3, CF refers to the non-linear cross flow drag as described in section 4.3, P refers to propulsion as described in section 4.5, V refers to the horse-shoe vortex as described in section 4.4 and C refers to the control surfaces as described in section 4.6. Putting the external forces and moments together with the rigid body dynamics would then give the six EOM’s. The following six equations are based

on the equations given by Feldman[4] and with some updates given by Lennart Boss´er, senior scientist, FOI.

These equations are for a submarine vessel moving with a forward speed. In the following paragraph all coefficients are introduced in dimensionless form, indicated with a prime index. In order to give the correct

(30)

(31)

(32)

(33)

Roll moment equation, K

Ixxp + (I˙ zz− Iyy)qr − ( ˙r + pq)Izx+ (r2− q2)Iyz+ (pr − ˙q)Ixy m[yG( ˙w − uq + vp) − zG( ˙v − wp + ur)] = +ρ 2l 5 [Kp0˙p + K˙ ˙r0˙r + Kqr0 qr + Kp|p|0 p|p| + K 0 r|r|r|r|] +ρ 2l 4_[K0 pup + Kr0ur + Kv0˙˙v + Kwp0 wp] +ρ 2l 3_[K0 ∗u2+ KvR0 uv + K 0 iuvFW(t − τT) + Kv|v|0 v|v|] +ρ 2l 3 [Kδ0_r1u 2 δr1− K 0 δ_r2u 2 δr2− K 0 δ_r3u 2 δr3+ K 0 δ_r4u 2 δr4] +ρ 2l 3_[K0 δ_r1|δ_r1|u 2_δ r1|δr1| − K 0 δ_r2|δ_r2|u 2_δ r2|δr2| − K 0 δ_r3|δ_r3|u 2_δ r3|δr3| + K 0 δ_r4|δ_r4|u 2_δ r4|δr4|] +ρ 2l 2_z 1Cl Z x1 x2 w(x)¯vFW(t − τ [x])dx

+ (yGW − yBB) cos(θ) cos(φ) − (zGW − zBB) cos(θ) sin(φ)

(34)

Pitch moment equation, M

Iyyq + (I˙ xx− Izz)rp − ( ˙p + qr)Ixy+ (p2− r2)Izx+ (qp − ˙r)Iyz m[zG( ˙u − vr + wq) − xG( ˙w − uq + vp)] = +ρ 2l 5 [Mq0˙q + M˙ rp0 rp + Mq|q|0 q|q| + (M 0 rr+ Mrrξ0 ξ)r 2 ] +ρ 2l 4_[M0 ˙ ww + M˙ q0uq + Mvr0 vr + Mwq0 wq + Mw|q|0 w|q| + M|q|δ0 su|q|δs] +ρ 2l 3_[M0 ∗+ Mξ0ξ + M 0 ξθξθ + M 0 ξ|θ|ξ|θ|u 2_{+ M}0 wuw + (M 0 vv+ M 0 vvξξ)v 2_] +ρ 2l 3 [Mvw0 vw + M|w|0 u|w| + Mww0 |w(v 2 + w2)1/2|] +ρ 2l 3_[M0 δsu 2_δ s+ Mδ0s|δs|u 2_δ s|δs| + Mδ0bu 2_δ b+ Mδ0b|δb|u 2_δ b|δb| + Mδ0sηu 2_δ s(η − 1 C)C] +ρ 2Cd Z l xb(x)w(x)[w(x)2+ v(x)2]1/2dx −ρ 2lCl Z x1 x2 xv(x)¯vFW(t − τ [x])dx

− (xGW − xBB) cos(θ) cos(φ) − (zGW − zBB) sin(θ)

(35)

Yaw moment equation, N

Izz˙v + (Iyy− Ixx)pq − ( ˙q + rp)Iyz+ (q2− p2)Ixy+ (rp − ˙p)Izx m[xG( ˙v − wp + ur) − yG( ˙u − vr + wq)] = +ρ 2l 5_[N0 ˙r˙r + N 0 r|r|r|r| + N 0 ˙ pp + N˙ 0 pqpq] +ρ 2l 4_[N0 pup + (Nr0+ Nrξ0 ξ)ur + Nv0˙˙v + Nv|r|0 v|r| + N 0 |v|r|v|r + Nwr0 wr + N|r|δ0 ru|r|δr] +ρ 2l 3_[N0 ∗u2+ (Nv0 + Nvξ0 ξ)uv + Nv|v|R0 v|(v 2_{+ w}2₎1/2_|] +ρ 2l 3_[N0 δru 2_δ r+ Nδ0r|δr|u 2_δ r|δr| + Nδ0rηu 2_δ r(η − 1 C)C] −ρ 2Cd Z l xh(x)v(x)[w(x)2+ v(x)2]1/2dx −ρ 2lCl Z x1 x2 xw(x)¯vFW(t − τ [x])dx

(36)

Chapter 5

Simulation method

The mathematical model for a submarine is a system of six differential equations, Equations 4.47 - 4.52, that describes the motion in six degrees of freedom. The different coefficients and dimensionless derivatives in the EOM’s are given in a submarine coefficient text file from a program provided by FOI. The following chapter describes the method of how the equations are used to extract the velocity and position of a submarine.

5.1 Assumptions

The following assumptions are used for this method: • The submarine is modelled as a single rigid body • The submarine is fully submerged.

• The mass and mass distribution are constant.

• The submarine is modelled in a calm environment i.e. this model is neglecting wave interactions. • Rotational effects of the earth are negligible, making the transformation between earth-fixed and

body-fixed reference frame valid.

• Fluid properties are considered constant.

• Local variances of the gravitational field are ignored.

• One propeller propulsion system mounted in the body frame x-axis only contributing with XP =

T (1 − tdf) and KP= −Q

• The submarine is moving with enough velocity and is large enough to yield Reynolds number over 5.5 · 105.

• The drag coefficient Cdis considered constant and not dependent on orientation

• The lift coefficient Cl is considered constant and not dependant on orientation

(37)

5.2 Non-linear cross flow drag

The non-linear cross flow is, as described in section 4.3, calculated using the strip method where the total force is computed by integrating each two dimensional section over the length of the submarine. The drag force in the y-direction is

YCF= − ρ 2Cd Z l h(x)v(x)pw(x)2_{+ v(x)}2 | {z } f(x) dx.

The integral is solved by using the numerical implementation of the trapezoidal rule[3]:

I = h 2 N X k=1 −1 (f (xk+1) + f (xk)) , h = L N − 1. (5.1)

The widths, heights and longitudinal locations x of each two dimensional segment are given in the coefficient file. However, the longitudinal location is given with its origin at the aft of the ship. The longitudinal

locations must be transformed to the body reference system by subtracting xap, see in Figure 5.1

Figure 5.1: Location of the aft longitudinal position with respect to centre of rotation, xap.

5.3 Horse-shoe vortex

The horse-shoe vortex is, as described in section 4.4, calculated by integrating the lift over a span between

the vortex starting point x1to the position x2where the vortex leaves the body. If the drift angle β is small,

the vortex will not leave the body until it reaches the aft of the vessel. However, as the drift angle increases, detachment from the hull occurs earlier. The location of the detachment from the hull is calculated according to Feldman[4] and Bohlmann[2] as

x2=

(

xap for|β| ≤ βST

x1− (x1− xap)(S1+ S2|β|) for|β| > βST,

(5.2) where the drift angle is calculated as

β = −tan−1v

u

, (5.3)

where S1, S2 and βST are constants that are obtained from experimental measurements. βST is the drift

angle where the detachment of the vortex begins to travel from the aft towards the tower. As explained

in section 4.4, the local lateral velocity ¯vFW at the sail generates a vortex, which travels down along the

(38)

velocity at the sail some time offset τ back in time. This offset can be calculated with Equation 4.31. Since the velocity u is constant over the span of the hull the time offset can be expressed as

τ =x1− x

u(t) , (5.4)

where x is a vector of locations between x1 and x2, giving a range of time offsets for each location. The

local lateral velocity for a location x is then calculated as stated in Equation 4.33, which is dependent on the lateral velocity component v and the rotational velocities r and p. However, the velocity components are needed at specific times t−τ , and the integrator is solving with a time-step that is either constant or adaptive depending on the preferred algorithm and error limit. A one-dimensional spline is used to interpolate the previous time-step solutions to render the lateral velocity as a continuous function of time. When the correct

lateral velocity is interpolated for all the locations between x1and the detachment at x2, the integral is solved

using the numerical trapezoidal method in Equation 5.1.

5.4 Control surfaces

In the coefficient file, the stern hydroplanes are calculated assuming a cross-configuration (×−conf iguration), and each coefficient is considered for one hydroplane. Therefore, the coefficients regarding the stern hy-droplanes must be multiplied with 4 before using in Equations 4.47 - 4.52. As presented in section 4.6, when using the cross-configuration all four stern hydroplanes are used to perform a manoeuvre, but for a cruci-configuration only two are used as shown in Figure 5.2, where the submarine is initiating a horizontal turn. Positive hydroplane deflection is with its trailing edge towards port side as the figure is portraying.

(a) +-configuration (cruci) seen from behind. (b) ×-configuration (cross) seen from behind.

Figure 5.2: The cruci- and cross-form hydroplane forces due to their deflection Fδand their force components

(39)

Assuming that the cross-configuration has a mounting of 45 degrees from the submarines vertical line, the hydroplane force for one blade for each configuration will be

Fδy×=

Fδ

√ 2,

Fδy+= Fδ. (5.5)

One blade on the cruci-configuration will thus create a force of√2 times larger than a blade with a sweep of

45 degrees. If a cruci-configuration model submarine is used, each of the stern hydroplane coefficients should

be multiplied with a factor of 2√2 since the cruci-configuration only uses two hydroplanes instead of four

for a turn or a dive.

In the EOM’s, Equations 4.47 - 4.52, the stern hydroplanes are treated with two variables, except in the rolling moment equation, Equation 4.50, where each plane has its own variable. The horizontal manoeuvring

variable is the rudder deflection δrcalculated as

δr×= (δ1+ δ2+ δ3+ δ4) 4 , (5.6) δr+= (δ2+ δ4) 2 , (5.7)

where δ1− δ4are the angles of the four different hydroplanes as seen in Figure 5.2. The vertical manoeuvring

variable is the stern-plane deflection δs calculated as

δs×= (−δ1+ δ2− δ3+ δ4) 4 , (5.8) δs+= (−δ1− δ3) 2 . (5.9)

For this thesis, the deflection of the stern hydroplanes are defined as positive when the trailing edge is to-wards port side of the submarine and negative toto-wards starboard.

On a submarine there are some dynamic movements that has to be considered to be able to get as good a model as possible. Such dynamics is considered with the hydroplanes. As the controller orders a deflection of the hydroplanes to execute a turn or dive, the planes will be rotating with a limited rate until the desired deflection is reached. This rotation dynamic is described with a first order differential equation:

˙δi=

(δc− δi)

Tr

, i = 1, ..., 5, (5.10)

where i = 1 − 4 are the stern rudders and i = 5 is for the bow plane. δc is the commanded deflection and

δi is the current deflection. Tr is a time constant for the change of the hydroplanes. To control that the

rotational velocity is not changing to fast and that the deflection does not overshoot, they are controlled with the following two statements:

| ˙δ| ≤ | ˙δ|max,

|δ| ≤ δmax. (5.11)

5.5 Propulsion

(40)

J = u(1 − wf)

nD ,

where wf and D are given in the coefficient file. When the advance ratio is known the thrust coefficient is

calculated as1

KT= kT0+ kTJJ + kTJ2J2+ k_TJ3J3+ k_TJ4J4. (5.12)

The torque coefficient is calculated with the same method as

KQ= kQ0+ kQJJ + kQJ2J2+ k_QJ3J3+ k_QJ4J4. (5.13)

Here, kQ0, kQJi, kT0 and kTJi are polynomial multipliers from the input file. With Equation 4.34 and 4.35

the thrust and torque are calculated, respectively, as

T = ρKTn2D4,

Q = ρKQn2D5.

Just like the hydroplanes cannot move instantly to the commanded deflection, the engine needs some time to change to a new desired rotational speed. Therefore, the rotational speed rate is controlled by a first order differential equation

˙n = (nc− n)

Tn

, (5.14)

where nc and n are the commanded and current revolutions per second, respectively. Tnis a time constant

for the change of the rotational speed. The rotational speed and its time derivative are limited according the following two statements:

| ˙n| ≤ | ˙n|max,

|n| ≤ nmax. (5.15)

The propeller moment must not exceed the maximum allowed propeller moment. If, during an acceleration, the propeller moment exceeds the maximum allowed moment, the rate of change of the propeller rotational speed is reduced by setting the time derivative to

˙n = Qmax

Q − 1

n. (5.16)

If the propeller is going in reverse the torque will be negative instead and therefore the moment must not

transcend the negative maximum allowed propeller moment, −Qmax. If this occurs, the rate of change of

the propeller rotational speed is limited to

˙n = −Qmax

Q − 1

n. (5.17)

5.6 Model

The mathematical model can, as described before, be expressed with Newton’s second law of motion

A ˙ν = F , (5.18)

(41)

      m−ρ₂l3X_u0_˙ 0 0 0 zGm −yGm 0 m−ρ₂l3_Y0 ˙ v 0 −yGm−ρ2l 4_Y0 ˙ p 0 xGm−ρ2l 4_Y0 ˙ r 0 0 m−ρ₂l3_Z0 ˙ w yGm −xGm−ρ2l 4_Z0 ˙ q 0 0 −zGm−ρ2l 4_K0 ˙ v yGm Ixx−ρ2l 5_K0 ˙ p −Izy −Izx−ρ2l 5_K0 ˙ r zGm 0 −xGm−ρ2l 4_M0 ˙ w −Ixy Iyy−ρ2l 5_M0 ˙ q −Iyz −yGm xGm−ρ2l 4_N0 ˙ v 0 −Izx−ρ2l 5_N0 ˙ p −Iyz Izz−ρ2l 5_N0 ˙ r               ˙ u ˙v ˙ w ˙ p ˙ q ˙r         =         X Y Z K M N         . (5.19)

At this form the equation can be re-written to be solved for the linear and angular accelerations, i.e. ˙

ν = A−1F . (5.20)

This form will thus give six ordinary differential equations that can be written as

˙ u = A−1₁₁X + A−1₁₂Y + A−1₁₃Z + A−1₁₄K + A−1₁₅M + A−1₁₆N, ˙v = A−1₂₁X + A−1₂₂Y + A−1₂₃Z + A−1₂₄K + A−1₂₅M + A−1₂₆N, ˙ w = A−1₃₁X + A−1₃₂Y + A−1₃₃Z + A−1₃₄K + A−1₃₅M + A−1₃₆N, ˙ p = A−1₄₁X + A−1₄₂Y + A−1₄₃Z + A−1₄₄K + A−1₄₅M + A−1₄₆N, ˙ q = A−1₅₁X + A−1₅₂Y + A−1₅₃Z + A−1₅₄K + A−1₅₅M + A−1₅₆N, ˙r = A−1₆₁X + A−1₆₂Y + A−1₆₃Z + A−1₆₄K + A−1₆₅M + A−1₆₆N. (5.21)

To get the translation and rotational motions, the ordinary differential equations in Equation 5.21 are solved by numerical integration from an initial state vector, using the Adams’ method[3]. The Adams’ method is based on both forward and backward Euler method, using a multi-step approach with Adams-Bashforth method as a predictor and Adams-Moulton as a corrector. The Adams’ method is described in Appendix B. With the same method, the earth-fixed coordinates and Euler-angles are extracted. From Equation 4.2 the following differential equations are extracted:

˙

φ = p + q sin(φ) tan(θ) + r cos(φ) tan(θ), ˙

θ = q cos(φ) − r sin(φ), ˙

ψ = [q sin(φ) + r cos(φ)] / cos(θ), ˙

x0= u cos(θ) cos(ψ) + v(sin(φ) sin(θ) cos(ψ))

+ w(sin(φ) sin(ψ) + cos(φ) sin(θ) cos(ψ)), ˙

y0= u cos(θ) sin(ψ) + v(cos(φ) cos(ψ) + sin(φ) sin(θ) sin(ψ))

+ w(cos(φ) sin(θ) sin(ψ) − sin(φ) cos(ψ)), ˙

w0= −u sin(θ) + v cos(θ) sin(φ) + w cos(θ) cos(φ). (5.22)

In this model there are eighteen derivatives and therefore the same number of state variables. Meaning that there are eighteen variables that change with time. So, after integrating the differential equations for a certain time interval, the following state vector will be updated at each time step:

x =

(42)

5.7 Post Processing

With the state variables calculated for each time step, more components can be calculated for presenting, analysing and verifying purposes. The following section is covering all the post processing variables that are calculated. Many of the variables and equations in this section has been shown earlier in the report, but are mentioned again for the readers convenience.

From the translational velocity components u, v and w, the translation speed V is calculated as

V =pu2_{+ v}2_{+ w}2_. _(5.23)

The composite stern hydroplane angles δr (turn) and δs (dive) for an ×-configuration are calculated from

individual hydroplane angles with Equations 5.6 and 5.8 as

δr= (δ1+ δ2+ δ3+ δ4) 4 , δs= (−δ1+ δ2− δ3+ δ4) 4 .

The drift angle β is a key component during turns and has a significant influence of the vessels turning ability. The drift angle is calculated with Equation 5.3 as

β = −tan−1v

u

.

From post processing, it is of importance that the propulsion variables can be analysed. Therefore, the advance ratio J, the thrust T and torque Q are calculated with Equations 4.34, 4.35 and 4.36 as

J = u(1 − wf)

nD ,

T = ρKTn2D4,

Q = ρKQn2D5.

Adding these variables and the time steps to the integrated state variables will give a matrix with the following 26 components to be used for analysis of the submarine:

x =

t u v w p q r φ θ ψ x0 y0 z0 δ1 δ2 δ3 δ4 δb n V δr δs β J T Q

(43)

5.8 Tests

With the model and post processing described in the previous sections, different tests can be simulated to verify if the submarine design is fulfilling the requirements or if the design has to be altered. The following tests are currently implemented in the program:

• Acceleration test. • Stern plane deflection. • Turning circle test. • Meander test. • Straight line test. • Hydroplane failure. • Control-system failure. • Zig-zag test.

• Step response test.

(44)

Chapter 6

Validation and verification

Validation is the process of assuring that a product or system meets the stakeholders’ needs. In this case the complete method, with towing-tank experiments with scale models to measure the hydrodynamic deriva-tives and applying these in a system of linearized differential equations, should be compared to the dynamic behaviour of a full-scale submarine. Such validation is beyond the scope of this thesis. The fact that the tank-test and simulation approach has been used for decades as a standard method for analysing of both surface and submarine ships is real-life evidence that the method is useful for its purpose.

Verification of the results from the tests listed in section 5.8 is a very important step, due to the fact that these are the results that the design of the submarine is verified on. The best way to verify the model is by taking full-scale submarine test data and compare it to the simulated results. However, due to the fact that submarine designs are sensitive information as are their manoeuvring capabilities, this verification process cannot be done. The new software is therefore compared to FOI’s current simulation software. This software creates an output file containing a matrix with the following vectors:

x =

t x0 y0 z0 ψ δ1 δ2 δ3 δ4 u v φ δs δr V −ψ −δs β θ r δb Q n J

This means that every result from the current simulation software can be directly compared to results simulated in this thesis, see section 5.7. However, due to the fact that the programs are using different integration methods, the time steps will differ and a visual comparison of the time solutions is deemed sufficient. The verification can be done by comparing all variables for each and every test. However, this is a bit time consuming and not very effective. During development of the software a structured approach to verification was needed. This approach, described below, focused on a subset of the available test cases, namely:

• Straight line test. • Bow-plane step response. • Acceleration test. • Turning Circle test. • Stern-plane deflection test.

(45)

The hydroplane force from the bow plane creates the least moment on the hull since its location is close to the CG, restricting the coupling to be mostly between the X - and Z -equations. The bow-plane step response also gives a chance to see that the hydroplane dynamics works. The acceleration test has the same coupling as a straight line test and gives an opportunity for further verification of the propulsion dynamics. After these tests, the amount of coupling is increased by starting to use the stern hydroplanes. The turning circle test will increase the coupling mainly between the X -, Y - and N -equations, while the stern-plane deflection is exercising the coupling between the Z - and M -equations.

Special concern is used when isolating the effects of the non-linear cross flow drag and the horse-shoe vortex.

This is done by manually manipulating the coefficient file and setting CL= CD= 0. The verification is thus

done in two steps for the above described tests. The two steps being:

1. Without contribution from non-linear cross flow drag and horse-shoe vortex. 2. With contribution from non-linear cross flow drag and horse-shoe vortex.

With this two-step approach it is possible to first check the linear part of the equations to verify a correct implementation and then proceed to verification of the non-linear contributions. As mentioned in sections 4.3 and 4.4, the non-linear cross flow drag and horse-shoe vortex are calculated with strip theory, which integrates two dimensional strips along the hull to get the three dimensional effect. For the non-linear cross flow drag, the number of two dimensional sections is given in the coefficient file. For the horse-shoe vortex the number of sections for integration is set by the designer of the program. If the data points are too few, the trapezoidal method might over/under estimate the force on the hull. However, too many points will increase the computational time. Therefore, a verification of the necessary number of locations along the hull to be integrated is performed. Increasing the amount of locations until the change is low enough in the result yields the accurate number of data points needed. The effect of the horse-shoe vortex is determined to be with 20 data points spaced with equal distance along the hull. The determination of the number of data points is shown in Appendix C.

The following sections presents the five chosen tests and their two-step comparison between the program created for this thesis and the program currently used by FOI. The current program is referenced to as ”Current”, while the program written for this thesis is referenced as ”New”. Each test is simulated with

the Flipper Topolino submarine1_{. The submarine general arrangement is shown in Appendix D and its}

derivative, mass properties and miscellaneous inputs are listed in Appendix E. The coded program steps and implementations are explained in the program structure in Appendix F.

6.1 Straight line test

The aim of the test is to see how the vessel is behaving without any hydroplane induced forces. The

per-formance of the submarine is measured by its loss/gain in depth z0, its pitch θ and yaw ψ angles. The

simulation begins with the vessel moving at a constant velocity of u0= ucfor a user defined simulation time.

The test is simulated with the following input variables:

uc= 5 knots

tsim= 200 s

The simulated result is shown in Figure 6.1.

1_{The Flipper Topolino is a submarine concept developed as a student task during a course in submarine design at FOI}

(46)

0 20 40 60 80 100 120 140 160 180 200 −0.6 −0.55 −0.5 −0.45 −0.4 −0.35 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 Ya w & P it c h [ ◦] t ime [s ] 0 20 40 60 80 100 120 140 160 180 200 197.2 197.4 197.6 197.8 198 198.2 198.4 198.6 198.8 199 199.2 199.4 199.6 199.8 200 De p th [m ] Yaw Pit ch De pt h

Figure 6.1: Yaw, pitch and depth as a function of time resulting from constant speed and no hydroplane deflection.

The figure shows how the vessel behaves without any hydroplane activities. Initially the pitch is decreasing

to a stable value of -0.45◦. This is a rather small angle, but as the depth plot shows, it has the effect of

decreasing the depth with almost 3 m. This shows that the vessel has a natural pitch moment as it moves with a forward velocity, which is a result of the drag from the sail and bow-plane hydroplanes. The vessel also has a tendency to ”pull” towards one side as seen from the yaw angle that is increasing with a almost constant slope. This means that the vessel will perform a slow turn in the horizontal plane as well. This is probably a result of the propeller torque giving the vessel a slight roll angle combined with the effects of drag on the sail.

Step 1.

The linear comparison between the two software’s using the same inputs as in the previous simulation are shown in Figure 6.2. 0 20 40 60 80 100 120 140 160 180 200 9.227 9.228 9.229 9.23 9.231 t ime [s ] To rq u e , Q [kN m ] 0 20 40 60 80 100 120 140 160 180 200 0.6355 0.6356 0.6356 0.6357 0.6357 t ime [s ] Ad v an c e rat io, J [-] Cur r ent New (a) Propulsion. 0 50 100 150 200 250 300 350 400 450 500 −0.4 −0.3 −0.2 −0.1 0 0.1 x0[m] y0 [m ] Cur r ent New 0 50 100 150 200 250 300 350 400 450 500 197 197.5 198 198.5 199 199.5 200 x0[m] z0 [m ] (b) Path.

(47)

From the comparison in the figures it can be seen that the forward motion correlates very well with the current program. There is a small difference in the torque that is constant over the whole simulation. This might be due to different precision and libraries for floating-point calculations. For an example, the current program and the new have small difference in the initial velocity. The EOM’s are handling its variables according the

metric system so as the speed is transformed from knots to m/s the rounding gives: u0,Current= 2.5720 m/s

while u0,New= 2.5722 m/s. There is also a small uncorrelated phenomena in the vertical plane at the end of

the simulation. This can be due to floating-point arithmetic’s, but can also indicate that there is something wrong with the Z - or M - equations.

Step 2.

Adding the non-linear terms for the simulation gives the results shown in Figure 6.3.

0 20 40 60 80 100 120 140 160 180 200 9.227 9.228 9.229 9.23 9.231 t ime [s ] To rq u e , Q [kN m ] 0 20 40 60 80 100 120 140 160 180 200 0.6355 0.6356 0.6356 0.6357 0.6357 t ime [s ] Ad v an c e rat io, J [-] Cur r ent New (a) Propulsion. 0 50 100 150 200 250 300 350 400 450 500 −0.4 −0.3 −0.2 −0.1 0 0.1 x0[m] y0 [m ] Cur r ent New 0 50 100 150 200 250 300 350 400 450 500 197 197.5 198 198.5 199 199.5 200 x0[m] z0 [m ] (b) Path.

Figure 6.3: Straight line test comparison between Current (blue) and New (red) simulation programs with non-linear terms. The left two graphs present a comparison in the propulsion terms J and Q. The right two graphs show a comparison of the vessel path both in the horizontal (top) and vertical plane (bottom).

From the comparison in the figures it can be seen that the forward motion correlates very well with the current program. The same difference that was discussed in the previous test, see Figure 6.2, is seen in this test as well. This is expected since the effect from the non-linear terms are insignificant due to the low lateral velocities and drift angle in a straight line test. The error in the vertical path is however smaller with the non-linear terms as shown in Figure 6.3(b).

6.2 Bow-plane step response

The step-response is a test where the rudder is set at a certain deflection at certain times. Assuming that the vessel has a pure velocity in the surge direction and that the initial value of the velocity is equal to the

commanded velocity, u0 = uc. All the other initial values are set to zero. The first hydroplane execution

δc = δ1 at t = 0 is held until the next time step, t = t1, where the second execution takes place, δc = δ2.

(48)

Figure 6.4: Example of a step-response test with rudder angle as a function of time.

The bow-plane step response is simulated with the following inputs:

δb0= 10◦ at t0= 0 s δb1= 5◦ at t1= 100 s δb2= −5◦ at t2= 200 s δb3= −10◦ at t3= 300 s u0= uc= 5 knots tsim= 400 s

Results form the simulation is shown in Figure 6.5.

0 50 100 150 200 250 300 350 400 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 1 2 3 4 5 6 7 8 9 10 11 δb &P it c h [ ◦] t ime [s ] 0 50 100 150 200 250 300 350 400 183 185 187 189 191 193 195 197 199 201 De p th [m ] δb Pit ch De pt h

Figure 6.5: Bow-plane deflection, pitch and depth as a function of time during a bow-plane step response test.

(49)

Step 1.

Using the same inputs as the previous simulation, a comparison of the two software linear terms is done and illustrated in Figure 6.6. 0 50 100 150 200 250 300 350 400 −15 −10 −5 0 5 10 15 t ime [s ] δb [ ◦] Cur r ent New

(a) Bow-plane deflection.

0 100 200 300 400 500 600 700 800 900 1000 −2 −1.5 −1 −0.5 0 0.5 x0[m] y0 [m ] Cur r ent New 0 100 200 300 400 500 600 700 800 900 1000 180 185 190 195 200 x0[m] z0 [m ] (b) Path comparison.

Figure 6.6: Bow-plane step response comparison between Current (blue) and New (red) simulation programs. The left graph show the bow-plane deflection over the simulation. The right two graphs show the horizontal (top) and vertical (bottom) paths of the submarine.

(50)

Step 2.

Adding the non-linear terms gives the result shown in Figure 6.7.

0 50 100 150 200 250 300 350 400 −15 −10 −5 0 5 10 15 t ime [s ] δb [ ◦] Cur r ent New

(a) Bow-plane deflection

0 100 200 300 400 500 600 700 800 900 1000 −2 −1.5 −1 −0.5 0 0.5 x0[m] y0 [m ] Cur r ent New 0 100 200 300 400 500 600 700 800 900 1000 180 185 190 195 200 x0[m] z0 [m ] (b) Path comparison

Figure 6.7: Bow-plane step response comparison between Current (blue) and New (red) simulation programs with the non-linear terms. The left graph show the bow-plane deflection over the simulation. The right two graphs show the horizontal (top) and vertical (bottom) paths of the submarine.

From the figure it is visible that the calculations correlates well, both in the horizontal and vertical plane. The small error in the vertical plane, see Figure 6.6(b), is reduced as the non-linear terms are introduced, see Figure 6.7(b).

6.3 Acceleration test

This test illustrates how the vessel moves during acceleration. It can be seen how different propeller dynamics can affect the time it takes for the vessel to reach a constant velocity. The vessel starts with initial velocity

u0= v0 = w0 = 0 and a commanded propeller revolution speed nc and will be simulated for a simulation

time tsim. To control the propellers dynamics, the maximum torque Qmax and rotational acceleration ˙nmax

is set by the user. The resulting revolution velocity n and the submarine’s velocity V are then plotted as functions of time.

The result from the test shown in Figure 6.8. was simulated with the following input variables:

(51)

0 20 40 60 80 100 120 140 160 180 200 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 t im e [s ] Re v o lu ti o n s [r p m ] 0 20 40 60 80 100 120 140 160 180 2000 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5 Ve lo c it y [k n o ts ] R e volut ion Ve loc ity

Figure 6.8: Revolution speed and submarine velocity as a function of time from an acceleration test simulated for 200 s.

In the figure the revolution speed of the propeller and the submarines velocity is plotted as a function of time. The submarine velocity is increasing towards a constant velocity, 9.7-9.8 knots. The revolution speed has a steep slope during the first few seconds of the simulation following the maximum revolution rate limit. From 5 to 80 seconds the speed of the revolution is limited by the maximum allowed torque of the propulsion, which is achieved at 85 rpm.

Step 1.

The linear terms are compared to the current software simulating with the same inputs as previous. Figure 6.9 presents the results of the comparison.

0 20 40 60 80 100 120 140 160 180 200 0 1 2 3 4 5 6 t ime [s ] Sur g e ve l. u [kn ] Cur r ent New 0 20 40 60 80 100 120 140 160 180 200 0 20 40 60 80 100 120 t ime [s ] Re v o lu ti o n v e l. n [r p m ]

(a) Velocity comparison

0 20 40 60 80 100 120 140 160 180 200 0 0.2 0.4 0.6 0.8 t ime [s ] Ad v an c e rat io J [-] Cur r ent New 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 t ime [s ] To rq u e , Q [kN m ] (b) Propulsion comparison

Figure 6.9: Acceleration test comparison between Current (blue) and New (red) simulation programs. The left two graphs show the surge velocity (top) and the propeller rotational velocity (bottom). The right two graphs show the propulsion terms J (top) and Q (bottom).

(52)

Step 2.

Adding the non-linear terms results in the comparison shown in Figure 6.10.

0 20 40 60 80 100 120 140 160 180 200 0 1 2 3 4 5 6 t ime [s ] Sur g e ve l. u [kn ] Cur r ent New 0 20 40 60 80 100 120 140 160 180 200 0 20 40 60 80 100 120 t ime [s ] Re v o lu ti o n v e l. n [r p m ]

(a) Velocity comparison.

0 20 40 60 80 100 120 140 160 180 200 0 0.2 0.4 0.6 0.8 t ime [s ] Ad v an c e rat io J [-] Cur r ent New 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 t ime [s ] To rq u e , Q [kN m ] (b) Propulsion comparison.

Figure 6.10: Acceleration test comparison between Current (blue) and New (red) simulation programs with the non-linear terms. The left two graphs show the surge velocity (top) and the propeller rotational velocity (bottom). The right two graphs show the propulsion terms J (top) and Q (bottom).

The non-linear terms have no significant influence on the acceleration test and the good correlation from Figure 6.9 is maintained.

6.4 Turning circle

The turning circle test is a measure of the ability to turn the vessel. The test is simulated for a certain

commanded rudder deflection δrc and commanded velocity uc. At the start of the turn, the submarine is

moving with constant forward speed u0= uc. All the other initial conditions are set to zero. The test results

are then plotted in a graph in the earth-fixed horizontal reference frame, meaning that the coordinates for

the x0(t) and y0(t) are plotted for each time step to illustrate the turning circle of the submarine, see Figure

Simulation of Submarine Manoeuvring Sebastian Thun´e A thesis presented for the degree of Master of Science