Technical Report:
INVESTIGATION OF THE NONLINEAR DYNAMIC BEHAVIOR OF A
“TUMBLEHOME-TOP” SHIP THAT UNDERGOES SURF-RIDING AND
BROACHING IN STERN QUARTERING SEAS
Authored by:
K. J. Spyrou1, K. Weems2 & V. Belenky3
1
National Technical University of Athens
Science Applications International Corporation
3
American Bureau of Shipping
2
The presented investigation was carried out during the first author’s stay at the American Bureau of Shipping
(ABS) in Houston, Texas, on sabbatical leave from the National Technical University of Athens. The work is part
of a project on the probabilistic assessment of ship stability that is currently funded by the Office of Naval
Research. The project partners are: Science Applications International Corporation and the American Bureau of
Shipping.
June 2008
American Bureau of Shipping
Corporate Technology Division
Department of Research and Product Development
ABS Plaza, 16855 Northchase Drive
Houston, TX 77060 USA
EXECUTIVE SUMMARY
Project objectives
The general objective of the current work has been to prepare the ground for a subsequent
implementation of a probabilistic methodology of ship assessment with respect to
broaching. The specific objectives were:
To determine whether the numerical code LAMP-3 could capture phenomena of
nonlinear dynamic behavior of a ship in extreme stern quartering seas that are identified
in the literature as “surf-riding” and “broaching”.
To implement, in an appropriately adapted version of LAMP, advanced numerical
techniques that would enhance the efficiency of LAMP-based investigations of surfriding and broaching. Automated stability analysis of surf-riding in quartering seas
should also be implemented.
Achievement with respect to above objectives
The findings support the argument that, against current theory, LAMP reproduces
reasonably the qualitative pattern of the capture of a ship into surfriding and the
instability of broaching. However for quantitative predictions the model will need to be
verified further.
Continuation analysis of surf-riding in quartering seas in six degrees of freedom, with or
without active rudder control, was implemented successfully.
Project recommendations
The compatibility of hydrodynamic modelling assumptions between LAMP-0 and LAMP3 needs to be investigated further in order to ensure the consistency of predictions by
the two variations of the code.
The maneuvering model used in LAMP should be calibrated against experimental data
before being used for the prediction of surf-riding and broaching when there are
aspirations of quantitative accuracy.
The possibility of developing a guideline for avoiding the encounter of these phenomena
seems feasible and it is worthy of further consideration.
2
TABLE OF CONTENTS
1. INTRODUCTION
4
2. APPROACH ADOPTED
5
3. THE EMPLOYED NUMERICAL CODES
6
4. THE INVESTIGATED SHIP
7
5. NONLINEAR SHIP DYNAMICS CAPTURED BY FOCUSED SIMULATIONS WITH LAMP-3
8
Main findings for Fn=0.30 (RPS=2.0)
9
Findings for Fn=0.33 (RPS=2.2)
17
Fn=0.36 (RPS=2.4)
23
Fn=0.38 (RPS=2.6) and Fn=0.41 (RPS=2.8)
23
Further parametric study
23
Summary of findings of simulation studies
24
6. INTRODUCTION TO CONTINUATION METHODS
32
7. MATHEMATICAL BASIS OF THE APPLIED CONTINUATION SCHEME
33
8. LAMP-BASED IMPLEMENTATION OF CONTINUATION
35
What about the hydrodynamics?
39
9. SAMPLE LAMPCONT RESULTS
40
3-DOF surf-riding in following seas
41
6-DOF surf-riding in stern oblique seas
42
10. LAMPCONT: SUMMARY OF STATUS AND FUTURE WORK
47
11. CONCLUSIONS
48
12. REFERENCES
50
3
1. INTRODUCTION
Surf-riding and broaching represent extreme types of ship behavior that are likely to occur
in steep following/quartering seas, at moderate to high speeds. Several factors have not
allowed, at least up to now, the development of straightforward design or operational
guidelines that would rule out the possibility these phenomena to be encountered during
operation. For commercial ships they are currently dealt with rather superficially, through
the recently revised IMO’s Guidance to the Master (IMO 2007).
The subject has a long history as a research topic. Key earlier references to be consulted are:
for surf-riding, Grim (1951 & 1983); Makov (1969) and Kan (1990). For broaching one should
read, at least, Weinblum & St Denis (1950); DuCane & Goodrich (1962); Wahab & Swan
(1964) and Motora et al (1982). In the last 15 years however there has been some significant
progress towards understanding the dynamical origins of surf-riding and broaching. Detailed
descriptions of these new findings have been presented, for example, in Spyrou (1996a, b, c;
1997). Also, they were neatly summarised in the SNAME book of Belenky and Sevastianov
(2007). Nevertheless, our ability to perform successfully quantitative predictions, especially
of broaching, still lags behind our level of understanding of the phenomena and it seems to
depend critically on progress in at least two areas:
On the one hand, there is need to expand towards incorporating more advanced
numerical techniques (to put it plainly, there is need to go beyond simulation) for
exploring multi-dimensional ship dynamics and for identifying efficiently all possible
critical types of behavior that might be hidden in the dynamical system under
investigation.
Secondly, confidence on the capabilities of the numerical hydrodynamic models to
reproducing reality at a reasonable level is needed to be enhanced further.
More deeply, one should realise that in fact these two areas are interconnected:
improvement of our efficiency in exploring ship dynamics will have undoubtedly an effect on
our ability to improve also the employed hydrodynamic model itself.
Looking at the wider picture, among current developments in the international field of ship
stability one should single out the recent placement of emphasis on the establishment of a
sound probabilistic method for stability assessment. Actually, the current work is part of an
effort of this type, although the specific interest has been exclusively on tumblehome type
ship forms. The hydrodynamic model that was available for this investigation was the well
known code LAMP.
From the foregoing it appears that a task to be inevitably faced soon (and in fact not only
within this running project but internationally), is the smooth integration of nonlinear
dynamical analysis of surf-riding and broaching within a practical probabilistic framework of
4
ship stability analysis. In preparing the ground for this, it is believed that there are tasks of
deterministic dynamics that need to be addressed appropriately before the non-trivial step
of integration is undertaken. In the following will be discussed these pending tasks and the
work that has accrued from them. Furthermore, the progress achieved during this fourmonth investigation will be presented.
2. APPROACH ADOPTED
The questions that governed the directions of this investigation are spelt out below:
Is the existing theory of broaching and surf-riding a reliable basis for practical
probabilistic assessment?
How does the “tumblehome” ship perform with respect to the theory, when assessed
through LAMP?
Is it feasible to implement (in LAMP) advanced numerical techniques for more efficient
investigation of surf-riding and broaching?
Can we specify the critical conditions of surf-riding and broaching that should betargeted
for probability calculation?
Accordingly, the devised investigation strategy specified two areas where work would be
deemed as essential and shows potential to producing immediately definitive
improvements:
a) The ability of LAMP to capturing the tendency of the tumblehome for surf-riding and
broaching should be clarified. This could be accomplished through running a campaign
of carefully selected simulation scenarios. Thereupon, the degree of coincidence with
the qualitative pattern predicted by the current theory of broaching could be assessed.
b) Steps should be taken towards implementing in LAMP advanced numerical techniques.
The one that seems feasible to implement at this stage is Continuation of stationary
states, interfaced with some “light” version of LAMP; specifically a version that could be
treated like a system of ordinary differential equations. Combined with simultaneous
stability analysis, this should expedite tremendously the exploration of surf-riding in
quartering seas which is closely linked with broaching.
5
3. THE EMPLOYED NUMERICAL CODES
LAMP incorporates several different formulations for the solution of the wave-body
hydrodynamic interaction problem. Two of these formulations were used for this work:
LAMP-3 uses an approximate body-nonlinear 3-D, 6 d.o.f. potential flow hydrodynamic
solver specially formulated for large lateral motions. Incident wave forcing and hydrostatic
restoring forces are calculated of the instantaneous wetted hull surface. Radiation and
diffraction forces are calculated over the mean wetted surface; however there are no
assumptions on constant forward speed and small lateral motions.
Hull lift forces are calculated considering the hull as a lifting surface of extremely low aspect
ratio or using lift coefficients evaluated using another potential code VoLaR (Vortex Lattice
Rationale). Vortex-shedding induced drag is modeled in a similar fashion to hull lift forces.
Wave component comes out naturally as a part of potential flow formulation.
Viscous drag force as well as yaw and sway viscous damping are implemented in a
conventional way using empirical coefficients. Conventional coefficient-based model of
propeller is used for thrust. As a result, forward speed is a result of calculations rather than
input figure, with number of revolutions of a propeller as an input. A more detailed
description of force model can be found in (Lin et al 2006).
While the LAMP-3 formulation represent a very reasonable compromise between realistic
modeling and computation efficiency, the memory effect related to radiation forces does
not allow to a vessel to be modeled as a dynamical system described by ordinary differential
equations. Nonlinear dynamics is essential in making a connection between numerical
simulation and the contemporary theory of broaching (Spyrou 1996a, 1996c, 1997). While
some connection is being made by a qualitative comparison of theoretical results and
numerical simulation, the application of Continuation method, at the present state-of-theart, does require that the dynamical system is represented as a system of ordinary
differential equations. The most direct way for dealing with that for the LAMP
implementation is, in the first instance, to remove the memory effect. It is noted that, for
the stationary states of surf-riding that are targeted by the Continuation method (to be
discussed in detail later in this report), hydrodynamic memory is not expected to change the
position of stationary states in state-space, it may however have some influence (one
expects not very substantial) on their stability properties. Therefore, a version of LAMP
called “LAMP-0” was developed (specifically for running the Continuation method –
simulations were still run with LAMP-3), in which the hydrostatic and Froude-Krylov forces
are evaluated on the instantaneous submerged body, while radiation forces are ignored and
added mass is implemented as constant coefficient. Wave-related damping and drag forces
now have to be included via external models or coefficients, in a similar fashion to viscous
effects, as they are no longer evaluated internally. The maneuvering model of LAMP-3
appears intact in the LAMP-0 implementation.
6
4. THE INVESTIGATED SHIP
The ship configuration used for testing and demonstrating the surf-riding and broaching
phenomena is the “tumblehome top” ship form (Fig. 1), derived from the ONR Topsides
Study (Bishop, et al. 2005). This as well as other variant forms have been tested at the Navy
Surface Warfare Center/Carderock Division (NSWC/CD). It is typically referred to as number
5613-1.
Fig. 1: ONR topside tumblehome hull (5613-1).
Compared to a more ordinary “flared top” vessel that has also been used in the apst, the
two have nearly identical underwater geometry, differing only by the angle of the transom
surface, which is similar to the existing DDG-51 combatant. The above-waterline shape of of
the tumblehome hull form is representative of the type of low-observable ship designs that
have given rise to the Navy’s concern about dynamic stability and capsizing. Table 1 lists the
principal dimensions and hydrostatics of the ONR Topsides tumblehome hull form (5613-1)
as specified by the NSWC/CD model test report and as realized in the LAMP geometry
model.
Table 1. Principal Dimensions of ONR Topside Hull Form
Model test results are available for this ship form, including roll decay tests and response in
regular waves.
7
5. NONLINEAR SHIP DYNAMICS CAPTURED BY FOCUSED SIMULATIONS WITH LAMP-3
From a preliminary numerical investigation it was established that: for a wave-length-toship-length-ratio equal to 1.0, wave steepness 1/20 and an exactly following sea, the
tumblehome vessel begins to exhibit surf-riding behavior from approximately Fn=0.30.
However, this represents a necessary but not a sufficient condition; because the adopted
type of response depends also on initial conditions. From that Fn value and to some
distance upwards, one should anticipate the occurrence of interaction phenomena arising
from the coexistence of at least two different attracting states: the periodic state where the
waves overtake the ship; and the stationary state of surf-riding (stationary is meant here
with respect to a wave crest). In general, a higher initial surge velocity makes capture into
surf—riding more likely. Why this happen could be grasped on the basis of the phase plane
diagram of surging and surf-riding that is shown in Fig. 2.
In the following we present our findings for Fn=0.30 and higher, until the Froude number
that corresponds to wave celerity is reached. It ought to be said however that, ship
behavior at lower speeds (specifically at Fn=0.21 and Fn=0.25) has also been examined.
Nothing different from an ordinary surge response has been found at those speeds, for
commanded heading angles up to 20 deg. Furthermore, vastly different initial conditions
seemed to converge always towards the same response, which basically seemed linear.
Fig. 2: Sub-regions of phase-plane leading to periodic surging and to surf-riding. The two
are separated by the inset of the unstable equilibrium, located near to the crest. Given
some assumed position of the ship on the wave, the diagram shows the range of initial
velocities that lead to periodic surging (dotted line) and those (dashed line) that lead to
stationary condition [modification of a diagram that has appeared in Spyrou (1996a)].
8
Main findings for Fn=0.30 (RPS=2.0)
Possible change of character of the periodic response near to the lower surf-riding threshold
has been targeted. An intension here was to look for a possible change of the ‘overtakingwave’ periodic motion in a way manifesting the so-called “cumulative” type of broaching.
Such a type of broaching is occasionally mentioned in the literature (Conolly 1972). For a
Japanese purse-seiner, it has been reproduced also through simulation (Spyrou 1997).
However in general it has not received much attention despite the fact there have been
reports referring to it for over 40 years. The study of the purse-seiner had been based on a
simpler mathematical model, which was taking into account motions in surge, sway, yaw
and roll. Moreover, that ship had been assumed to be contouring the wave in condition of
hydrostatic balance as far as heave and pitch are concerned. No account of “hydrodynamic
memory” due to waves radiated from the ship had been considered. It remains therefore to
be verified whether such a “periodic” type of broaching is “generic”. Apparently, the
possibility to be reproduced by LAMP should add a lot of confidence about it.
The tumblehome’s centre of gravity has been set by 1.0 m higher than the originally
specified value, in order to rule out the possibility of capsize or even the interference of
nonlinear rolling. Another choice that had to be made was the setting of gain values for the
rudder’s controller. In the first instance the proportional as well as the differential gain were
given the value of 3.0, in their respective units. No integral gain was used. One should note
here that for the assumed wave-length-to ship-length-ratio, wave celerity’s Froude number
is just about 0.4. Therefore one expects the waves to be clearly overtaking the ship.
Commanded heading ψr was selected as the control parameter, varied successively up to 24
deg. For each ψr value a 20 min (real-time) simulation was performed as a minimum.
Time history and phase plots of a typical periodic response are shown in Fig. 3 for ψr=14
deg. Increase of the commanded heading ψr to 19 deg induced however substantial change
of pattern (Fig.4). As revealed from the time history of surge velocity, the ship has transited
from the ordinary response to one that is definitely different in a qualitative sense. The
following characteristics are singled out from this transition:
A gradual drop of the mean speed was realised and eventually a quasi-periodic
pattern of response has emerged.
The power spectrum produced on the basis of surge velocity’s time series (steadystate part) shows clearly the existence of a second incommensurate period which is
a very long one (over 400 s!). In a practical context such a period is of course
insignificant, however what is important here is the awareness of the realised
change of state per se.
The emerging response appears like an erratic zig-zag oscillation that one might
interpret as a succession of incomplete turns. The true trajectory of the ship in global
X-Y axes is also shown.
9
Yaw was augmented, bringing about much larger rudder oscillations that reached
the max deflection (assumed to be 35 deg). Although the transition might not be
labelled as too abrupt, it definitely manifests the occurrence of a jump phenomenon
and it should be classified as occurrence of broaching.
Heave and pitch amplitudes have grown very considerably. In fact, the possibility of
considerable transfer of energy into these modes is something observed for the first
time due to the inclusion of all six degrees of freedom in the mathematical model.
Hence, the observed significant drop of mean speed could be attributed to the
emerging severe horizontal and vertical oscillations.
To understand the sensitivity of the observed phenomena to control selections, the above
investigation was repeated for a higher proportional gain (ψ = 5.0). In Fig. 5 are shown key
features of the observed system dynamics, for ψr =15 deg. The frequency of response and
the pattern of surging indicate that the system has undergone a period-doubling i.e. the
response has turned into sub-harmonic. This is verified by the power spectrum of the time
record of the steady surge motion. This qualitative change should have been invoked by a so
called “flip” bifurcation. Simultaneously, substantial increase in the amplitude of the
performed yaw/rudder oscillation is realised. Indeed, a consequence of the flip bifurcation is
known to be that, the amplitude of the generated sub-harmonic response is increasing fast
as the control parameter departs from its flip bifurcation point value.
Slight increase of the commanded heading, to ψr =16 deg, gave rise to the earlier discussed
jump phenomenon towards the quasi-periodic response pattern (Fig. 6). However, in this
case the yaw-rudder oscillation had already become large before realisation of the jump, as
consequence of period doubling. It is conjectured that this occurrence is a feature that
resides primarily in the horizontal plane dynamics of the ship. The realised drop of mean
speed is in this case associated mostly with transfer of energy into heave and pitch with
subsequent growth of these vertical plane motions. Obviously the higher gain value
restricted the yaw/sway oscillations (adding to this, sway velocity did not show any
significant increase before and after the jump).
In Fig. 7 is summarised the variation of steady yaw amplitude as function of ψr, for Fn=0.30.
One should note some substantial increase in this amplitude, once the surge motion pattern
has turned subharmonic. The jump towards an erratic quasi-periodic pattern of motion,
which happens at higher ψr, can also be observed in this diagram. Rather counterintuitively
however, it comes about as slight reduction of the already very large yaw oscillation
amplitude. On the other hand, this jump phenomenon is characterised by substantial energy
transfer into vertical motions (heave and pitch). Surge time histories just before and right
after the occurrence of period-doubling are shown in Fig. 8.
10
Fig. 3: Fn=0.30 (RPS=2.0), ψr=14 deg, aψ=3
pitch (deg)
2
1
0
1
2
0.9
1
1.1
1.2
1.3
heave (m)
11
1.4
1.5
1.6
1.7
Fig. 4: Fn=0.30 (RPS=2.0), ψr=19 deg, aψ=3
surge velocity (m/s)
12
10
8
0
200
400
600
800
1000
1200
1400
1600
1800
time (s)
1350
1000
Sfft j fft
500
0.00019
0
0
10
20
1
30
40
jfft
45
240
200
150
Y (m)
6
2 100
P
50
0
40
300
2000
4000
6000
6
4
12
2
P5
4
1 10
1
P
X (m)
8000
4
1.2 10
4
1.4 10
4
1.6 10
17000
Fig 4 (cont.): Fn=0.30 (RPS=2.0), ψr=19 deg, aψ=3
18
10
rudder angle (deg)
R1
0
0
10
R1
10
20
20
30
30
40 40
40
20
20
15
20 15
10
4.5
5
0
P6 0
yaw angle (deg)
P6
yaw angle (deg)
10
5
12.94
5
5
10
10
10
4
pitch (deg)
2
R2
P5
0
2
35
3.81599
3.3469
4
0.5
0
0
0.46377
0.5
1
1.5P 2
P3
heave (m)
2
2.5
3
3.5
119.7
3.31663
4
3
heave (m)
rudder angle (deg)
10
2
3
P
1
0
0.46806
1
400
600
800
0
P
time (s)
300
13
1000
1200
1400
1400
Fig. 5: Fn=0.30 (RPS=2.0), ψr=15 deg, aψ=5
12
11.5186
surge velocity (m/s )
11
1
V
10
9
8.67294
8
0
100
200
0
12
300
0
( R)
time (s)
400
500
600
600
1300
surge velocity (m/s )
rudder angle (deg)
11
10
R1
10
1000
1
V0
Sfft j fft
10
500
9
20
7
1.000510
30
0
0
8
0
10100
200
20
1
300
30
0
( Rjfft
)
time (s)
40040
500 50
600 60
60
30
40
20
15
rudder angle (deg)
20
10
5
0
5
10
P6
yaw angle (deg)
0
R1
20
40 40
10
5
0
14
P6
yaw angle (deg)
3
2
1
14
5
10
15
15
Fig. 5 (cont.): Fn=0.30 (RPS=2.0), ψr=15 deg, aψ=5
2
2
heave (m/s )
1.5
3
P
1
0.5 0.5
100
200
300
400
500
600
0
( R)
time (s)
40
600
Fig. 6: Fn=0.30 (RPS=2.0), ψr =16 deg, aψ=5
3.4
20
rudder angle (deg)
heave (m)
3
0
2
R1
3
P
201
0.4
0
40
0
10
200
5
400
600
P6
0
yaw
( R) angle (deg)
time (s)
0
11.3827
0
5
800
10
15
900
12
(m/s )
velocity
surgeangle
(deg)
rudder
20
11
10
0
1
V
R1
9
20 8
7.23992
7
0
40 0
200
10
400
600
0
( R) 0
5
time (s)
P6
yaw angle (deg)
5
15
g)
20
800
10
900
15
Fig. 7: Summary of behaviour for Fn=0.30 (RPS=2.0), ψr varied from 4 to 20 deg and aψ=5
Fig. 8: Simulations just before (upper diagram) and just after (lower) the change towards
subharmonic motion (Fn=0.30, ψr=8 deg, aψ=5)
12.383
12.5
surge velocity (m/s )
12
11.5
V1
11
10.5
10
9.788 9.5
460
480
500
520
450
12.7
540
560
580
V0
time (s)
600
599.8
12.5
surge velocity (m/s )
12
11.5
V1
11
10.5
10
9.762 9.5
1060
1080
1100
1120
1050
V0
time (s)
16
1140
1160
1180
1200
3
1.210
Findings for Fn=0.33 (RPS=2.2)
At this stage propeller’s revolutions were increased in order to make nominal speed exceed
the threshold RPS value for surf-riding. The term “global” implies here that, for an exactly
following seaway, the ship should be attracted towards surf-riding from any initial condition
of surge. More generally of course, the type of response adopted should also be dependent
on the direction of motion; which is controlled through the setting of commanded heading.
In this particular study the commanded heading was varied in a stepped manner, within the
range between 0 and 20 deg. In the first batch of simulations, the initial surge velocity was
set at a very low value in order to ensure that we don’t miss out on any coexisting periodic
motions. Summary of the obtained responses, characterised by their steady amplitude, is
shown in Fig. 9. At ψr just above 10 deg, the “edge” of the surf-riding domain is reached.
The escape from surf-riding that inevitably occurs once this boundary is crossed (resulting
simply from the fact that there is not sufficient magnitude of surge wave force to fill the
difference between effective thrust and resistance at a speed equal to c/cosψ) is essentially
a jump phenomenon, leading abruptly back to the “overtaking-waves” periodic motions.
Energy wise, this corresponds to a jump from a high potential energy level to a lower one
where the difference is the higher kinetic energy in the latter condition. Obviously, the
higher the difference of potential energy levels the higher the kinetic energy of the periodic
response.
The observed evolution of periodic responses after the escape from surf-riding are
described next: initially, these appear as ordinary asymmetric (nonlinear) responses at the
frequency of encounter. By further increase however of the commanded heading they turn
into subharmonic according to the pattern discussed earlier for the lower nominal speed. A
typical time history of subharmonic steady surging can be seen in Fig. 10, exhibiting also the
characteristic longer stay in the vicinity of wave crests. It is remarkable that, further increase
of ψr invokes return to the ordinary periodic response pattern which persists for
commanded headings at least up to ψr=20 deg. Several of the above runs have been
repeated from a high initial surge velocity, in order to capture the range of coexistence of
surf-riding and periodic responses.
The above procedure was repeated for higher proportional gain (αψ=5.0). A very interesting
observation was that, if the initial velocity had been set at a very low value, the higher gain
produced a periodic surge response; whereas the lower gain led always to surf-riding. A
comparison of initial transients for the two discussed cases of control, as observed in the
time domain for the various state variables of the system, is shown in Fig. 11 (ψr = 4 deg).
Considering especially the time histories of surge and yaw, it is quite intriguing that despite
the very little difference at the initial stage, qualitatively different responses were
eventually adopted.
17
Fig. 9: Fn=0.33, ψr from 0 to 20 deg, αψ = 3
surf-riding only
coexistence of surf-riding with ordinary periodic surging
coexistence of surf-riding with subharmonic surging
coexistence of surf-riding with ordinary periodic surging
ordinary periodic only
Fig. 10: Fn=0.33, ψr =12 deg, αψ = 3:
subharmonic surging with apparent longer stay near crests
14.2 14
surge velocity (m/s )
I:
II:
III:
IV:
V:
1
V
12
10 10
840
860
880
900
920
0
R
time (s)
840
18
940
960
980
1000
1000
Fig. 11: The proportional gain influences whether the ship remains on surf-riding or
returns to periodic response. To be noted that the rudder ‘hits’ the max as the ship is on
its way to either of patterns:
a) Fn=0.33, ψr=4 deg, αψ=3 -> surf-riding
R1 ( R)
j
N N1 j 1 1
30
surge velocity (m/s )
yaw angle (m/s)
20
0
P6
20
20
040
6
P
60 2060
0
50
100
150
00
40
60
40
0
200
250
P0
time (s)
50
100
150
0
P
time (s)
300
300
200
250
300
40
( R)
1
0
20
35
0
50
100
150
0
P
time (s)
0
2.6
pitch angle (deg)
rudder angle (deg)
20
P5
19
1.9
200
250
300
300
rudder an
1
( R)
0
20
Fig. 11 (cont.)
0
50
100
150
0
P
time (s)
200
250
300
17
surge velocity (m/s )
15
10
1
V
15
surge velocity (m/s )
5
0
10
0 1
V 0
0
50
100
150
0
P
time (s)
250
300
300
5
9.68
0
0
50
100
R1
R2
200
40
150
0
P
time (s)
200
250
300
40
roll angle (deg)
20
11.31
4
0
P
3.097
P 1 P 2
1.595
20
40 40
0
50
100
150
0
P
time (s)
0
20
200
250
300
300
Fig. 11 (cont.)
b) Fn=0.33 (RPS=2.2), ψr=4 deg, αψ=5 -> periodic surging
R1 ( R)
N N1 j 1 1
j
30
yaw angle (m/s)
20
0
6
P
20
40
60 60
0
50
100
150
0
40
200
250
300
P0
time (s)
300
40
1
R
0
20
35
0
50
100
150
0
P
time (s)
0
200
250
300
300
2.6
2.4
pitch angle (deg)
rudder angle (deg)
20
P5
2.2
21
2
8
6
4
2
0
2
4
6
8
rudder angle (deg)
1
R
0
20
Fig. 11 (cont.)
0
100
150
0
P
time (s)
200
250
300
15
10
15
1
V
surge velocity (m/s )
surge velocity (m/s )
16
50
5
10
1
V
0
0
0
0 5
50
100
150
0
V
time (s)
200
250
300
300
9.68
0
32.1242
R1
0
50
100
150
0
V
time (s)
200
250
300
40
R2
roll angle (deg)
20
4
P
0
11.31
3.097
20
32.1385 40
P 1 P 2
0
50
100
1.595
150
0
P
time (s)
0
22
200
250
300
300
Fn=0.36 (RPS=2.4)
A remarkable new feature that appeared at this higher propeller rate was a stable
oscillatory type of surf-riding, residing at the outskirts (in terms of commanded heading) of
the domain of surf-riding. As the ship is carried along by a single wave, it is also oscillating in
all directions on the wave’s down-slope. Through stability analysis of the stationary surfriding states, this quite fascinating occurrence has been explained in the past to be the
result of a Hopf bifurcation (Spyrou 1996a). Moreover, it was indicated that the differential
gain plays a very influential role, controlling the onset of such behavior. In Fig. 12 is shown
the inception and further metamorphosis of surge velocity’s oscillation, by varying gradually
the commanded heading from lower to higher values. The oscillations emerge as a stable
pattern of response at ψr about 9 deg. The frequency content of the surge response at some
characteristic commanded heading angle can be deduced from the power spectrum of Fig.
13 (some distortion of the time history of surge velocity from the harmonic pattern is
noticed, due to the presence of higher harmonics, most prevalently the 2ω harmonic).
Interesting is also the coexistence of oscillatory surf-riding with the more ordinary
overtaking wave periodic motion, for a range of commanded headings (Fig. 14).
Fn=0.38 (RPS=2.6) and Fn=0.41 (RPS=2.8)
The two finally tested speeds were selected slightly below and slightly above wave celerity
(for a wave-length-to-ship-length-ratio equal to 1.0 this corresponds to approximately
Fn=0.4). A common characteristic is that the character of oscillatory surf-riding becomes
more complex and in some case period-doubled oscillations are noticed (Fig. 15). These
matters as well as the disappearance of surf-riding oscillations have been discussed in
Spyrou (1996b). A characteristic transition from oscillatory surf-riding towards the ordinary
overtaking wave response can be seen in Fig. 16; at such a large commanded heading value
oscillatory surf-riding cannot be sustained and the ship inevitably escapes.
Further parametric study
The above investigation was repeated for extremely different values (much higher) of the
yaw and pitch radii of gyration. However the above described pattern was found to persist
and thus it is believed that it has quite a generic character.
23
Summary of findings of simulation studies
The range of control parameters’ values (nominal Froude number and commanded
heading) where surf-riding appears according to the above investigation is shown in Fig.
17.
The phenomena encountered comply qualitatively very well with the expectations
according to the current theory of broaching.
Critical conditions of surf-riding and broaching that could be the target of probability
calculations can be rationally specified.
24
Fig. 12: Gradual emergence of oscillatory surf-riding; Fn=0.36, αψ=3
a) ψr=4 deg
surge velocity (m/s )
18
18
16
V1
14
11.9 12
50
100
150
30
200
250
V0
time (s)
300
300
b) ψr=8 deg
surge velocity (m/s )
18
18
16
V1
14
12
11
50
100
150
200
30
V0
time (s)
25
250
300
350
400
400
Fig.12 (cont)
c) ψr=12 deg
surge velocity (m/s )
18
18
16
V1
14
12
10.5
100
200
30
300
400
V0
time (s)
500
500
d) ψr=16 deg: capture into oscillatory surf-riding after long transient
18
18
surge velocity (m/s )
16
14
V1
12
10
9
200
400
30
600
V0
time (s)
26
800
1000
1000
Fig. 13: Fn=0.36, Ψr=17.9 deg, αΨ=3:
a) Typical oscillatory surf-riding with presence of higher harmonic apparent.
17
surge velocity (m/s )
17
16.5
V1
16
15.5 15.5
620
640
660
610
680
700
V0
time (s)
703
b) Power spectrum of above time series
270
200
Sfft j fft
100
0
0
0
10
20
30
1
40
jfft
27
50
60
70
70
Fig. 13 (cont.)
c) Phase-plane projections of steady-state orbits.
30
30
rudder angle (deg)
20
30
10
R1
20
rudder angle (deg)
0
10
R1
10
15
6
0
4
2
0
5.26544
2
4
6
8
P6
yaw angle (deg)
7.19312
19.86
10
6
4
2
0
2
4
6
8
P6
yaw angle (deg)
2.6
2.6
R2
pitch angle (deg)
2.4
P5
2.2
14.09
0
P2
119.7
2
1.9
8
8
6
4
2
0
2
4
6
P4
roll angle (deg)
8
9.68
R1
R2
28
11.31
3.097
P 1 P 2
8
1.595
Fig. 14: Fn=0.36, ψr=20 deg, αψ=3; coexistence of oscillatory surf-riding with ordinary
“overtaking wave” periodic motion. The adopted pattern depends on initial conditions
a) initiated from very high initial speed -> oscillatory surf-riding
R1 ( R)
j
surge velocity (m/s )
18.5
N N1 j 1 1
18
16
1
V
14
12 12
0
100
200
0
V
time (s)
000
300
400
400
30
b) initiated from very low initial speed -> periodic surging after long transient
surge velocity (m/s ) rudder angle (deg)
20
R1 ( R)
j
16.9
R1
N N1 j 1 1
10
15
0
10
1
V
10
7.06
7.062
7.064
7.066
7.068
7.07
7.072
7.074
P6
yaw angle (deg)
5
2
0
200
400
600
800
0
V
time (s)
0
30
rudder angle (deg)
20
10
R1
0
29
1000
1200
1400
1400
7.076
Fig. 15: Fn =0.41, ψr=24, aψ=3: period-doubled oscillatory surf-riding
R1 ( R)
j
N N1 j 1 1
18.5
surge velocity (m/s )
18
1
V
17
16
15.5
400
600
800
0
V
time (s)
400
1000
1200
1200
30
20
surge velocity (m/s )
rudder angle (deg)
Fig. 16: Fn =0.41, ψr=24, aψ=3: transition from oscillatory surf-riding to ordinary periodic
motion
10
R1
R1 ( R)
j
20
20
0
N N1 j 1 1
18
10
1
V
16 5.36
5.358
5.356
5.354
5.352
5.35
5.348
5.346
P6
yaw angle (deg)
14
12
11
0
200
400
0
V
time (s)
000
30
30
(deg)
20
600
800
900
Fig. 17: Boudaries of surf-riding on the plane of control parameters. In the diagram are
shown: the domain of surf-riding from some initial conditions (filled circles); the domain
of surf-riding from all initial conditions (empty circles); and the boundary of oscillatory
surf-riding (triangles).
31
6. INTRODUCTION TO CONTINUATION METHODS
For nonlinear and multi-dimensional systems, simulation alone is usually ineffective. For
such systems, behavior can depend critically on initial conditions; and there are too many of
these to be examined. Quick and reliable elicitation of key dynamical information through
implementation of more advanced numerical techniques is thus quite essential. An
enhanced exploration potential is offered by the so-called “Continuation” methods. The
main utility of Continuation is for identifying how the steady-states of a dynamical system
move in state-space when one or more control parameters are varied. Unstable solutions
can be traced by a Continuation scheme without problems. Steady-states represent of
course only the long-term behavior of a dynamical system. Nevertheless, knowing the type
and arrangement of all stable and unstable steady states is essential for effective
performance assessment of a dynamical system because steady-states dictate state-space’s
topology. In general, Continuation currently could work for stationary as well as for periodic
states. Continuation can be used also for tracing the locus of bifurcation points as some
parameters are varied.
A mathematical model brought into the canonical form of a system of ordinary differential
equations (possibly combined with algebraic equations) represents the usual input.
Mathematical details about the principle of Continuation for such systems can be found in
Spyrou et al. (2007). Efforts are noted also towards implementing continuation schemes for
dynamical systems described by partial differential equations (Shroff & Keller, 1993). In
Continuation the length of the solution curve is often used as a parameter of the numerical
scheme in order to pass over bifurcation points, where for example, the solution curve
might fold back or continue in more than one direction. This technique is referred-to in the
literature as ‘pseudo-arc-length’ continuation. The concept seems to have been introduced
by Keller (1978) and it was soon taken up by several people, initially in the dynamics
research community, like Kubicek & Marek (1982), Rheinboldt (1983) and others. An
accessible introduction to the mathematics involved in popular continuation analysis
algorithms can be found, for example, in Seydel (1994).
A few packages currently exist for conducting Continuation analysis. DsTool (Back et al.,
1992) from Cornell University for example, is a user friendly dynamics investigation tool with
some limited continuation capability. It worked initially in UNIX but currently is operational
on several platforms. It also incorporates a simple Graphic User Interface (GUI) capability.
The package LOCBIF can perform Continuation of equilibrium and periodic solutions of lowdimensional (below order 10) dynamical systems (Khibnik et al., 1993). CONTENT is designed
to perform simulation, Continuation, and normal form analysis of dynamical systems
(Kuznetsov & Levitin, 1995-1997). It is able to predict bifurcations of equilibria. It operates in
UNIX, Windows, Linux and other environments. The most computationally advanced
algorithm until recently was AUTO 86/87 (Doedel et al., 1997). It can be used for performing
Continuation of equilibrium, periodic and homoclinic solutions. Also, it can locate most
types of bifurcation points. It operates in UNIX and it was initially written in Fortran. AUTO
32
2000 is a more recent version written in C programming language. Unfortunately, the
interface is not particularly suitable for a non-specialist and this imposes difficulties when
analyzing systems showing complex behavior. A more friendly version of this is MatCont
that exploits MATLAB’S powerful capabilities including its graphics inerface (Dhooge et al.,
2003). MatCont is capable to trace equilibria and identify associated bifurcations such as
saddle-node, pitchfork and Hopf points. Once such singular points of the system are
detected, there is possibility to perform the so-called “codimension-2” continuation in order
to discover bifurcation points that could arise when two control parameters of the system
are altered simultaneously. This does not mean of course that the parameters are varied
simultaneously in a physical sense, but rather that regions of system’s parameter space are
associated with specific types of behavior, separated by the obtained continuation curves.
Continuation of homoclinic orbits has recently being introduced in MatCont by using
HOMCONT (Kuznetsov & Levitin, 1995-1997) toolbox which was formerly developed in
AUTO (Doedel et al., 1997).
7. MATHEMATICAL BASIS OF THE APPLIED CONTINUATION SCHEME
Continuation has been applied in the current project in order to identify the states of surfriding in quartering seas and for checking their stability properties through calculation of
their eigenvalues and eigenvectors. Consider for example our dynamical system of ship
motion is quartering seas, presented in generic o.d.e. form:
x f x; α;t
(1)
where x is the state vector of dimension say n; α is the control parameters’ vector, say of
dimension m, which could include real controls like the rudder angle and the commanded
heading, parameters that represent the environment, like the wave height, length etc. or
even ship design parameters. As usual, t is the time. Strictly speaking, bringing our ship
motion equations to the form of equation (1) may not be possible without having to deal
essentially with an infinite number of ordinary differential equations. However, often a
reasonable low dimensional approximation is possible to be specified. Stationary states can
be identified if one could remove from (1) the explicit time dependence and then request all
components xi of the velocity vector to be zero. Then (1) would become:
f x; α 0
(2)
Removal of the direct time dependence is possible for the pressure-related terms for which
ship’s position in waves, rather than time, is the key factor in the calculation; but it may not
be trivial for perturbation/radiation terms which are in fact responsible for the infinite
dimensional nature of the ship motions problem. The applicability of Continuation seems
thus to be dependent on how important these terms are, for the specific scenarios that are
33
investigated; and if they are, whether some reasonable low-dimensional approximations of
these could be produced.
Differentiation of (2) with respect to some selected control variable αk leads to having to
solve equivalently the following system of algebraic equations:
J 1 x;
dx
f
,
d
x 0 x0
(3)
The scheme can be applied for all control variables; but for the sake of simplicity let’s
assume here that we are dealing with only one (at a time) varying control .
Unfortunately, the Jacobian J is singular at bifurcation points and once we approach such
points the scheme would break down. Therefore, it would be impossible to go over
bifurcation points, especially over foldings of the steady-states locus (Fig. 17). To circumvent
this, the length z of the solution curve is then introduced as an additional parameter – by
this, from here on one follows the classical procedure known in the field as “arc-length
continuation” or its variant “pseudo-arc length continuation” [e.g. Keller (1978); Kubicek &
Marek (1982)]:
dfi
dz
n
j 1
fi dx j fi d
0
x j dz dz
(4)
Obviously an additional equation should be introduced giving the arc length, based on Pythagoras’
theorem:
n 1
2
dx j
1
dz
j 1
(5)
For convenience, one could name αk as the xn+1 variable. Then (4) and (5) can be brought
into the more general form of an algebraic system of n linear, plus 1 nonlinear, equations.
dxi
More specifically, from (4) the terms
can be expressed with respect to a single derivative of
dz
dx
some arbitrarily selected state variable k :
dz
dxi
dx
i k , i 1, 2,..., k 1, k 1,..., n 1
dz
dz
(6)
The coefficients i , i 1, 2,..., k 1, k 1,..., n 1 are identified by applying Gauss elimination on
the n x n Jacobian matrix:
34
f1 f1
f1
f1
x ,..., x , x ,..., x
k 1
k 1
n 1
1
J k ...
f n
f
f
f
,..., n , n ,..., n
xk 1 xk 1
xn 1
x1
Then
(7)
dxk
can be calculated from the transformed version of (5) after inserting equation (6) into
dz
(5):
dxk
dz 1
2
n 1
i 1
ik
2
i
1
The sign of the derivative
(8)
dxk
is specified by the way one moves along the solution curve (to the
dz
right or to the left).
Automated stability analysis at each state should accompany continuation. This is done in
the standard way by means of local linearisation, calculation of the characteristic polynomial
from the Jacobian matrix, and from the polynomial’s roots derivation of eigenvalues and
eigenvectors (see for example Strogatz 2001).
8. LAMP-BASED IMPLEMENTATION OF CONTINUATION
Continuation may play several roles in the development and application of the procedure of
capsizing probability calculation, considering surf-riding and broaching. First and foremost,
the dynamic equilibria of the ship in following or quartering seas describe not only stable
solutions of surf-riding and broaching but also the unstable solutions that define the
boundaries of transition between regions of regular oscillatory motion in waves and surfriding and/or broaching. The analysis will eventually be used to define the boundary
between the “non-rare” and “rare” problems associated with these phenomena and help to
develop the form of the “rare” problem that must solved and the upcrossing states that will
define the specific “rare” problem for each upcrossing. To some extent, the results of the
Continuation analysis plays a roll with regard to broaching and surf-riding that the roll
restoring (GZ) curve plays in capsizing itself.
In addition, Continuation’s ability to look into details of the ship’s dynamic model allows it
to play a key part in the evaluation of the ship motion simulation system’s capability to
properly capture the physics of the phenomena and the importance of different models or
terms of the simulation on the evaluation of dynamic stability. This ability will be used for
35
the analysis, verification, and development of the LAMP-based force modeling itself as well
as the analysis, verification, and development of simulation model input for specific ship
configurations.
The present formulation of the Continuation analysis for equilibria is based on a dynamic
system defined by a set of ordinary differential equations (ODEs) of the form of equation
(1). In our case, x is a state variable vector consisting of the ship position and velocity and
is a control vector consisting of such quantities as propeller speed, command heading, etc.
The right hand side (RHS) vector f is a function of the external forces acting on the ship, the
mass and mass moments of inertia, the kinematic relationships between the derivatives of
the position vector and velocity, and (if necessary) control or servo models. Continuation
searches for the locus of solutions x and that satisfy the relationship expressed by
equation (2). These points represent solutions where the time derivatives of the state
variables are zero, so they are equilibria, although not necessarily stable ones. These
equilibria are found by performing an initial search for a single equilibrium point and then
tracking the curves of equilibria through the solution domain. The name Continuation
derives from the “continuous” nature of this tracking and the locus of solution points. The
key quantity in this search is the Jacobian matrix J, which represents the derivatives f with
f (x; α),
f ( x; α) . One of the key
respect to both the state and control vectors,
xi
a j
aspects of the Continuation formulation is the that RHS vector f and the Jacobian J must be
explicitly computed from the state vector x, control setting , and data such as the sea
conditions, hull geometry, etc. The analysis does not allow for a “memory” effect where the
forces acting on the ship are dependent on the history of how the ship got to its current
position.
Previous implementations of the Continuation approach have been primarily been based on
ODE descriptions of the motion of a ship in a seaway in which the RHS has been explicitly
calculated using terms involving powers of the state and control variables and fixed
coefficients.
One of the principal objectives of the present effort is to integrate
Continuation with LAMP’s time-domain simulation of a ship moving in waves. This was
done by integrating NTUA’s Continuation code with a subroutine-based version of LAMP’s
force calculation routines, under a new top level supervisory code. The resulting code was
named LAMPCont.
Three key parts of creating LAMPCont were:
re-implementing LAMP’s EOM in diagonalized form suitable for Continuation;
evaluating the RHS vector f using the LAMP force calculation;
evaluating the Jacobian matrix J.
36
The first of these parts was based on LAMP’s standard 6-DOF formulation of the equations
of motion, which is summarized by this equation for the LAMP User’s Manual:
(9)
Refer to Section 8 of Volume XII: LAMP Theory and Implementation for details. For use with
Continuation, a specialized set of equations was derived for each equilibrium problem that
was examined. The specialization was required in order to properly define f(x; )=0 as the
desired equilibrium condition. For example, the surf-riding problems had to be formulated
so that the state variables associated with the direction of wave propagation were defined
as the position wrt the wave crest and the ship velocity wrt the wave celerity. For cases
with autopilot and rudder servo modeling, the rudder deflection R - which LAMP computes
using a simple explicit expression decoupled from the EOM - had to be implemented as part
of the EOM using the equation (Ψ is the real heading, ΨC is the commanded heading, r is
the yaw rate and the rest are constants):
TE R R K R ( C ) K RTD r
(10)
The second part of the implementation was to evaluate the RHS vector f using the LAMP
force calculation. This involved three steps:
1. Transform the Continuation state vector x to the LAMP state variables PMGG(6) and
VMG(6), and transform the control vector to the LAMP variables RUDDEF,
PROP_RPS, etc.
2. Evaluate forces and moments acting on the ship using a subroutine implementation
of LAMP.
3. Calculate f from the EOM using the force returned by the LAMP subroutine, the
mass and mass moment of inertia, and an estimate for the added mass.
The subroutine-based implementation of LAMP uses exactly the same subroutines as a
regular LAMP simulation at each time step. This force calculation includes the bodynonlinear hydrostatics and Froude-Krylov pressure forces; appendage forces due to rudders,
bulge keels, and skegs; propeller force; maneuvering forces such as hull lift; and
approximate viscous force models. Because the force must be computed explicitly in terms
of the state and control variables, this force does not include forces associated with the
hydrodynamic perturbation potential of the wave-body interaction problem such as
radiation or diffraction. As a result, the implementation corresponds to the LAMP-0 or
“hydrostatics-only” model.
37
A third key part of the implementation is the calculation of the Jacobian matrix J. In the
present implementation, this is done by setting up “variant” state vectors with a small
perturbation of each state or control variable in turn, computing the RHS vector for the
variant states, and estimating the derivatives via finite difference:
xV x xi
f (x; α)
xi
f (xV ; α) f (x; α)
xi
(11)
An additional piece of the LAMPCont implementation was to integrate standard math
library routines for the Eigenvalue and Eigenvector calculation.
Figure 18: LAMPCont Program Structure
The structure of the LAMPCont program is show in Fig. 18. The input to LAMPCont includes:
LAMP input control file defining the geometry, appendages, and other LAMP data.
Continuation problem to run.
Initial values, upper and lower limits, and perturbation increments for a Jacobian
calculation for each state and control variable.
38
Additional problem-dependent definitions such as wave height and length, and
settings for control parameters that are not being varied in the Continuation
analysis.
Continuation controls such as step size limits and number of points.
LAMPCont’s principal output consists of the locus of equilibria solution points and the
Eigenvalues and Eigenvectors of the Jacobian at the solution points.
The following problems have been implemented in LAMPCont:
Surf-riding in following seas, 3-DOF (surge, heave, pitch) vs. propeller RPS.
Turn in calm water, 3-DOF (surge, heave, yaw) vs. rudder angle.
Turn in calm water, 6-DOF vs. rudder angle.
Surf-riding in following/quartering seas, 6-DOF vs. rudder angle (unsteered).
Surf-riding in following/quartering seas, 6-DOF vs. commanded heading (autopilot).
Surf-riding in following/quartering seas, 6-DOF vs. commanded heading (autopilot +
servo equation).
The turning problems were initially implemented primarily as successively more
complicated tests of LAMPCont, but may prove to be very useful in analyzing the
characteristics of LAMP’s “maneuvering” force models, both in an overall modeling sense
and for ship-specific model input.
As part of the LAMPCont development, a specialized graphics-tool called PltCont has been
developed to plot LAMPCont results, including projections of the solution locus and
Eigenvalue.
What about the hydrodynamics?
As mentioned above, the present formulation of the Continuation analysis requires a
“point” evaluation of f(x; ), which precludes LAMP’s regular time-domain evaluation of the
wave-body hydrodynamic disturbance, including the effects of radiation, diffraction, etc. As
a result, this initial implementation corresponds to LAMP-0, or a “hydrostatics-only”
modeling, with relatively minimal accounting for additional hydrodynamic effects, such as a
constant added mass coefficient. The impact of this issue cannot yet be evaluated and will
be a principal objective as the project continues. However, some points have been
identified.
First, even with the LAMP-0 approximation, the LAMP-based implementation of the surfriding and broaching problem should provide an adequate characterization of the problem
39
to initiate the application of the split-time approach to these phenomena, in a similar way
that the piecewise linear (PWL) model of the roll equation was the basis of the development
of the split-time approach for capsizing in beam seas. This development will include the
definition of the upcrossing scenarios related to surf-riding and broaching, the
characterization of the conditions at upcrossing, and the formulation of the rare problem(s)
for these phenomena.
Within a LAMP-3 simulation, the current LAMP-0 implementation may well be adequate for
identifying the “distance” of the ship to the equilibrium boundaries, although this is
primarily conjecture at this point. However, the hydrodynamic effects related to radiation
should be small at the point of equilibrium, since the ship motion must match the motion of
the wave. Diffraction and forward speed effects may well be important, but can, in
principle, be estimate for the a) values from the current simulation (i.e. point from which we
are measuring the distance) or b) a pre-computed approximate hydrodynamic solution
based on the ship’s traveling at an equilibrium point.
In addition, the LAMP-0 implementation of LAMPCont is likely to be more than adequate for
assessing force models, especially those related to propulsion and maneuvering.
If the disturbance hydrodynamics do seem to be important in this problem, several
approaches can be investigated to better incorporate them, including:
IRF-like expressions for disturbance potentials (especially diffraction), “linearized”
about u = Vwave .
Variation to instantaneous time-domain solution (for distance to boundary).
Expanded Continuation scheme with memory effect using some kind of state space
approximation, developed from theoretical models or characterization of LAMP
simulation results.
Expanded ODE-like terms with constant or state-dependent coefficients.
9. SAMPLE LAMPCONT RESULTS
As described above, the LAMP-based Continuation analysis has now been implemented for
several surf-riding and ship turning problems. Preliminary testing and evaluation of
LAMPCont appears to show that results are:
self-consistent, in that the equilibrium solutions, or at least stable equilibrium
solutions, can be reproduced via direct simulation;
reasonable, in that they qualitatively match the results of previous Continuation
analysis of theoretical models.
40
A few sample results are shown here.
3-DOF surf-riding in following seas
The simplest problem that has been implemented in LAMPCont is surf-riding in following
seas, solving for an equilibrium in surge, heave, and pitch as a function of propeller speed
(in revolutions per second). Fig. 19 shows a sample result for the ONR tumblehome-top ship
in a wave that is 200 m long and 4 m high. In this analysis, the ship has 2 propellers which
are modeled using a specified KT curve derived from a generic B-series propeller with
properties similar to that used for U.S. Navy test of the model 5415 destroyer hull form.
The main left hand plot is a projection of the equilibria solution locus showing the position
of the ship’s center of gravity relative to the wave crest. In this plot, X=0 and 200
corresponds the wave crest, X=100 is the wave trough, and X=50 is halfway down the face of
the wave. The curve shows that, for this wave, there is both a lower and upper limit of the
propeller speed for which a surf-riding equilibrium can be found, below which the wave
cannot supply enough extra force to propel the ship at the wave celerity and above which
the wave cannot provide enough retarding force to slow the ship to wave celerity. For
points in between, there are two equilibria.
Figure 19: LAMPCont Results for Surf-Riding in Following Seas
The two equilibria points are marked for a propeller speed of 3.0 RPS, which corresponds to
a calm water self-propulsion speed just below the wave celerity. The two equilibria points
are on the face of the wave just below the crest and just above the trough. The right hand
41
plot shows projections for the heave motion (+ down relative to calm water flotation) and
pitch (+ bow up) of the locus of equilibrium points. The inset plots on the left show the
Eignvalues at the two points, which indicate that the point near the crest is an unstable
equilibrium while the point near the trough is stable. These results seem quite reasonable
and agree well with previous theoretical analysis of the surf-riding in following seas.
As a “consistency” check of these results, a series of LAMP-0 simulations were made for the
ship operating in this wave condition with different propeller speeds. Figure 20 shows time
histories of ship velocity in the global X direction (i.e. in the direction of the wave
propagation) and heave motion (vertical position of the center of gravity wrt the mean
water surface) for two propeller speeds. In both cases, the ship starts with a speed close to
wave celerity, so a surf-riding equilibrium is likely to be reached if one exists.
Figure 20: Simulation Check for Surf-Riding in Following Seas
At the lower propeller speed (2.0 RPS), for which the Continuation analysis indicated there
was no stable surf-riding equilibrium, the ship slows to oscillatory surging a speed will below
the wave celerity. At the higher propeller speed (3.0 RPS), the ship quickly transitions to
surf-riding at a point near the wave trough, with steady heave and pitch values that match
those predicted by LAMPCont.
6-DOF surf-riding in stern oblique seas
Three different versions of the problem of 6-DOF surf-riding in long-crested following or
stern oblique seas have been implemented:
Unsteered ship with specified rudder deflection.
Steered ship with specified commanded heading and PD autopilot.
Steered ship with specified command heading, PD autopilot and rudder servo.
42
In the LAMPCont implementation of the 6-DOF surf-riding problems, the displacement and
velocities of all 6 rigid body motions are solved for except for ship’s position along the crest
of the wave (global Y coordinate), which can have a non-zero but constant velocity at
equilibrium. This unconstrained Y motion is theoretically reasonable but does prevent a
challenge to using simulation to check the predicted equilibria.
For the unsteered ship, the rudder deflection angle is the control parameter in the
Continuation calculation. For the steered ship, the command heading angle is the control
parameter and PD rudder control is applied. Two steered ship models have been
implemented. In the first, the rudder deflection is explicitly set based on the heading angle
and yaw rate via the equation:
R K R ( C ) K RTD r
(12)
This model assumes instantaneous rudder response. In the second, a servo lag term in
introduced and the deflection of the rudder is now computed using a servo equation that is
solved simultaneously with the equations of motion [see equation (10) presented earlier].
Note that the rudder deflection must be constant at the equilibrium points (by the
definition of equilibrium, its derivative must be zero like the state variables), so the
introduction of the servo equation should not change the computed equilibria. It can,
however, change the force derivatives and hence the stability of the equilibria, which can be
seen by the Eigenvalues of the Jacobian matrix. In a similar fashion, the proportional gain
of the PD rudder control, KRTD in the equations above, will not affect the equilibria but will
affect the stability of the equilibrium points.
In the present effort, the 6-DOF surf-riding Continuation analysis was performed in stern
seas for the ONR tumblehome hull form with varying degrees of directional stability. The
stability of the ship was varied by changing the longitudinal position of hull lift (effectively
changing the turning moment due to hull lift) and the size of the bilge keels, skeg, and
rudders on the computational model.
Fig. 21 shows sample results for the stability variants with a fixed propeller speed of 3.0 RPS
in a wave with that is 200 m long and 4 m high. The upper plot shows the ship’s yaw angle
relative to the wave direction (in radians) vs. the rudder deflection angle (in degrees). A
yaw angle of 0.0 corresponds to pure following seas. The lower plot shows the position of
the ship on the wave as the distance from the wave crest to the CG. As the ship gets less
and less stable, a larger rudder deflection is required to maintain a large angle relative to
the wave.
While a few self-consistency check have been performed for this case, the primary check of
these results have been compare them with the theoretical evaluations described in
previous Continuation work. Fig. 22 shows the yaw vs. rudder equilibrium solution loci for
two of the stability variations along with a drawing describing the stability characteristics of
43
a similar plot derived from a theoretical analysis of the problem. The results are
qualitatively very similar including the presence of a region of Hopf bifurcation, which can
be identified by the position of the Eigenvalues in the inset plots of the LAMPCont results.
Figure 23 show some results for a similar set of calculations for a larger of wave of 6m
height. The left hand plot shows the expected result that a larger rudder deflection is now
required to maintain an oblique angle to the wave, but also show an incomplete solution
locus. The right hand plot shows the solution locus plot vs. the point index, which
effectively is the “path” that Continuation takes to track the solution. It shows points where
the Continuation solution gets stuck or doubles back on itself. This is a technical issue in the
Continuation code that will need to be resolved in future work.
Figure 21: LAMPCont Results For Surf-Riding in Quartering Seas (Unsteered)
44
Figure 22: LAMPCont Comparison to Theoretical Analysis of Surf-Riding for an Unsteered
Ship
Figure 23: LAMPCont Results for Unsteered Ship in 6m Wave
45
Fig. 24 shows initial results for ship with autopilot in the 4m wave, again at a propeller
speed of 3.0 RPS. These results show that the surf-riding equilibrium for the stable ship is
very close to the command heading, while the less stable variants are seeing a larger and
larger error in the realized vs. command course.
Figure 24: LAMPCont Results for Steered Ship
46
Fig. 25 shows similar results for the larger wave height of 6m. For the least stable variants, a
small commanded heading can induce a very large yaw angle. The smaller plot shows the
rudder angles associated with the equilibrium solutions. There is a slight “hitch” in the
curve for the most unstable case, which appears to be a result of the stall of the hull or
appendage lift due to very large sideslip angle. By default, the LAMP lift model is
discontinuous at this point, with the lift expression abruptly replaced by an empirical “eddymaking” model. Despite this discontinuity in a principal force mode, the Continuation
method was able to track the equilibrium solution locus through this region.
This calculation was made with the servo model. However, the solution locus was found to
be the same with and without the servo model, which is exactly what was expected.
Figure 25: LAMPCont Results for Steered Ship in 6m
10. LAMPCONT: SUMMARY OF STATUS AND FUTURE WORK
In summary, the Continuation analysis procedure developed at NTUA has been successfully
integrated with a subroutine-based force calculation to produce a LAMP-based Continuation
analysis code called LAMPCont. While some technical issues in the equilibria tracking
remain to be resolved, the basic implementation has been shown to be consistent with
direct simulation of the stable equilibria as well as the results of previous Continuation
research of a more theoretical nature. However, the formulation of the Continuation
approach does not allow for memory in the calculation and has therefore prevented the
inclusion of the full LAMP-based hydrodynamic force calculation. As a result, the current
47
implementation is based on the LAMP-0 or “hydrostatic-only” model of wave-body
hydrodynamics.
With this basic implementation of the Continuation approach completed, future work will
focus primarily on investigating and possibly mitigating the effects of the hydrodynamic
model approximations and the application of the Continuation approach to the
characterization of surf-riding and broaching, with the goal being to predict a probability of
surf-riding and broaching in irregular waves and an evaluation of the risk of such
phenomena leading to extreme roll motion including capsizing.
11. CONCLUSIONS
Behavior that is consistent with the existing theory of broaching and surf-riding has been
reproduced by use of LAMP for the tumblehome-top ship.
Capture into surf-riding in quartering seas, as well as escape from it, that takes into account
all six degrees of freedom of ship motion has been studied for the first time.
Continuation has been introduced into LAMP for the investigation of surf-riding in
quartering seas, with or without active control. Several example scenarios were run
successfully.
The developed approach that combines focused simulations with Continuation is suitable
for developing specifications of critical conditions for probability calculation.
48
REFERENCES
Belenky, V.L. & Sevastianov, N.B. (2007) Stability and safety of ships: risk of capsizing,
SNAME, Jersey City, ISBN 0-939773-61-9.
Bishop, B, Belknap, W., Turner, C., Simon, B., Kim, J. (2005) Parametric Investigation on the
Influence of GM, Roll Damping, and Above-Water Form on the Roll Response of Model 5613,
Report NSWCCD-50-TR-2005/027, Naval Surface Warfare Center/Carderock Division, West
Bethesda, Maryland.
Conolly, J.E. (1972) Stability and control in waves: a survey of the problem. Proceedings of
International Symposium on Directional Stability and Control of Bodies Moving In Water.
Journal of Mechanical Science 14(7), (Supplementary Issue), 186–193.
Dhooge, A, Govaerts, W, Kuznetsov, Y.A., Mestrom, W., Riet, A.M., and Sautois, B., (2003).
"MATCONT and CL_MATCONT: Continuation Toolboxes for MATLAB," Report of Gent
(Belgium) and Utrecht (Netherlands) Universities.
Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y., Sandstede, B., and Wang, X.J.
(1997). "AUTO97-00: Continuation and Bifurcation software for Ordinary Differential
Equations (with HomCont), " User’s Guide, Concordia University, Montreal, Canada.
Du Cane, P. & Goodrich G.J. (1962) The following sea, broaching and surging, Transactions of
the Royal Institution of Naval Architects, 104.
Grim, O. (1951) Das Schiff in von Achtern Auflaufender See, Jahrbuch der
Schiffbautechnischen Gesellschaft, 45.
Grim, O. (1962) Surging motion and broaching tendencies in a severe irregular sea. Report
929, Davidson Laboratory.
Grim, O . (1983) Das Schiff in von achtern kommendem Seegang. Schiffstechnik 30, 84-94.
IMO (2007) Revised Guidance to the Master for avoiding dangerous situations in adverse
weather and sea conditions, MSC.1/Circ.1228
Kan, M. (1990) Surging of large amplitude and surf-riding of ships in following seas. Selected
papers in Naval Architecture and Ocean Engineering, no. 28. Tokyo: Society of Naval
Architects of Japan.
Keller, H.B. (1977) Numerical solution of bifurcation and nonlinear eigenvalue problems,
Applications of Bifurcation Theory, P.H. Rabinowitz, ed., New York, London, Academic Press,
359-384.
Khibnik, A., Kuznetsov, Y.A., Levitin, V.V., and Nikolaev, E.V. (1993) Continuation Techniques
and Interactive Software for Bifurcation Analysis of ODEs and Iterated Maps, Physica D, 62,
No 1-4, 360-371.
49
Kubicek, M. & Marek, M. (1983) Computational methods in bifurcation theory and
dissipative structures. Springer.
Kuzntsov, Y.A. & Levitin, V.V. (1995-1997). CONTENT: A Multiplatform Environment for
Analyzing Dynamical Systems, Dynamical Systems Laboratory, CWI, Amsterdam, The
Netherlands.
Lin, W.M. Zhang, S., Weems, K. and D. Luit, (2006) Numerical Simulations of Ship Maneuvering in
Waves, Proceedings of 26th Symposium on Naval Hydrodynamics, Rome, Italy.
Makov, Y. (1969) Some results of theoretical analysis of surf-riding in following seas (in
Russian), Transactions of the Krylov Society 126, 124-128.
Motora, S., Fujino, M., Koyonagi, M., Ishida, S., Shimada, K. & Maki, T. (1981) A
consideration on the mechanism of occurrence of broaching-to phenomena. Selected
Papers in Naval Architecture and Ocean Engineering, no. 150. Tokyo: Society of Naval
Architects of Japan.
Rheinboldt, W.C. & Burkardt, J.V. (1983) A locally parameterized continuation process, ACM
Transactions on Mathematical Software, 9, 2, 215-235.
Shroff, G.M. & Keller, H.B. (1993) Stabilization of unstable procedures: the recursive
projection method. SIAM Journal of numerical analysis, 30, 4, 1099-1120.
Spyrou, K.J. (1996a) Dynamic instability in quartering seas: The behavior of a ship during
broaching, Journal of Ship Research, SNAME, 40, No 1, 46-59.
Spyrou, K.J. (1996b) Homoclinic connections and period doublings of a ship advancing in
quartering waves (1996). CHAOS: An Interdisciplinary Journal of Nonlinear Science, American
Institute of Physics (AIP), 6, No. 2, 209-218.
Spyrou, K.J. (1996c) Dynamic instability in quartering seas-Part II: Analysis of ship roll and
capsize for broaching, Journal of Ship Research, SNAME, 40, No 4, 326-336.
Spyrou, K.J. (1997) Dynamic instability in quartering seas-Part III: Nonlinear effects on
periodic motions. Journal of Ship Research, SNAME, 41, No 3, 210-223.
Strogatz, S. (2001). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology,
Chemistry, and Engineering (Studies in Nonlinearity). Perseus Books Group. ISBN 0-73820453-6.
Wahab, R. & Swaan, W. A. 1964 Coursekeeping and broaching of ships in following seas.
Journal of Ship Research, 7, 1-15.
Weinblum, G. & St. Denis, M. (1950) On the motions of ships at sea. SNAME Transactions,
58, 184-248.
50
View publication stats