Drift Forces and Moments
TECHNISCHE UNIVERSITEIT Scheepshydramechanica Archief Mekelweg 2, 2628 CD Delft Tel:015-786873/Fax:7818369 In ,plaats van de nadruk te leggen op de effecten van het westerse ko-lonialisme op de europese volkeren, zou bij het geschiedenis onderwijs op middelbare scholen veel meer aandacht geschonken dienen te worden
aan de onthutsende gevolgen die dit kolonialisme heeft op de overheerste volkeren
10. Het verplicht stellen van boeken van Jan Wolkers voor de literatuurlijst van middelbare scholieren komt voort uit schroom sexuele voorlichting te geven.
U. In besprekingen van het oeuvre van Heinrich Heine wordt het 'Such der Lieder' ten onrechte aangemerkt al's een romantische jeugdzonde.
eelten habt Ihr mich verstanden, eelten auch verstand ich Euch, Mur wenn wit im Kot uns fanden,
ea verstanden wir uns gleich. Buda der Lieder Die Heimkehr LXXVIII
behorend bij het proefschrift
Time-Domain Calculations of Drift Forces
and Moments
van
Henk. J. Prins
11. Het grootste resterende taboe onder tedmische wetenschappers is het geloof. Opvallend genoeg bestaat er onder deze baanbrekende mensen geen enkele behoefte dit taboe teldoorbreken,
2. De correlatie tussen geleerdheid en egocentrisme ligt vervaarlijk dicht biji
1.
Gezien de tegenstrijdige belangen die er spelen bij de verschillende ver-kiezingen, is het uit staatsrechtelijk oogpunt onjuist meer dan een soort verkiezing te houden in een tijdsbestek van een jaw-.
4.. Bij het recentelijk ingevoerde, plaatselijk corrigerend referendum, moet de tegenpartij naast het ophalen van voldoende handtekeningen, ook zorg dragen voor een minimale opkomst en een meerderheid van stemmen. Dit lijkt op het gevecht van David tegen Goliath. Het zou een wonder zijn als de 'David' dit zou winnen.
Rechtvaardiger zou zijn de locale overheid te laten aantonen dat haar besluit op genoeg steun kan rekenen van de bevolking.
5. Het summum van politieke macht is het kunnen beslissen tot het aan-leggen van een ongewenste en onrendabele spoorlijn. Het summum van
politick onbenul is het daadwerkelijk nemen van dat besluit.
'6. Het summum van politick vernuft is het maken van een politieke blunder aan het eind van een ambtstermijn. Dit vergroot de kans op het krijgen van een erebaan.
7. Een rechtstaat staat of valt met handhaving van de openbare orde en de wet.. Bij de handhaving van de wet wordt echter te veel rekening gehouden met de kans op, escalatie, en dus met verstoring van de openbare orde..
of Drift Forces and
Moments
Tijdsdomein Berekeningen van Drift Krachten en
Moinenten
TECHNISCHE UNIVERSITET
laboratorium vow Scheepshydromechanka
Archlet
PROEFSCIIRIFT Mekehveg 2.2828 CD Delft
Tat: 015 - 788873 . Fait 015 . 781838
ter verkrijging van de graad van doctor aan de Technische Universiteit
op gezag van de Rector Magnilicus Prof. ft. K.F. Walker, in het openbaar te verdedigen ten overstaan van een commissie,
door hei college van Dekanen aangewezen, op vitidag 31 ma.art 1995 te 16.00 liar
door
Henk Jirrnien PRINS,
wiskundig
geborenitte Amsterdam.
Delft,
Prof. dr. ir. A.J. Hermans
Samenstelling promotiecommissie:
Prof. dr. O.M. Faltinsen (Norwegian Institute of Technology) Prof. dr. ir. A.W. Heemink (TU Delft)
Prof. dr. J.N. Newman (MIT)
Prof. dr. ir. J.A. Pinkster (TU Delft) Prof. dr. ir. P..). Zandbergen (U Twente) Dr. ir. J.E.W. Wichers (Mann, Wageningen)
Copyright @I995 by H.J. Prins
CIP_DATA KONINKLIJKE BIBLIOTHEEK, DEN HAAG
Prins, Henk Jurriiin
Time-domain calculations of drift forces and moments / Henk Jurrien Prins. - : s.n.1.
Thesis Technische Universiteit Delft. With index, ref. -With summary in Dutch.
ISBN 90-9007986-6
Nail 811
Subject headings: ship-hydrodynamics / time-domain simulations / drift forces.
Introduction'
List of symbols 11
1
Problem Definition and Mathematical Model
131.1 Physical problem 13
1.2 Mathematical model 15
1.3 Linearization of the boundary condition 18
1.3.1 Free-surface COD dition 19
1.3.2 Body boundary condition 20
1.4 Forces and moments 20
1.4.1 Perturbation. series 22
1.4.2 Hydrodynamic coefficients and equation of motion 24
1.4.3 Drift forces and moments 26
1.5 Double-body potential 28
2
Numerical Method
292.1 Introduction 29
2.2 Boundary-integral method 30
2.2.1 Discretization of the boundary conditions 32
2.3 Time integration 32
2.3.1 Algorithm for zero speed 33
2.3.2 New algorithm for forward speeds 37
2.3.3 Suggestion concerning time integration of non-linear
equations 39
2.4 Matrix solver 40
3
Two-dimensional test problem: cylinder of infinite length 43
3.1 Introduction 43
3.2 Simplified model 44
3.3 Numerical aspects 46
3.3.1 Convergence and stability 48
.
...
....v
....
....
. . .3.3.2 Effectiveness of absorbing boundary condition . . . . 49
3.4 Results 50
3.4.1 Comparison of results found by using the
double-body potential and the undisturbed-flow potential . 56
3.5 Conclusions 60
4 Three-dimensional test problem: floating sphere
614.1 Introduction 61
4.2 Mathematical model 62
4.2.1 Absorbing boundary condition 63
4.3 Numerical aspects 65
4.4 Results 67
4.5 Conclusions 73
5
Results for a 200kDWT tanker
755.1 Introduction 75
5.2 Mathematical model 76
5.3 Numerical aspects 77
5.3.1 Numerical differentiation on the hull 78
5.3.2 Double-body potential 80
5.4 Results 84
5.4.1 Hydrodynamic coefficients 86
5.4.2 Improvement of roll damping . .
....
895.4.3 Drift forces and moments 93
5.5 Conclusions and recommendations qg
6
General time signals
1016.1 Introduction 101
6.2 General absorbing-boundary condition 102
6.3 Step-response functions 10.4
6.3.1 Fitting of the step-response function using Lag ecre
polynomials 107
6.3.2 Results 108
6.3.3 Forward speed effects 111
6.4 Slow-drift forces 114
6.4.1 Results for infinite depth 116
6.4.2 Results for finite depth 117
6.5 Conclusions 118 Conclusions 119 . . . . . . . . . . .... . . .
...
...
.A Asymptotic analysis of the Green's function
123A .1 Green's function 123
A.2 Asymptotic analyses of the function tP 125
A.2.1 Complex integration 125
A.2.2 Method of stationary phase 127
A.2.3 First correction 129
B Error estimation of the absorbing-boundary condition
133BA Introduction 133
13.2 Error in the potential 134
13.3 Error in the normal derivative 135
BA Integral over the Green's function 136
B.5 Total error 138
Summary
139Samenvatting
141 Acknowledgements 145Curriculum Vitae
147 Bibliography 149 Index 153 . . .When a ship is sailing on an ocean, it is exposed to several forces. These forces are due to wind, waves, current, and speed of the vessel. Especially the first three may cause the ship to loose its course, i.e. to drift away. Drifting forces can be very large, in particular the force due to the incom-ing waves. In this thesis this part of the driftincom-ing effect is studied, and the effect; of forward speed of the vessel is taken into account.
The drifting of a ship due to waves may be surprising. In the case of har-monic waves, the motion of the ship may be expected to be harhar-monic as
well, thus no drifting would be expected. However, drifting does occurs, even in harmonic waves. This may be illustrated by two practical examples. When kids play football, it is very common, especially in Holland, that the ball ends up in a ditch or pool. A very practical way to get the ball out
of the water, is throwing stones or dirt at the ball. The stones generate
waves, which cause the ball to drift to the shore. This drifting even occurs when the waves are small. A second example is a. ship lying at anchor in waves. If the waves do not come head on, the ship wants to drift away into the direction of wave propagation. Because this is prevented by the anchor, the ship will rotate until it lies in head waves.
The examples show that drifting effects can be large and practical. They also show that the effect of drifting is important for both moored systems and freely floating objects. In the case of moored systems the importance of drifting may seem to lessen, when the system lies in head waves, as the example illustrates. I Iowever, the drift force is not a steady force, but may also include slowly varying components. The frequency of thesecomponents may be close to the eigen-frequency of the mooring system, thus possibly
yielding severe damage to the system. This damage might for instance be breaking of anchor lines or the destruction of the mooring system.
In the case of a freely floating ship a major contribution to the drifting ef-fect is caused by incoming waves. However, from measurements it appears that the forward speed of the ship increases the drift forces considerably. Thus the forward speed has tubetaken into account, together with existing
Figure 1.1: Wave pattern generated at the stern
ocean currents. But also in harbour circumstances the drift forces may be very important. Although the waves are in general not very high inside harbours, drift forces become large due to shallow water. Thus besides speed and current, bottom effects have to be considered.
Before studying the combination of waves and forward speed, it is impor-tant to understand the forces due to waves and forward speed separately. The forces due to waves only are a special case of the problem we are going to consider, i.e. the forward speed is zero. Thus insight in these forces will be gained within this thesis. The forces in the case when no waves are present, however, are not part of this study, but knowledge about these forces is assumed. Therefore a short explanation will be given about the forces solely due to the forward speed of the ship.
When a ship sails through otherwise calm water, it will generate a steady ship-fixed wave pattern. The generation of the waves leads to a power loss of the ship: the wave resistance. At low speeds, these waves have two
dis-tinct sources: the bow and the stern. The bow waves travel in a V-like shape, the so-called Kelvin-wave pattern. These waves can be observed when watching a water-bird swimming. This Kelvin-wave can be predicted very well using analytical methods. The second kind of waves are gener-ated at the stern, see Figure 1.1. They are very pronounced for sailing-boats at reasonably high speeds. These stern waves can not be predicted very easily; even numerical techniques experience difficulties trying to calculate these waves. To compute the wave resistance at finite speed, it is shown by Raven [25] that the complete non-linear free surface has to be taken into account. The results obtained by his method RAPID are very promising. For an extensive study of the wave resistance, the reader is referred to We-hausen [32].
Figure 1.2: Double-body flow.
Unfortunately, the drift forces we want to study are strongly influenced by the total fluid flow due to the forward speed of the ship (i.e. waves
gener-ated at both bow and stern, and the flow around the hull). Thus it may
seem important to know the stationary wave field. As mentioned before, this wave field is difficult to calculate. Thereforeone has the choice either to solve an extra problem for these waves, or to make an assumption about the influence of these waves. In this study we will make an assumption, because we are interested in forces acting on the ship, and not in the waves
themselves. The most important part of the fluid flow to be taken into
account, is the fluid velocity around the ship. This velocity can be approx-imated using the 'double-body' flow: the .flow past the double body, see
Figure 1.2. The ship and the fluid domain are mirrored in the undisturbed water surface. Thus note that Figure 1.2 still represents a side-view of the geometry. Up to the stern the double-body fluid velocity is a good approx-imation of the general flow which included wavemaking. Because no steady wave pattern is taken into account, the flow is much easier to calculate than the general flow including waves. The double-body flow is therefore a practical approximation of the general flow.
If the speed of the ship is constant, the wave-resistance problem is a time-independent problem. The drift-force problem, however, is time-dependent because the waves and thus the motion of the ship will be a function of time. This complicates the problem considerably. In the past the only possible way to solve these time-dependent equations was to assume the waves and all other quantities to be perfectly harmonic. This assumption reduces the equations to much simpler equations in the so-called frequency domain. At first this may seem a reasonable assumption, because most waves indeed
look rather harmonic. However, at full sea a ship may encounter rather .n,m0.,Innntsvrn
strange and large waves. These waves are definitely not harmonic. Until recently, the only way of predicting the ship's behaviour under these severe seastates was by performing model tests.
With the recent development of large computers, however, it became pos-sible to study the time-dependent equations themselves. The waves are then allowed to be non-harmonic and more realistic sea-states can be simu-lated. Thus much attention to the time-dependent equations has been paid in recent years, and is still being paid, with the ultimate goal of making expensive model tests superfluous.
Besides the time-dependency of the hydrodynamic equations, another ma-jor problem arises when studying ship motions. The equations are highly non-linear, even if forward speed is absent. The equation which has to
be solved is the Laplace-equation, which is linear itself. The necessary boundary conditions, however, are non-linear. On the hull of the ship, for example, we have that the velocity of the fluid should match the velocity of the ship. In other words, the ship should always stay in contact with the ocean. Unfortunately, the velocity of the ship is not known beforehand: it is part of the solution. Besides, the position of the hull is also unknown, thus it is not known where exactly the condition should be applied. A sec-ond non-linearity-arises from the water surface. The position of the surface is unknown, and the condition to be imposed is non-linear as well.
In recent years some attention has been paid to solving the complete non-linear equations. First attempts have been made by for instance Ho-mate [26], Broeze [2] and Van Daalen [1]. The first two of these studies were focussed on the propagation and breaking of water waves, without the presence of an object. The most recent study by Van Daalen concerned the two-dimensional water flow induced by the motion of an object. How-ever. no incoming waves are simulated and current or forward speed are not included in the mathematical model. Because it is not known yet, what the influence of forward speed is on the behaviour of the equations and on the numerical algorithm to be used, we chose to study a linear problem including forward speed. As a next step the non-linear approach and the knowledge of forward-speed effects can be combined in order to allow nu-merical simulation of the non-linear motion of sailing ships.
To obtain a linear problem, we linearize the boundary conditions of the differential equation. We will assume that all waves and motions arising in the simulations are small. Note that this assumption is also made for the stationary waveheight due to the double-body flow; this waveheight is there-fore neglected. Then it is allowed to linearize the water surface around the undisturbed surface, which is a flat plane. The non-linear condition which
has to be applied on the water surface is linearized using the fact that the waves are small. The condition on the hull is linearized around the averaged position of the ship, instead of imposing it on the actual position of the ship.
This linearization process is described in chapter 1. The non-linear equa-lions are summarized briefly, without going into too much detail. (For a more thorough discussion of the principles of water waves and ship mo-tions the reader is referred to Newman [15]. ) Then a perturbation series is used to derive the equations for the forces and moments and the equa-tion of moequa-tion. This derivaequa-tion is an extension of the derivaequa-tion given by Pinkster [20]. The last part of chapter 1 will concern the approximation of the stationary potential by the double-body potential.
To solve the mathematical model a numerical algorithm is needed. This numerical algorithm is discussed in chapter 2. The three-dimensional prob-lem will be reduced to two dimensions by means of the boundary-eprob-lement method. The resulting integral equation is discretized assuming constant quantities on each element. Then it will be shown that the numerical al-gorithm suggested in literature is unstable when forward speed is included, and thus unsuitable for our numerical simulations. A new algorithm will be derived which is stable for all simulations carried out in this study. A suggestion is given how to extend our algorithm to the non-linear problem.
The numerical algorithm is first tested for the case of an infinitely long cylinder, see chapter 3. The integral equation can then be reduced to only one dimension. Convergence of the time-integration will be shown, and the results will be compared with results given by Vugts [31] and Zhao [36]. Further testing is done in chapter 4. The three-dimensional case of a float-ing hemisphere is studied, thus eliminatfloat-ing rotations and moments. An absorbing-boundary condition will be derived by an asymptotical analysis of the Green's function satisfying the Neumann-Kelvin free-surface
condi-tion. Results will be presented for the hydrodynamic coefficients and the drift forces. The results will be compared with Pinkster [20] for zero for-ward speed, and with Nossen [17] for the coupling coefficients when forfor-ward speed is included.
After these two tests, the numerical method is used in chapter 5 to calculate the hydrodynamic coefficients and the drift forces and moments on a com-mercial super tanker. The results for zero forward speed will be compared with Pinkster [21]; the results for non-zero forward speed are new and can not be compared with literature.
time signals. This enables us to calculate the hydrodynamic coefficients in one single simulation. Furthermore it becomes possible to calculate the slow-drift motion of an object, in deep and shallow water. The calculations are performed for the case of the infinitely long cylinder. To ensure accuracy of the simulations an improved absorbing boundary condition, developed and programmed by Sierevogel [28], will be used.
Finally, conclusions will be drawn and recommendations will be made for improvements of the numerical algorithm and the computer programs.
added-mass coefficients [kg]
A asymptotic added mass [kg]
Awl waterline area [m2]
no added-damping coefficients [kg s-1]
: asymptotic added-damping [kg s-I]
Cij -. restoring-force coefficients [kg S-2]
T; asymptotic restoring-force [kg s-2]
P
force [N]Pi i-th order force [N]
Fn Fronde number [-]
En!, grid Fronde number [-]
Green's function [in-i]
11 momentary wetted hull surface [m2]
Ho averaged wetted hull surface [in2]
Kii step-response functions [kg s-2]
Lp Laguerre polynomials of the order p [-]
J17/ moment [N in]
.cli i-th order moment [N m]
Mz3 mass matrix [kg]
R -. radius [m]
T transfer function [kg in-' s-2]
forward speed [in s-1]
V momentary velocity of the hull [m s-1]
V displacement volume [1113]
ij
translation displacement [m]Ci wave velocity [in s-1]
d diameter [in]
g : gravitational constant [m s-2]
h fluid depth [m ]
k wave number [m-1]
ii normal vector [-]
P pressure [Pa]
Po atmospheric pressure [Pa]
-I time [s]
i
local fluid velocity [m s-1]il coordinates [m]
(x, y, z) coordinates [m]
ib point of buoyancy [nil
ig point of gravity [nil
41, total potential [m2 s-1]
fi rotational displacement
total displacement [in]
0 angle of incidence
small parameter [-]
(
wave height [in]Ca wave height of incoming wave [In]
A wave length Tug
field coordinate in Green's function [m]
P density [kg iii-3]
T non-dimension parameter [-]
(I) time-dependent potential [m2 s-'1
71. stationary potential [in2 s-1]
frequency of encounter [s-1 ]
w
co scaled frequency of encounter [-]
wo frequency [s-1]
wo scaled frequency
time step [s]
Ax grid size [in]
Problem Definition and
Mathematical Model
In this chapter the physical problem will be defined and formulated in a mathematical model. Then the non-linear boundary conditions arising, will be linearized by assuming that the quantities involved are small. Finally, equations for the forces and moments will be derived using it. perturbation!
series for both the position of and the pressure on the hull.
1.1
Physical problem
in sailing on oceans, a lot of physical processes influence the behaviour of a ship. Some of these are waves, wind (especially under stormy conditions) and currents. Also the speed of the ship itself is of great importance. The problem we will consider, is the interaction between a sailing ship and the ocean it is sailing at. Especially the influence of waves and current on the behaviour of the ship is considered in detail.
The influence of current and the speed of the ship, is evident. They will cause a resistance, which has to be overcome by the ships engines. Side currents will of course make the vessel drift aside, but this effect can be calculated rather easily, or measured in towing tanks, i.e. the equivalence of wind tunnels.
Model tests show that the influence of waves can be very large. It may cause drifting of the ship, i.e. loosing its course. This effect is widely known and used by kids, trying to get a football out of the water. Throwing stones into the water generates waves which cause the ball to drift slowly into the direction of wave propagation. This effect is measured in model tests, performed in large water basins. Although difficult to perform, these mea-surements are reliable.
However, it appears from measurement that the coupling between waves and current causes a considerable increase in the drift effects. Unfortu-nately, the drift measurements become even more difficult when current is involved. Therefore this effect is subject of mathematical modelling and calculations in this thesis.
Drifting effects are not only important for ships sailing at the ocean, but also for ships manoeuvring close to or in a harbour. Although waves are generally not very high inside harbours, the drift forces may still he very large due to bottom effects. The influence of the bottom on the drift forces will therefore be studied for a commercial super tanker.
Ocean waves can in general be modelled reasonably well by a sine function of one or two frequencies. If two frequencies are involved, the object (ship or a -floating platform) will not only react to those frequencies, but also responses to the sum and difference frequencies. Especially this difference frequency may be in the order of magnitude of the resonance frequency, and may therefore cause significant damage to the structure. This effect
should also be included in the mathematical model.
For the convenience of readers who are not familiar with naval hydrody-namics, a short summary of the most important jargon will be given:
surge horizontal translational motion of the ship in
the sailing direction, i.e. in the direction of the length of the ship.
sway : horizontal translational motion ofthe ship,
per-pendicular to surge.
heave
,.vertical translational motion of the ship.
roll rotational motion round the longitudinal axis,
or the 'surge'-axis
pitch : rotational motion round the 'sway'-axis.
yaw rotational motion round the 'heave'-axis.
point of gravity point of gravity of the whole ship, including load and structures above the water-line.
point of buoyancy : mathematical point of gravity of the hull shape;
point of application of the buoyancy force. drifting : being carried along by current and waves; also
used for moments to indicate that these
mo-ments have the same physical cause as the forces which cause the actual drifting.
1.2
Mathematical model
We will consider an object of general shape in three dimensions, sailing at sea in the presence of waves and current. The object is free to move in all directions or to rotate round any of its axis.
The most obvious choice for a coordinate system would be a system which is fixed relatively to the earth. However, this will cause the object to sail out of the coordinate system. Because we are only interested in processes which act on the object, we would like to have our object within the frame of our coordinate system. Therefore we define a coordinate system moving with the average speed of the body, such that the undisturbed free surface coincides with the plane z = 0 and that the centre of gravity of the body is on the z-axis. The z-axis points upwards, which implies that the fluid domain is in the half space with negative z-values. Note that this coordinate system moves relatively to the earth. This changes the frequency of the waves as described in the model. Therefore we will distinguish between the frequency in the space-fixed reference frame, ok, and the frequency of
encounter, co.
The forward speed is translated into a constant undisturbed horizontal
velocity U of the fluid the object is floating in. The fluid depth is supposed to be constant, h. In order to study drift effects, we assume that regular incoming waves are travelling in the water-surface in a direction which makes an angle B with the positive x-direction, see also Figure 1.1.
LI
Figure 1.1: Aerial view of the geometry.
The goal of this study, as stated in section 1.1, is to calculate the interaction between the ocean and an object. This can be done by calculating the forces
and moments acting on the object. The force is given by pildS
With the normal on the objects surface II. The pressure p can be
cal-culated using the Navier-Stokes equations or Bernoulli's equation, which couple the local fluid velocity with the pressure.
Because the fluid under consideration is water, we assume our fluid to be incompressible. If we furthermore assume the flow to be irrotational at all times, we can introduce a velocity potential (1) given by
= . ( I. )
Note that the latter assumption implies that we neglect the effects of vis-cosity. Using this velocity potential, the Euler equations can be integrated into the equation of Bernoulli:
P Po 84)
1
-=t
+ g z + co n s t . (1.2)This equation will be used calculating the pressure distribution on the ob-ject and the wave heights at the water surface.
In order to use Bernoulli's equation, we need a set of differential equations for the potential including proper boundary conditions. The potential has to satisfy the Laplace equation:
and
which can be derived easily from the continuity equation. This equation represents the physical principle of conservation of mass.
At the water surface we have two boundary conditions:
the pressure should equal the atmospheric pressure (p0), and
a fluid particle can not leave the surface. In mathematical form these conditions are given by
(1.3)
Tt
:v(1) v4) +
+ const = 0 + 54) az P= 2 2. +both at z = (,the unknown momentary position of the water surface, which is given by
1 ( ail> 1-.
--2(1
V) V(I) +
const)These conditions are well known as the dynamic and kinematic conditions. Because there are no exterior restrictions to the water surface, this bound-ary will be referred to as the free surface. To obtain one condition, the
above two are combined, leading to
a24, -a4)
+2V(1)VT-t -2V(1)V V41) V(I)
+9Oz
=0 at
z = (41.4)at2
The floating object is assumed to be impermeable, so no fluid particles can cross this boundary. The normal velocity of the fluid should therefore equal the normal velocity of the object:
a4)
-V on H . (1.5)
Note that slip may occur along the hull, due to the assumption of an invis-cid fluid. The velocity V is the objects velocity in the coordinate system fixed to the average position of the object. In the earth-fixed coordinate system, this velocity has to be superimposed on the undisturbed velocity U.
On the bottom of the fluid domain we have, like on the boundary of the floating object, that no fluid particles may cross this boundary. As the
bottom is fixed, we have
a(I)
= 0 at
z = -h
. (1.6)On
As our free-surface boundary condition is a second-order differential equa-tion in time, we have to provide two initial condiequa-tions. We assume that at t = 0 the fluid is undisturbed, i.e.
41(i, 0) = 0 , (1.7)
and
04)(40)
at = 0 (1.8) To make the solution of our mathematical model unique, we have to impose an extra condition, called the radiation condition. This condition states that waves can only be generated by the body, except for possible incoming waves, and that they should travel away from the body. When studying the equations analytically, this condition can be imposed at infinity. Innumerical studies however, this is not possible due to limited computer time and memory. Therefore this condition has to be given on an artificial boundary.
There are several possibilities for the radiation condition to be imposed, extensively reviewed by Romate [26]. Most commonly used are a damping zone in the free surface or partial differential equations on the artificial boundaries. A damping zone has as advantages that it is easy to implement, and that it has good reflection properties for a wide range of frequencies. The disadvantage is that a large domain is needed to have good absorption. Partial differential equations, however, can be applied closer to the body. The radiation condition would then read
IN>
-57 + 0 ,
t=i
with ci the local phase velocity of the wave to be absorbed and N the num-ber of different waves in the problem. Romate showed that this condition leaves us with a well posed problem, or at worst a weakly-ill posed one. A disadvantage of these equations is that it absorbs only those waves whose
wave velocity is included. Other waves will be partly absorbed, partly re-flected. Recently, a new boundary condition has been developed by Sierevo-gel [28], which has good reflection properties for a wide range of frequencies and can still be chosen very close to the object, thus combining the advan-tages of the previously mentioned methods.
In major parts of this thesis, partial differential equations will be used, be-cause in these problems only one wave has to be absorbed at the boundary. In the last chapter, concerning the investigation of general time signals, the new absorbing boundary condition developed by Sierevogel will be applied.
1.3
Linearization of the boundary conditions
The problem formulated in the previous section contains several non-li near-ities. Unfortunately it is not possible yet to solve these equations, especially in the case of an object in current or an object with a forward speed. In the case of zero current, or no forward speed, first attempts have been made by Romate [26], Broeze [2] and van Daalen [4]. However, major problems arise in solving these non-linear equations, even without forward speed. Including forward speed probably worsens these problems. Therefore we will linearize the boundary conditions, in order to study the problems aris-ing from introducaris-ing the forward speed. Our results and experiences may then be used solving the non-linear problem including forward speed or current. Linearizing the boundary conditions restricts our results to small1.3.1
Free-surface condition
The non-linear free-surface condition is given by (1.4):
a24)
-a4)
(-
ad)+ 2v4)
.vat
+ -2\74) -V VS- VS + = 0 at=(.
at2
This condition is non-linear in two ways. The first obvious non-linearity is in the potential 4) itself. The second non-linearity is hidden in the fact that the position of the boundary, C, is unknown and a function of the
potential 4). To linearize this equation we assume the potential to be a linear combination of a time-independent part and a time-dependent part:
4)(i, t) = d)(i) + t) (1.9)
Furthermore we assume that both the time-dependent and the time-in-dependent potential are small. We assume that for 0 linear terms are sufficient, and.for 0 quadratic terms. Using this in (1.4), we get
.920 + + (tTb t71).) + +
at2 at 2
act) a(j)
gaz + az 0
at z=(
. (1.10)For ç we linearize Bernoulli:
=
(a0
+ -1f'7') '776 + t75. --\,0-
It72)g
at
2Because at infinity there will be no wave elevation, the constant in the equation of Bernoulli is given by --11/2.
To overcome the second non-linearity, the unknown position of the free surface, we expand the condition into a Taylor series around z
. 0. Thus
we assume ( to be small, which is in accordance with both 0and being
small. We get for the time-dependent problem:
a2°
\la°
++t76, t (to t.(7) +
at2 at 2
ao u2)
( 020
1 a30g az+ 2
a"'2 ga2taz)
a2-6 (t"-c-6- \70 2-a = 0 at z = 0 . (1.12)
az2 at
For the time-independent problem we have, up to the second order of the stationary potential, ya715 0 at z = 0 (1.13) az at z = 0 . (1.11) z =
1.3.2
Body boundary condition
On the hull we have, as stated earlier, that the normal velocity of the fluid should equal the normal velocity of the body. Therefore we have
Timman [29] showed that this condition can be linearized into
'On
and
aa
671= n
[(V°. V)
(6. t./.)-(-biBoth conditions now hold on the average position of the hull. Here, d Is the total displacement vector, given by
= + x
41.19
with I the translational, and Si the rotational motion of the body relative to the centre of gravity
1.4
Forces and moments
The goal of this thesis as stated in the Physical problem, is to calculate the interaction between an object and the presence of waves and current. To calculate this interaction, we have to calculate the forces and moments acting on the object due to incoming waves and current. They can be
calculated using Bernoulli's equation. The pressure on the hull is given by
a4)
- p t.(I) gz
u2)
on r.dt 2 2 Kir)
The forces and moments then follow from integrating the pressure ever the body's wetted surface:
fS
(.38)
and 'On =f
p-
g X [(IS . KI.19) (a.a4) (1.15) V ( + =Note that the moments are calculated relative to the centre of gravity, not relative to the origin, see also (1.16).
Unfortunately, the actual position of the body, Ii, is not known, and due to the linearization all quantities are known on the average position. Besides the position of the hull, the size of the integration domain, the wetted part of the bull, is also unknown. Thus some kind of an approximation has to be made to evaluate the above integrals. We will derive approximate formulae extending the method outlined by Pinkster [20]. He used the same kind of perturbation series as we will in the following, but forgot to expand one of the variables. This forced him to correct his expressions along the way, which will not be necessary in the correct derivation. Furthermore an extra term will appear in the drift forces and moments, which was overseen by Pinkster.
We previously assumed that the time-dependent potential is small. Because this potential causes the movements of the hull, the latter is assumed to be small as well. The pressure at the actual position of the body can therefore be expanded into a Taylor series around the average position, Ho:
PH = PH. + PH. + 1
a V PH0
o(ii)
Furthermore the size of the actual wetted surface can be estimated by the sum of the average wetted surface and an oscillatory disturbance of this average surface. We therefore have
ff (g)dS =
f f (Z)d5+ f (E)dS .ho Hose
The integral over the oscillatory wetted surface can be approximated by
ff (i)dS
I f
(i)dzdl . (1.21)Hosc W1C"
Combining these two approximations results in an estimate of the total forces and moments. As may be seen from the Taylor expansion of the pressure, the forces and moments will be proportional to powers of the movement of the object. However, the dependency is more complicated than suggested by this expansion. The pressure on the hull is dependent on the potentials, which in its turn are dependent on the displacement, see equations (1.17) and (1.15). The displacements are on the other hand dependent on the pressure through the forces and moments. Because the movements and the potentials were assumed to be small, the obvious way to proceed is introducing a perturbation series.
(1.20)
1.4.1
Perturbation series
To evaluate the integrals resulting from the estimates made above, we as-sume that all quantities can be written in a perturbation series. The small parameter involved in these series could for instance be the wave height of the incoming wave. This wave height is indeed small, because we assumed that the time-dependent potential is small as well.
The perturbation series for the quantities involved are: Ec(1) + 62((2) + o(E3)
EOM E20(2) + 0 (el
PHOI P-F Ep(i) E2 p(2) (63) = Ey(i) E2)-e(2)+ (7.3) = + 6.211(2) +o (E3) = ;77
-
I ) X (T. E21i(2) X Vi-
+ (E3)5 =
E I x E 2 fi( 2 ) x+ 0 (El
The first terms in the perturbation series are time independent, so
;=
V_U2
is the stationary wave-height. These stationary contributions are assumed to be small, in order to be able to linearize the free-surface condition as has been done in equation (1.12). For the same reason the time-dependent parts of the series are assumed to be small as well. Note, however, that these assumptions are independent from each other and that no difference in order is assumed between the stationary and first-order contributions. Substituting these series into the Taylor series of the pressure on the actual wetted surface, (1.20), and collecting equal powers off, we get:
PH
= P
(1) PH = (1) + {j? (1) + fi (1) x3;)} V P
1 P =4 =
+ + = + +0
+ ={[y(i) + fin)
_ Eg)] v12-73 2From Bernoulli and the perturbation series for the potential we get for the components of the pressure series on the average wetted surface:
= (gzo V-46 *V-46 U2)
IP)
_p (V;;T) Vo(i) 64,1(:))302) 1
p(2) + vo(1) vo(2) vt.))
at 2
The floating bodies of interest for our calculations sail at low speed or float in moderate currents. Therefore, terms of 0 (U2) in p will be neglected in the calculations of the forces and moments.
We substitute these expansions of the pressure into the equations for the force and moment, (1.18) and (1.19). This results in a power series of the forces and moment in E. The stationary contribution is then given by:
0 I1-5HTidS =
0)
(1.22) Ho pgV and = pH ig) X 771 d,9 = pgV Yb (1.23) Ho 0Equation (1.22) represents Archimedes' law. The stationary resistance of the object is of the order U2, and has therefore been neglected. At the end of this chapter, the stationary potential will be approximated by the
double-body potential. For this potential, the stationary resistance is zero due to the paradox of d'Alembert, as shown for instance by Meyer [12]. The zero-forward-speed moment, equation (1.23), is zero if the buoyancy point and the point of gravity of the body are on the same vertical axis. If this is not the case, the object will rotate until this requirement has been fulfilled. Therefore we assume the object to be in equilibrium as far as buoyancy is concerned. (2) PH = P(2) + {y(2) + 11(2) x 5g) fig) [(')x (Tn.
4)] }
vT) {1(1) /1(1) xg)} v
+ + I p + =1.4.2
Hydrodynamic coefficients and equation of motion
When we continue expanding the forces and moments along the way indi-cated above, we find formulae for the first-order quantities:
= (p(111)-#, + rt( )) and
Pg
Pg
, X3(1) Awl S12(1) f o xdS + 121(')Df ydS i(
Df (a0(1)
P G'-c; Vck' ' ridS at HoM'
fpcp;.) /Jo ) X 7pT(1) x Pig) x T-7]
d S =AT) f ydA + t2r) [f g2c1A + (zb zg)1] 1-21) xydA
xdA f xydA-1-52(21) f x2 dA+(zb zg)V _11(11) sl(21) f ao(1) ).t.(1)) at + Ho (1.24) 279) x TidS (1.25)
with D the water-line surface of the object, and Awl the area of D. Note that the integral over the oscillatory wetted surface does not contribute to the first-order quantities. However, if we would have included second-order effects in the velocity in the derivation, the first-order quantities would have been affected. For instance, the first-order force would be corrected with
pg cdl.
wl
As stated before, these terms are not important in the applications we are interested in. dS F +
f
-,As can be seen from both (1.24) and (1.25), only few terms depend on the first-order potential directly. Most terms depend on the first-order motions in combination with stationary quantities. However, these motions can only be calculated when the first-order forces and moments are known. This problem is very well known in calculating the displacement of a mass connected to a string. The internal force of the string is proportional to the displacement of the mass. To overcome this problem, the restoring force is taken as part of the differential equation. Only external forces, independent of the displacement, are considered as forcing terms. Analogously, the terms depending on the first-order motion are considered to be restoring forces and moments in the equations of motion, instead of being terms of the first-order forces and moments. If we assume the object to be symmetrical in the y-plane, the restoring-force coefficients are given by:
C33 pg
I
dAC35
pg xdA
C44 pg{(zb 2.9)V + y2dA
C53 xdA
C55 pg{(zb zg)V x2dA}
The other elements of the restoring force matrix C are zero, some due to the symmetry assumption. The restoring coefficients are independent of the forward speed, up to the second order. For ships sailing at higher velocities, these coefficients must be adapted.
Unfortunately, the terms depending directly on the first-order potential are not entirely independent of the motions. This may be observed from (1.15), where the potentia cb is coupled with the motion To calculate the forces on the object these dependencies have to be removed. To achieve this we introduce added mass (Aii) and added damping (Bii) coefficients. These coefficients represent the dependency of the force on the motion, or vice versa. For an extensive explanation of the principle of added hydrodynamic coefficients, the reader is referred to Ogilvie [18].
Because the forces and movements depend on each other, we first force the object to oscillate in a given direction. The forces thus calculated depend directly on the motion of the object. For purely sinusoidal motions the
+
=
added mass and damping coefficients are given by lifting the forces to the acceleration and velocity of the body:
a2.
ax;= =A.. ot2' 1323 at
Here j is the direction of the motion, either translational or rotational, and i the direction of the force. Tins can, of course, also be done for the moments. The coefficients are then still called Ai; and Ilj but now with
i E [4, 5, 6].
When these coefficients are known, we calculate the forces due to the incom-ing wave field, for a fixed object. Then the motion of the object is calculated using the equations of motion and the hydrodynamic coefficients:
i(M + A)
a2f
aC + CY =-
(P.Inc
Aline)
at2 at (1.26)
with (XI, )(2, X 31Cli, 92,113)T. The mass matrix Al is diagonal and consists of the mass and the relevant moments of inertia.
For the calculation of the hydrodynamic coefficients in the case of general time signals, the reader is referred to chapter 6.
1.4.3
Drift forces and moments
In studying drift effects, only the time-average values of the forces and mo-ments are of interest, not the oscillatory forces themselves. If we assume the incoming wave to be a sine function with zero mean value, it may be seen from (1.24) and (1.25) that the mean value of the lirst-order quantities are zero. This is. of course, caused by the fact that these forces and mo-ments have a linear dependency on the time-dependent quantities 0, X and SI Thus second-order forces and moments will be needed to calculate the drifting of an object. In general sea keeping, however, the incoming wave may be seen as a finite sum of sine functions. Still, the first-order forces and moments will have zero mean value and second-order effects have to be taken into account.
By continuing our perturbation series, we arrive at formulae for the mean values of the second-order forces and moments:
(p(2)) ( 001) at VO(1) V71)7i:d5 (Y(31))2Tidl
+0)
X M.V(1) p ( ) V-0(1)fidS 2J
Ho 0P9(
0 C2(1')/21)C3,5Note that several terms with zero mean value have been left out. This
formula is equivalent to the one given by Pinkster [20], except for the last vector which is missing in his thesis, and the forward speed contributions. The last vector represent the force generated by the displacement of the point of gravity of the hull deck due to two combined rotations.
For the average second-order moment we find
1
tO(1) f0(1)
zg) X TidS pI
(6(1) (150(t1) 0(I ) ) x 777.d,c Ha+I
pg (01) et:(31))2(z ig) X 2 wl _12(1)n(i),-, to-41)c,552 a1-3 4,1 )Again, the last vector is missing in Pin kster's thesis. As may be seen from the above equations, the formulae for the drift forces and drift moments are almost equivalent. Not regarding the last vector, translational motion is equivalent to rotational motion, and the normal
to ri
Eg) XNote that the integral over the oscillatory wetted surface does contribute to the second-order quantities, even now that we did not include second-order effects of the velocity. Of course, including those effects would give rise to much more contributions of this integral. However, these terms are not important for objects sailing at low to moderate speed.
(1.28)
)
)
(1.27) wil (Kti(2)) =I
pg1.5
Double-body potential
In the section concerning the linearization of the boundary conditions, sec-tion 1.3, the potential has been split up into a stasec-tionary and an instasec-tion- instation-airy part:
(D(E, t)= 0(i) (/)(E,t) .
The stationary potential includes the resistance of the object in calm wa-ter and the Kelvin wave patwa-tern; the instationary part includes all time-dependent effects. Unfortunately, the stationary potential is difficult to calculate. However, from (1.13) and (1.14) it can be seen that up to the second order in the velocity U, the normal derivatives of the stationary potential equal zero. Therefore, no waves are created up to this order in our linearization. This allows us to approximate the stationary potential with the double-body potential: the potential arising from the case when the whole computational domain is extended with its mirror image in the plane z = 0.
One would expect that taking into account higher-order terms in the ve-locity U would lead to the well-known Kelvin wave pattern. Unfortunately, this is not the case. To find this stationary wave pattern, the method has
to be extended in another way. A first possibility is solving the complete non-linear system. This has been done by Raven [25], using a Dawson-like approach. However, these calculations are very time-consuming. A faster way has been developed by Sakamoto [27]. He linearized the free-surface condition around the wave-height resulting from the double-body potential. He then found wave-like corrections to this double-body poten-tial. Hermans and Van Gernert [7] showed that Sakamoto's choice of the double-body potential was the only correct one to yield solutions resem-bling the Kelvin-wave pattern. The disadvantage of this linearization is
that it is only applicable to low speeds.
In the derivation of our linear free-surface condition we linearized around = 0. Thus no wave-like solution can be found for the stationary potential. This may only be correct when these waves are small, i.e. whenthe speed of the object is low. This is indeed the case for the commercial supertanker for which we will perform calculations in chapter 5.
For simplified geometries, the double-body potential can be calculated an-alytically, which means that the derivatives needed in the boundary con-dition of the instationary potential are exactly known. This enables us to eliminate the errors arising from numerical differentiation of the stationary potential.
Numerical Method
In this chapter a numerical method will be given to solve the equations as given in the previous chapter. The differential equation will be dis-cretized using a boundary-integral 'method. Then it is shown that the time-integration algorithm known in literature is unstable for problems including forward speed. A new algorithm will be given, which is stable for all time steps and space discretizations.
2.1
Introduction
In the previous chapter, we derived a mathematical model for the interac-tion between current and waves, and an object. The model consists of a governing equation with linearized boundary conditions. The wanted in-teractions, i.e. the forces, all follow from integrals over the wetted surface of the object, as soon as the potentials on this object are known. To calcu-late these potentials, we have to solve the Laplace equation. This equation itself is time-independent. However, the boundary conditions do depend on time, so sonic sort of time integration has to be carried out. Therefore it will be useful to choose a fast solver of the Laplace equation.
There are several ways of solving the Laplace equation. It can be dis-cretized straightforwardly by using a finite-difference or a finite-element method. These methods have the advantage of discretizing the differential equation directly, and generating a sparse matrix. However, since we are only interested in values of the potential or its derivatives on the surface of the object, there is no need to discretize the entire fluid. Still doing so will lead to very large, though sparse matrices. Especially for computations in
three dimensions, this can be a big problem.
An other method to solve the Laplace equation is using a boundary-integral method. In this method, only values of the potential on the boundary are
needed. This reduces the number of dimensions to two, and therewith
shrinks the matrices with a factor N, the amount of elements in one direc-tion. Unfortunately, this method has the disadvantage of generating full matrices, which in its turn will increase computer time needed to solve the matrix system.
Although the advantages and disadvantages of the methods may be bal-anced, we chose the boundary-integral method, as has been done by many authors in recent years. A major consideration was, that computer time is less limited than computer memory.
The time integration of the equations is a major problem in itself. Until
quite recently, all calculations were performed in the frequency domain, thus avoiding the time-integration problem. However, the limitations of the frequency-domain approach forced the investigations to be focussed on time-domain analyses. Most of the recent studies use a rather obvious and straight-forward algorithm to perform the time integration. Rotate [26], for instance, used this algorithm to calculate non-linear wave propagation. However, the algorithm used in literature is not suitable for problems in-cluding forward speed, as will be shown in this chapter. Thus a new algo-rithm has to be developed.
2.2
Boundary-integral method
The boundary-integral method is based upon Green's second identity. It
reformulates the differential equation into an integral over the boundary, thus reducing the number of dimensions of the problem. It can be shown that if the boundary of the fluid domain is twice continuous differentiable, and if the function (15 is twice continuous differentiable inside the domain and once continuous differentiable on the boundary of the domain, the
following identity holds, see Colton [3]:
G (7, a (-.)
(4'(±.)
I
(0
(t)
,;()
dS . (2.1)
dV
Here G("i,e) is the Green's function of the Laplacian, and the normal pointing out of the fluid domain. Note that in this representation the Green's function should represent a source, not a drain. The Green's func-tion depends on the number of dimensions of the problem, and will therefore be given in the consecutive chapters. Now the problem has been reduced to solving this integral equation together with the boundary conditions of the original differential equation.
Et should be noted that the boundary of the fluid domain used in the cal-culations is not twice continuous differentiable, thus violating the demands of Green's theorem. However, all collocation points are on allowed parts of the boundary. The used boundary can be seen as an approximation of the real boundary which is twice continuous differentiable.
To solve the integral equation, we discretize the boundary by dividing it into panels. On such a panel we assume the potentials to have constant values. This simplifies (2.1) to
N
=
f
E
G,i=1 ,..,N
On
"
= e3 J=1 3e,
with N the number of points in the discretization. The integralsover each
separate element are calculated analytically, using results from literature. The formulae used will be given along with the Green's function.
We could also have used a higher-order panel method, assuming the quan-tities on an element to be linear or quadratic. This has been done bymany authors, see for instance Romate [26] and Nakos [13]. Of course this will increase the accuracy, if the number of panels is chosen the same. However, the number of unknowns increases with a factor 3 or 6, due to the extra unknowns per element. This would mean that, given a maximal number of unknowns due to limited computer memory, the number of elements
has to be decreased by this factor. To overcome this problem, Romate
eliminated two of the three unknowns by using finite differences for the tangential derivatives prior to solving the matrix equation. If this is done accurately enough, the order of the panel method can be maintained. But the integrals involved in discretizing the integral equation are much more difficult to calculate and will absorb much more computer time. Therefore we chose to use constant elements and differentiate the spatial derivatives
n u tnerically afterwards.
Discretizing the boundary integral in the way outlined above, results in a matrix equation which can be writtenas
Aj)(t)
= BZi(t)
(2.2)where the potentials and the normal derivativesare still functions of time.
2.2.1
Discretization of the boundary conditions
The governing differential equation has been discretized in the previous section by using a boundary-element method. The boundary conditions belonging to this integral equation have to be discretized as well. This involves discretization of both time and spatial derivatives of the potential. The discretization of the time derivatives is part of the time-integration procedure and will be discussed in both sections 2.3. I and 2.3.2.
The spatial derivatives of the potential are discretized by using the closest neighbours of an element in the discretization of the free surface. These neighbours are known from the mesh-generation. Because the steady poten-tial has been approximated by the double-body potenpoten-tial, no z-derivatives of the potential are involved in the free-surface condition. Therefore we have to discretize derivatives in two dimensions only. We used five neigh-bours and the element itself to calculate the derivatives. This implies that the first-order derivatives were estimated up to second order in the grid-size, while the second-order derivatives were estimated up to first order.
Differentiating the potential numerically may seem to be contradicting the assumption of constant quantities over an element. However, these con-stant values should be seen as samples of an otherwise smooth function. This function is assumed to have continuous derivatives up to the neces-sary order. Thus the derivatives may be approximated by a finite difference scheme acting on the samples of this function.
The boundary conditions also contain derivatives of the steady potential, approximated by the double body potential, as mentioned in section 1.5. These derivatives are calculated analytically if possible, or numerically. When calculated numerically, the method outlined above has been used.
2.3
Time integration
The time integration has to be performed on a time-independent matrix equation combined with time-dependent boundary conditions. The matrix
equation is the discrete equivalent of the Laplace equation:
Aq(t) = Bic(t)
.We first consider the boundary condition on the free surface. This condition contains several derivatives of the double-body potential. However, rela-tively far away from the object this potential approaches the undisturbed-flow potential. Thus the free-surface condition may be studied using the
potential' Ut, under the assumption that the object is far away or no body is present at all. The free-surface condition can then be written as:
a2(i.) 820 a215 ao
at2 2Q- 2axatax2 + g az
o (2.3)
which can be discretized into
a2,7 ae-i;
+
C2-E
at2 at
In this condition the spatial derivatives have been discretized as outlined in the previous section; the vector f represents boundary effects in the dis-cretization.
The absorbing boundary condition contains time derivatives of the poten-tial as well and is treated in the same way as the free-surface condition. The boundary condition at the objects surface has no time derivatives and is therefore not involved in the time integration.
2.3.1
Algorithm for zero speed
In literature, several algorithms are suggested for integrating the time-differential equations, see for instance R.omate [26] and Nakos [14]. They all have the same underlying principle, which will be outlined briefly.
prescribe 0 or at the boundaries at given time level to solve the discretized Laplace equation (2.2) at this time level
integrate the time-dependent boundary conditions in time, to get new values of cb or On at the boundaries at time level to + At
repeat these steps until desired time level T has been reached
The difference between the algorithms of Romate and Nakos lies in the fact that they use different time-integration methods for the time integration of the boundary conditions.
As shown by Romate, this algorithm gives good results and is stable for wave problems without body or current. Van Daalen [1] showed that this is also the case for problems involving a floating object. However, this al-gorithm breaks down when using it in problems involving forward speedor
current. This may be seen from the following simple exercise concerning the eigen-values of the differential equation.
When there is no forward speed, the free-surface condition is given by
820
ao
o
at2 an
If we assume the bottom to be infinitely far away, the normal behaviour of the potential is e. Substituting this into the boundary condition gives
a20
gicch o at2
The analytical eigen-values of this equation are
A = +i\./gk .
10
I m(lambda)
Figure 2.1: Eigen-value of the discretized free-surface condition without forward
speed 1111111111111111111111111111111111111111 1111111111011111111121111111111111111111 1111111111111111111111111111111111111111 1111111111111111111111111111111111111111 1111111111111011111111111111111111111111 111111111111.111111111111111111111111111 1111111111111111011 11111111E1111111M 1111111111111111111111111111111111111111 11111111111111111111EINIENNININNENNEENNI 11111111111111 111111111
mum
IIN IIIIININIINIII 11/11/1H IIIIIENEN
11.1111111111111111111111111111111111111 11111111111111111111111111111111111111111 111111111111111111111111111111111111111111 usIII11111111111111111111111111111111111 1111111111111111111111111111111111111111 11111111111111111.11.1111111511111m1111 1111111111111111111IMINumENNINIMIN ENNummommIllom11111111ENNommulm 1 111111111111111 1111 I I II -2 0 2 4 Redambda 4
Figure 2.2: Eigen-values of the discretized free-surface condition including forward speed U = .03. 1111111111111111111111111111111111111111 1101111111111111111110111111111111111111 1111111111111111111111111111111111111111 121111111111111111111111.1111111111111111111111111111111111111.1111111111111111111 11111111111111111M1111111111111111111 1111111111111111111M111111111111111111 INNEN INN IN NEE NEVIN NEENNENNEENNIN
11111 III II III EMI 11111111111111
1111111111111 11111E U1111111111111111 111111111111111111MAINNINNEN11111111 1111111111111111111§ E111111111111111111 111111111111111111IN E1111111111111111111 11111111111112111111111111111111111111111 1111111111111111111111111111,01111111111 ININSIIIIM1111111111111110.1111110111 1111111111111111111111111111111111111111 11111EmNUMM111111111111111111111111110 1111111111111111111111111111111111111111 111111111111 11/11/11 11111 11/11/1111 -6 -4 -2 0 2 4 6 Reflambda1 0 f/e(lambda)I 4 6 6 -4
tr)
Figure 2.3: Eigen-values of the discretized free-surface condition including forward
speed IT = .045.
Now the water surface is discretized using an equidistant grid, thus yield-ing a matrix equation. Figure 2.1 shows the inverse of the characteristic determinant of this equation, with the peaks representing the eigen-values of the discretized free-surface condition. They were calculated for a one-dimensional free-surface; the wave-number k was based on w = 10, which is a middle frequency in chapter 3. It may be seen that the eigen-values are indeed purely imaginary, and close to the analytical value of A = 10. Figure 2.2 shows that in the case of forward speed, some eigen-values are shifted along the imaginary axis, as may have been expected, because the frequency of encounter depends on the velocity. These eigen-values were calculated assuming that the steady potential equals the undisturbed po-tential (Ix. The derivatives where approximated by their standard second-order difference schemes.
However, Figure 2.3 shows, that for higher velocities the eigen-values are shifting from the imaginary axis onto the real axis. Unfortunately, some of these eigen-values have positive real part. Note that this forward speed is not very high. For the two-dimensional problem considered in chapter 3, this would mean a Froude number of 0.06, which is quite low compared to the maximal Froude number used in that chapter Fu = 0.14.
The effect of the speed and the grid size on the eigen-values is shown in Table 2.1. When no forward speed is present, the effect of the grid size on the eigen- values is negligible. However, when forward speed is included, reducing the grid size causes the eigen-values to become real and positive. This implies that, whatever speed is used, reducing the grid size will even-tually make the differential equation unstable. Note that this conclusion may be drawn independent of the time step, because no time discretization
111111111111111111111111111111111111111111111111 11111111111 1111111111111111111111111111111 11111111111 1111111111111111211111111111 1111111111111111111111111111111111111111 1111111111111111111111151111111011111111111 111111111111111111111111111111111111111 1111111111111111111111111111111111111111 11111111111 111111111111111111 11 11111 11111111111 111111111111111111 11 11111 11111111111 111111111111111111 111111111 1111111111111111111111111111111111111111 11111111111111111111111111111111111111111 1111111111111111111111111111111111.11111 11011111111111111111111111111111111111111 4 INIMS11111111111111111111111111111111111 1111111111111111111111111111111111111111 111111111111111111111111111111211111111111 11111111111111111111111111111111111111111 1111111111111111111111188111111111111111
11 1111 5.1184 INA Off MINN 11/11/11 -6 -4 2 0 2 4 6
Table 2.1: Eigen-values of the free-surface boundary condition for several speeds and grid-sizes.
has taken place yet. These results suggest that the stability is a function of the grid Froude number, Fnh - 11Although the smallest grid size used in Table 2.1 may seem very small, this same grid size is used in the calculations in chapter 3 for higher frequencies.
The fact that some eigen-values have positive real part, implies that the differential equation has exponentially increasing solutions, which are
ob-viously unphysical. Numerically integrating this boundary condition will therefore give rise to instabilities, simply because the equation itself is un-stable. Thus another algorithm has to be found to integrate the time-dependent equations, without integrating the boundary conditions
them-selves.
It may seem surprising that condition 2.3 can not be integrated stably.
After all the condition can be rewritten as
a
0)2
(
aoT)i+u,
gz
= 0id x
The term between the brackets can be seen as a derivative in a moving
frame of reference, and condition 2.3 can therefore be transformed into its zero-speed equivalence. Thus it should be possible to integrate the Neumann-Kelvin condition stably. However, it should be noted that this is only possible using the substantial derivative explicitly. In the above analysis this has not been done, because the free-surface condition used in the actual computations of the next chapters, condition 1.12, can not
U = 0 U = .0.5
11 = .03 It = .015 h = .0075 h = .03 h = .015 h = .0075
10.04i 10.07i 10.14i 12.03i 13.28i 15.32i
10.04i 10.06i 10.12i 11.75i 12.65i 13.67i
10.03i 10.05i 10.09i 11.34i 11.71i 11.17i
9.98i 9.95i 9.88i 10.84i 10.60i 10.57i
9.98i 9.96i 9.90i 10.57i 10.57i 8.13i
9.99i 9.97i 9.93i 10.30i 9.46i 4.78i
10.02i 10.03i 10.05i 9.79i 8.41i +1.88
10.00i 9.99i 9.97i 8.81i 7.56i +4.21
10.01i 10.0-1i 10.01i 9.02i 6.55i +5.31
10.01i 10.01i 10.01i 9.36i 6.93i +4.99
+
'
be rewritten in terms of substantial derivatives. Thus the stability of the Neumann-Kelvin condition has been analysed using straight-forward dis-cretizations of the separate terms.
Nakos [14] still uses the above given algorithm, in spite of presence of for-ward speed and the use of condition 1.12. He states that he had troubles with instabilities, but has overcome them with a filtering technique. As may be concluded from the above, Nakos filtered out instabilities which where inherent to his equations and his way of integrating them. This leaves him with no guarantee that his solution is correct, because filtering affects the solution of a problem considerably. Furthermore he states that his algorithm can only be stable for decreasing grid sizes, if his time step is reduced simultaneously. We think, though, that his algorithm is only sta-ble, due to numerical damping of his time integration scheme. Moreover, a scheme which becomes unstable when reducing the grid size without re-ducing the time step, is highly undesirable. Therefore we oppose the use of this algorithm.
The failure of the above algorithm may be explained by the following : ea-soning. The Laplace equation has harmonic functions as solutions. The boundary conditions serve to determine what harmonic function is the over-all solution of the complete system. Solving the Laplace equation and the boundary conditions separately, can be viewed as satisfying the boundary condition with a harmonic forcing function. The solution will then be a superposition of a harmonic function and the eigen-functions of the differ-ential operator. In the absence of forward speed, these eigen-functions are harmonic as well, and determine the dispersion relation. However, in the presence of forward speed not all eigen-functions of the free-surface condi-tion are harmonic; some are exponentially increasing and some exponen-tially decreasing. This means that at every time step an error is introduced which disturbs the accuracy of the solution. This also explains why
fil-tering may be successful. It also suggests that a better algorithm may be found combining the Laplace equation and the boundary conditions, thus imposing harmonic solutions to the boundary condition.
2.3.2
New algorithm for forward speeds
As concluded in the previous section, the boundary condition at the free surface can not be integrated separately. However, assuming that our model is a good representation of physics, all equations together have a stable so-lution. This indicates that we have to combine the Laplace equation and the free-surface boundary condition to overcome non-physical instabilities.
!["Setnik oratorów polskich XIX i XX wieku", Antoni Bednarek, Lublin 2008 : [recenzja]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==)