Numerical prediction of the surf-riding threshold of a ship in stern quartering waves in the light of bifurcation theory

O R I G I N A L A R T I C L E

Numerical prediction of the surf-riding threshold of a ship in stern

quartering waves in the light of bifurcation theory

Atsuo MakiÆ Naoya Umeda

Received: 18 January 2008 / Accepted: 1 May 2008 / Published online: 12 August 2008 Ó JASNAOE 2008

Abstract In order to develop design and operational criteria to be used at the International Maritime Organi-zation (IMO), critical conditions for broaching are explored in the light of bifurcation analysis. Since surf-riding, which is a prerequisite to broaching, can be regarded as a het-eroclinic bifurcation, one of global bifurcations, of a surge-sway-yaw-roll model in quartering waves, the relevant bifurcation condition is formulated with a rigorous math-ematical background. Then an efficient numerical solution procedure suitable for tracing the surf-riding threshold hypersurface is presented with successful examples. This deals with all state and control variables in parallel, and excludes backward time integration and an orthogonal condition in the iteration process. The bifurcation condi-tions identified were compared with the results from a direct numerical simulation in the time domain. As a result, it was confirmed that the heteroclinic bifurcation provides a boundary between motions periodically overtaken by waves and nonperiodic motions such as surf-riding and broaching.

Keywords Bifurcation theory
Surf-riding threshold
Heteroclinic bifurcation
Newton method

Following and quartering seas List of symbols

c Wave celerity

Fn Nominal Froude number g Gravitational acceleration

GZ Righting arm H Wave height

Ixx Moment of inertia in roll

Izz Moment of inertia in yaw

Jxx Added moment of inertia in roll

Jzz Added moment of inertia in yaw

Kp Derivative of roll moment with respect to roll rate

Kr Derivative of roll moment with respect to yaw rate

KR Rudder gain

KT Thrust coefficient

Kv Derivative of roll moment with respect to sway

velocity

Kw Wave-induced roll moment

Kd Derivative of roll moment with respect to rudder

angle

K/ Derivative of roll moment with respect to roll angle

L Ship length between perpendiculars m Ship mass

mx Added mass in surge

my Added mass in sway

n Propeller rate

Nr Derivative of yaw moment with respect to yaw rate

Nv Derivative of yaw moment with respect to sway

velocity

Nw Wave-induced yaw moment

Nd Derivative of yaw moment with respect to rudder

angle

N/ Derivative of yaw moment with respect to roll angle

p Roll rate r Yaw rate R Ship resistance

t Time

T Propeller thrust

TD Time constant for differential control

TE Time constant for steering gear A. Maki
N. Umeda (&)

Department of Naval Architecture and Ocean Engineering, Graduate School of Engineering, Osaka University, 2-1 Yamadaoka, Suita, Osaka 565-0971, Japan e-mail: umeda@naoe.eng.osaka-u.ac.jp DOI 10.1007/s00773-008-0017-2

u Surge velocity v Sway velocity

Xw Wave-induced surge force

Yr Derivative of sway force with respect to yaw rate

Yv Derivative of sway force with respect to sway

velocity

Yw Wave-induced sway force

Yd Derivative of sway force with respect to rudder angle

Y/ Derivative of sway force with respect to roll angle

zH Vertical position of centre of sway force due to

lateral motions

v Heading angle from wave direction vc Desired heading angle for autopilot

d Rudder angle / Roll angle k Wavelength

s Time duration for forward integration

nG Longitudinal position of centre of gravity from a

wave trough

1 Introduction

Broaching is a phenomenon in which a ship cannot main-tain a constant course despite maximum steering effort. If the ship’s speed is high enough, the centrifugal force due to this uncontrollable yaw motion could result in capsizing. This phenomenon often occurs when a ship runs in fol-lowing and quartering seas with a relatively high forward speed, especially when it is surf-ridden. Thus it is impor-tant to estimate the surf-riding threshold.

Grim [1] explained that the surf-riding boundary coin-cides with the case where a trajectory from an unstable equilibrium of an uncoupled surge model on a wave is connected to another unstable equilibrium. This is a het-eroclinic bifurcation in nonlinear dynamical system theory. Makov [2] demonstrated the validity of Grim’s statement by using a phase plane analysis. Ananiev [3] obtained an analytical approximate solution by applying a perturbation technique. Umeda and Renilson [4] extended this approach from an uncoupled surge model to a coupled surge-sway-yaw model. Spyrou [5] and Umeda [6] numerically obtained the heteroclinic bifurcation for the uncoupled surge model and the coupled surge-sway-yaw-roll model with a PD autopilot, respectively. It is worth noting that the results from Umeda and Renilson’s mathematical model were systematically compared with those from free-running model experiments and that the model was consequently upgraded to give quantitative agreement with them [7].

Utilising this it is possible to directly predict the surf-riding threshold, which can then be used in design and/or operational criteria for ships.

For the uncoupled surge model in pure following waves, analytical formulae to calculate the surf-riding threshold were presented by Kan [8] and Spyrou [9,10].

For the case of coupled surge-sway-yaw-roll model in stern quartering waves, however, it is necessary to predict the surf-riding threshold numerically because no analytical solution is currently available. Numerical techniques for directly identifying heteroclinic bifurcation have been extensively investigated so far. Literature reviews on the numerical methods can be found in Kuznestov [11] and Kawakami et al. [12] For example, Kawakami [13] in 1981 mathematically formulated an equation set to determine heteroclinic bifurcation using the Runge Kutta method to numerically obtain the Jacobian in his mathematical for-mulation. Here the separatrix loop can be explicitly determined together with information in the locally linear field close to the unstable fixed point. In this technique, apart from the numerical accuracy of the Runge Kutta scheme, no further approximation exists.

Friedman and Doedel proposed a different approach [14,

15], which is now available as HomCont in the software package: AUTO97 [16]. Here, assuming that the state equation is differentiable, the Jacobian is numerically and approximately determined by numerical differentiation and a collocation method.

Despite the progress of these methods, no numerical technique has been utilised for practical purposes with the mathematical model of the coupled surge-sway-yaw-roll model in stern quartering waves using a proportional and differential (PD) autopilot, which is common in the field of ship dynamics. This is probably because the mathematical model validated using model experiments for the current phenomena seems to be too complicated for such theoretical methodology. Therefore, it is important to demonstrate the applicability of such methodology to surf-riding of ships as a real-world problem. Leaving the rudder saturation phenomena and capsizing, which could result in nondiffer-entiable equations, to be incorporated in the future, the authors attempt to apply Kawakami’s approach but with a differentiable mathematical model as the first step.

As a simplest application of Kawakami’s approach [12], Umeda et al. [17] proposed a theoretical formulation of this bifurcation problem with their coupled surge-sway-yaw-roll mathematical model. In their procedure by utilizing the Newton method focussing on the nominal Froude number (based on the Froude number obtained in calm water with the given propeller rpm) the heteroclinic bifurcation is directly determined. Here, the trajectory near a saddle-type surf-riding equilibrium point to reach another saddle is investigated using the Newton method.

However, this algorithm requires us to calculate the effect of many unknown variables sequentially in each iteration. In other words, only the nominal Froude number is

systematically changed in the iterations, and the other vari-ables are consequently changed. Thus, this method is not suitable for obtaining the heteroclinic bifurcation points, or its hypersurface, as functions of other control parameters. Hence, in this research, to explicitly utilise other unknown variables as well in the iterations, the authors proposed the multivalued Newton method in which unknown variables can be simultaneously determined for identifying the het-eroclinic bifurcation points. This approach can be regarded as a more direct examination of the applicability of Kawakami’s approach. As a result, with some practical improvements in the algorithm, it was confirmed that the surf-riding threshold in stern quartering waves can be estimated by one of the established numerical procedures using the mathematical model developed and validated by the authors.

2 Mathematical models

The mathematical model used in this paper is a manoeu-vring model of the surge-sway-yaw-roll motion with a PD autopilot for the prediction of broaching associated with surf-riding in regular following and quartering seas [6]. When a ship is running with a high forward velocity in following waves, the encounter frequency becomes much smaller than its natural heave and pitch frequencies. Therefore, these motions can be estimated quasistatically by using the stable equilibrium in the vertical plane.

The coordinate systems used are defined in Fig.1. O-ngf is a wave-fixed coordinate system with its origin fixed on the calm water level at a wave trough, and GS- xyz is an upright body-fixed coordinate system with its

origin at the centre of gravity, the x-axis pointing toward the bow, the y-axis to starboard, and the z-axis downward. The state vector, x [ R8, and control vector, b [ R2, of this system are defined as follows:

x nð G=k; u; v; v; r; /; p; dÞ T ð1Þ b n; vð cÞ T ð2Þ The dynamical system can be represented by the following state equation: x¼ F x; bð Þ  fð1ðx; bÞ; f2ðx; bÞ; . . .; f8ðx; bÞÞ T ; ð3Þ where f1ðx; bÞ  u cos v  v sin v  cð Þ=k ð4Þ f2ðx; bÞ  Tðu; nÞ  RðuÞ þ X½ wðnG=k; vÞ ð5Þ f3ðx; bÞ  ðm þ m½ xÞur þ Yvðu; nÞv þ Yrðu; nÞr

þ Y/ðuÞ/ þ YdðnG=k; u; v; nÞd

þYwðnG=k; u; v; nÞ=ðm þ myÞ ð6Þ

f4ðx; bÞ  r ð7Þ

f5ðx; bÞ  N
vðu; nÞv þ Nrðu; nÞr þ N/ðuÞ/

þ NdðnG=k; u; v; nÞdþNwðnG=k; u; v; nÞ =ðIZZþ JZZÞ ð8Þ f6ðx; bÞ  p ð9Þ f7ðx; bÞ  m½ xzHurþ Kvðu; nÞv þ Krðu; nÞr þ KPðuÞp þ K/ðuÞ/ þ KdðnG=k; u; v; nÞd þ KwðnG=k; u; v; nÞmgGZð/Þ=ðIxxþ JxxÞ ð10Þ f8ðx; bÞ  d  K½ Rðv  vcÞ  KRTDr=TE: ð11Þ Since the wave-induced forces and moments are periodic functions of the relative position of the ship to a wave, the system described by the above equations can be regarded as nonlinear and autonomous. The solution of the state equations are indicated by w = w(t, XI; b), where XIis the

y

z

O

G

φ

ζ

η

x

y

η

initial value of the state vector x. Based on the aforemen-tioned mathematical model, numerical calculations were carried out for a 135GT Japanese purse seiner used as a subject ship of the International Towing Tank Conference (ITTC) benchmark testing [18]. The principal characteris-tics and the body plan of this vessel are shown in Table1

and Fig.2, respectively. Hydrodynamic coefficients and other related coefficients can be found in the literature [22].

3 Heteroclinic bifurcation

A nonlinear dynamical system described by Eq.3 could have fixed point:

X0¼ Gð 0; u0; v0;v0; 0; /0; 0; d0Þ T

; ð12Þ

where

F Xð 0; bÞ ¼ 0: ð13Þ

In this paper, for the sake of brevity, nG/k is denoted G. The

fixed points which correspond to surf-riding, where the ship runs at wave speed in regular waves, F(x; b) are linearised at X0, putting x = X0+ DX to obtain the

following equation:

D _X¼ DF Xð 0; bÞDX; ð14Þ

where

DF Xð 0; bÞ  ofiðX0; bÞ=oXj: ð15Þ If an eigenvalue of DF(X0; b), ki, which is obtained by

DF Xð 0; bÞ  ki

½ 
DX ¼ 0 ð16Þ

has a positive real part, then the local asymptotic behaviour at x0is unstable.

Hartman’s theorem and the stable manifold theorem [23] enable us to investigate the local topological structure of the system using Eq.3. That is, there exist local stable and unstable manifolds, Ws

locðX0; bÞ and Wlocu ðX0; b Þ;tan-gent to eigenspaces, spanned by DF(X0; b) at X0. Then the

global stable and unstable manifolds Ws and Wu are obtained by letting points in Wlocs flow backward in time

and those in Wlocu flow forward.

The numerical survey [17] of the system described by Eq. 3applied to the subject ship was executed for all discrete operational and environmental parameters used in this paper. Here all the eigenvalues were calculated for all fixed points for the control parameters. This result indicates that there is normally one fixed point having only one eigenvalue with a positive real part, k1, if a fixed point exists [17]. There are

some exceptions which are not relevant to the surf-riding threshold discussed here. Thus, such a fixed point has a one-dimensional unstable invariant manifold and a seven-dimensional stable invariant manifold. For heteroclinic bifurcation the invariant unstable manifold Wu is set to coincide with the stable manifold Ws. As a simplified case shown in Fig.3, X1[ R

8

are also saddle-type equilibrium points defined by Eq. 17.

X1¼ Gð 1; u1; v1;v1; 0; /1; 0; d1Þ T

: ð17Þ

Here the vectors Xa Xu X0; Xx Xs X1 2 R8 represent the unstable eigenvector and the stable eigen-vector, respectively.

4 Calculation algorithm and numerical results

Kawakami [18] developed the procedure to determine the heteroclinic bifurcation point for an autonomous nonlinear dynamical system. As mentioned, normally at a

saddle-Table 1 Principal characteristics of the subject ship

Items Values

Length: LPP 34.5 m

Breadth: B 7.60 m

Depth: D 3.07 m

Draught at forward perpendicular: df 2.50 m

Mean draught: dm 2.65 m

Draught at aft perpendicular: da 2.80 m

Block coefficient: Cb 0.597

Longitudinal position of centre of gravity from the midship: xCG

1.31 m aft

Metacentric height: GM 1.00 m

Natural roll period: T/ 7.40 m

Rudder area: AR 3.49 m2

Rudder aspect ratio: K 1.84

Time constant of steering gear: TE 0.63 s

Proportional gain: KP 1.0

Time constant for differential control: TD 0.0 s

Maximum rudder angle: dMax ±35°

C_L

type equilibrium point, there is one eigenvector in the unstable direction and seven stable eigenvectors in stable directions. Therefore, it is not straightforward to determine a trajectory in the seven-dimensional stable invariant manifold spanned by the seven eigenvectors in stable directions. To overcome this difficulty the orthogonal condition between the left eigenvector of an eigenvalue and the right eigenvector for the other eigenvalues [24] is applied. By applying this condition we can judge whether Xx coincides with the local stable invariant manifold as

shown in Eq.18.

XT_x
h ¼ 0: ð18Þ

Here the left eigenvector, h, satisfies the following condition:

DFTðX1; bÞ  laI



h ¼ 0: ð19Þ

Figure 3 also shows a schematic view of the orthogonal relation between the stable invariant manifold and the left eigenvector of the eigenvalue with positive real part. Then, if we find the propeller rate, n0, which satisfies the

following relationships, this is a heteroclinic bifurcation point. FðX0; n0; vcÞ ¼ 0 2 R 8 _ð20Þ FðX1; n0; vcÞ ¼ 0 2 R 8 _ð21Þ det DFðXð 0; n0; vcÞ  laIÞ ¼ 0 2 R1 ð22Þ ½DFðX0; n0; vcÞ  laI
Xa¼ 0 2 R8 ð23Þ XT_a
Xa d21¼ 0 2 R 1 _ð24Þ XT_x
Xx d22¼ 0 2 R 1 _ð25Þ DFTðX1: n0; vcÞ  laI

h ¼ 0 2 R8 _ð26Þ hT
h  1 ¼ 0 2 R1 _ð27Þ XT_x
h ¼ 0 2 R1 _ð28Þ Wðs; X0þ Xa; n0; vcÞ  Wðs; X1þ Xx; n0; vcÞ ¼ 0 2 R 8 ð29Þ Here d1 1 and d2 1 are prescribed values. In this

equation set there are 43 unknown variables, as shown in Eq.29. Y XT 0; X T 1; la; XTa; X T x; h; s; n0  T 2 R43 _ð30Þ Since the above procedure is based on Kawakami’s work, it is necessary to obtain the values of Ws and Wu by numerically backward and forward integrating Eq.3 from the equilibrium point with a small perturbation along the stable and unstable invariant manifold, respectively. However, although Xa can be determined using

Kawakami’s procedure, numerical backward integration

along the stable invariant manifold is not always successful. This is because the numerical error increases using backward integration since the trajectory is converging to a stable invariant manifold around the neighbourhood of a fixed point, as explained theoretically by the authors [25]. Thus, an alternative is to determine a bifurcation parameter with which the unstable invariant manifold from a saddle coincides with the local stable invariant manifold of another saddle in its vicinity. Based on this the following relationship for a heteroclinic bifurcation point was obtained:

FðX0; n0; vcÞ ¼ 0 2 R8 ð31Þ FðX1; n0; vcÞ ¼ 0 2 R 8 _ð32Þ det DFðXð 0; n0; vcÞ  laIÞ ¼ 0 2 R 1 ð33Þ ½DFðX0; n0; vcÞ  laI
Xa¼ 0 2 R8 ð34Þ XT_a
Xa d21 ¼ 0 2 R 1 _ð35Þ DFTðX1; n0; vcÞ  laI

h ¼ 0 2 R8 ð36Þ hT
h  1 ¼ 0 2 R1 _ð37Þ WTðs; Xa; n0; vcÞ
h ¼ 0 2 R1 ð38Þ Wðs; X0þ Xa; n0; vcÞ  X1 k k2d2 2 ¼ 0 2 R 1_{: ð39Þ}

Here note that the orthogonal condition in the original algorithm, i.e. Eq.18, is replaced by Eq.38. By applying this equation we can judge whether a trajectory approaches another saddle-type equilibrium point along the stable invariant manifold.

Define the unknown variable, Y, to be solved as follows: Y XT 0; X T 1; la; XTa; h; s; n0  T 2 R35 _ð40Þ

Since dim Y is 35, 35 equations are required. Therefore the vector equation W(Y) [ R35 is defined, and allowed to be zero: W(Y) = 0

x

x

x

x

Fig. 3 Schematic view of heteroclinic bifurcation and orthogonal relation in a three dimensional phase space. Here S denotes the surface spanned by the stable eigenvectors

W Yð Þ  FðX0; n0;vcÞ FðX1; n0;vcÞ det DFðXð 0; n0;vcÞ  laIÞ DFðX0; n0;vcÞ  laI ð Þ
Xa ½ _i for i¼ 1; . . .; 7 XT_a
Xa d21 DFTðX1; n0;vcÞ  laI  
h
 i for i¼ 1; . . .; 7 hT
h  1 WTðs; Xa; n0;vcÞ
h Wðs; X0þ Xa; n0;vcÞ  X1 k k2d2 2 8 > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > : ð41Þ In this paper, the above equation set is numerically solved by the Newton method. The Jacobian matrix for the nonlinear dynamical system (Eq.41) evaluated at a point Y is defined as A(Y). Suppose Y(j) is already close to a solution, then the Newton iteration is given by following equations:

Yðjþ1Þ¼ Yðjþ1Þþ a
H ; j¼ 0 ; 1 ; 2 ; . . . ð42Þ A Y ð Þj
H ¼ W Y ð Þj_{: ð43Þ}

Here, let a [ R1be a positive number representing the step size. This step size should be small enough for stable convergence to the solution. Numerical calculations were carried out using the newly proposed algorithm, and an example of the results is shown in Fig.4. The wave steepness is 0.05 and the wavelength to the ship length ratio is 1.0. In this case the heteroclinic bifurcation point obtained is the nominal Froude number, Fn, of 0.3329 for the autopilot course of 5° from the wave direction. Below this propeller rate the ship is overtaken by waves and above this the ship is captured by a down wave slope. The iterations required by this algorithm to converge to Fn = 0.3329 in the above sea state are shown in Fig.5.

5 Further improvement

To determine the heteroclinic bifurcation point more effi-ciently, the new algorithm was further improved. The authors noticed that orthogonal condition Eq. 18is almost always satisfied when a trajectory approaches a saddle-type equilibrium point. Thus, this orthogonal condition is

-1 -0.5 0 0.5 1 ξG/ λ ξ_G/ λ 0.6 0.8 1 1.2 U C / Surge Velocity Fn = 0.3329 -1 -0.5 0 0.5 1 ξG/ λ 0.6 0.8 1 1.2 Uξ C /

Component of Velocity in the Wave Direction Fn = 0.3329 -1 -0.5 0 0.5 1 ξG/ λ -0.03 -0.02 -0.01 0 0.01 V C / Sway Velocity Fn = 0.3329 -1 -0.5 0 0.5 1 -4 -2 0 2 φ ). ge d( Roll Fn = 0.3329 -1 -0.5 0 0.5 1 ξG/ λ 0 10 20 30 χ ). ge d( Yaw Fn = 0.3329 -1 -0.5 0 0.5 1 ξG/ λ -25 -15 -5 5 δ ). ge d( Rudder Fn = 0.3329

Fig. 4 An example of the heteroclinic bifurcation for wave steepness of 0.05, wavelength-to-ship-length ratio of 1.0, autopilot course of 5°, and rudder gain of 1.0

excluded from the iteration and only checked after con-vergence. Additionally, X1was calculated from Eq.44for

the sake of brevity.

X1¼ Gð 0 1; u0; v0;v0; 0; /0; 0; d0Þ T

: ð44Þ

Then unknown variables are reduced to Eq.45. Y XT

0; la; XTa; n0

 T

2 R18 _ð45Þ

Since dim Y is now 18, 18 equations are required. Therefore the vector equation is defined as W(Y) [ R18, which is assumed to be zero: W(Y) = 0.

W Yð Þ  FðX0; n0;vcÞ det DFðXð 0; n0;vcÞ  laIÞ DFðX0; n0;vcÞ  laI ð Þ
Xa ½ _i for i¼ 1; . . .; 7 XT_a
Xa d21 Wðs; X0þ Xa; n0;vcÞ  X1 k k2d2 2 8 > > > > > > > < > > > > > > > : ð46Þ Figure5 shows the convergence behaviour for a nominal Froude number of 0.3329 utilising the above improved procedure with a wave steepness of 0.05, a wavelength-to-ship-length ratio of 1.0, and an autopilot course of 5°. The number of iterations required is much less than that shown in Fig.5, resulting in a method that is suitable for practical application.

This method was applied to different autopilot course and wave conditions. First the bifurcation point is estimated for the autopilot course of 0° and then the bifurcation points for other autopilot course are continuously traced as shown

in Figs.6,7,8and9. The bifurcation points obtained were compared with results obtained from the numerical simu-lation in the time domain based on the sudden change concept [19]. Because of the dependence of the initial value on nonlinear ship motions such as surf-riding, it is necessary to specify realistic initial conditions even for simple numerical simulation. Thus, in this concept, the initial value is selected to be a periodic state under certain and safe control parameters and then the control parameters are requested to change to different control parameters. This is based on the fact that a surf-riding zone that is dependent on the initial conditions cannot appear when the ship initial state is a periodic state for the critical nominal speed of heteroclinic bifurcation or less [20]. A similar method was utilised by Spyrou [21] focussing on transient phenomena from stable surf-riding equilibrium. In this numerical sim-ulation in the time domain the initial state is fixed at the periodic state of Fn = 0.1 and vc= 0.0° and then

numer-ical integration of Eq.3 was executed with the specified nominal Froude number and autopilot course for 1,000 s. The specified autopilot course changed from the initial value of 0° to 40° with increments of 0.5° and the specified nominal Froude number was 0.1 to 0.5 with increments of 0.01. Thus 6,400 numerical runs were executed. The time series obtained were categorised into: periodic motions; surf-riding; broaching; and capsizing using criteria devel-oped by Umeda and Hashimoto [22]. The heteroclinic bifurcation points obtained using the present method predict the boundary between the periodic motions, with the vessel being overtaken by the waves, and other motions such as surf-riding reasonably well, at least for smaller autopilot courses. Therefore, the present method for identifying the surf-riding threshold can be used as an alternative to the time-consuming numerical simulation.

When the autopilot course is increased, the Froude num-ber for the surf-riding threshold is also increased. When the wave steepness, H/k, is increased, the Froude number for the surf-riding threshold is decreased because of the increase in the wave-induced surge force. For a wave steepness of 0.05 or below, stable or oscillatory surf-riding occurs above the heteroclinic bifurcation points. Here stable surf-riding means that a ship is captured by a stable equilibrium point, where the real parts of all the eigenvalues are negative. If one eigenvalue has a positive real part and nonzero imaginary parts, the riding equilibrium could be oscillatory riding. The boundary between stable and oscillatory surf-riding can be identified as the Hopf bifurcation, as discussed by Umeda et al. [17]. This phenomenon was experimentally identified by Kan et al. [26], and then explained by Spyrou [27]. For much larger wave steepness, such as 0.1, broaching and/or capsizing occur. This is because broaching could occur once surf-riding happens under such wave conditions. Thus, the heteroclinic bifurcation can be used as a threshold

0 200 400 600 800 Iterations 0.32 0.324 0.328 0.332 0.336 Fn Original procedure Improved procedure

Fig. 5 Convergence process of the original procedure for nominal Froude number of 0.3329

for broaching. However, it is noteworthy that the heteroclinic bifurcation itself does not distinguish broaching from surf-riding.

Although the sudden change concept used here for the numerical simulation is designed to minimise the depen-dence on the initial conditions, the small disagreement between the bifurcation and the surf-riding threshold could be explained as follows. For larger wave steepness, a

saddle can be connected by a trajectory with another saddle that is more than one wavelength away in a special case. This problem could be solved if we formulated the het-eroclinic condition for all possible combinations of two saddles, but this could be prohibitively complicated. Moreover, broaching periodically occurs in a very limited region. Here, the ship is forced to surf-ride, and then broaches but does not capsize. Because of the large head-ing angle due to broachhead-ing, the ship is overtaken by a wave but is then forced to surf-ride again in the down slope of the following wave. Thus, this phenomenon can be categorised as a periodic state but is different from oscillatory surf-riding. Its period can be a multiple of the encounter period [28].

In Fig.8 the bifurcation points cannot be continuously traced at the autopilot course of 5°. By using the phase plane analysis, we determine a bifurcation point for an autopilot course of 10° and then utilised this as the initial value for tracing bifurcation points. As a result, an arc of bifurcation points from 10° is obtained. At this stage, the reason for failing to find the bifurcation points between 5° and 10° is not known, and this task remains for the future. Umeda et al. [17] also compared the results from direct numerical simulation in the time domain with those using bifurcation analysis for the same conditions, similar to Figs. 6,7,8,9in this paper. However, Figs.6,7,8,9differ from the figures in Umeda et al. [17]. For example, whereas the previous method can only be used to generate results up to an autopilot course of 2° for H/k = 0.04, k/L = 1.5, the present method can lead more points, as shown in Fig.8. In contrast, for the case of H/k = 0.1 and k/L = 1.637 shown in Fig.9 the number of points obtained from the bifurca-tion analysis is substantially smaller. The bifurcabifurca-tion points obtained by the previous method but not in the present method are mainly inside the broaching or capsizing area, so it would appear that the present results are more rea-sonable. Therefore, we can conclude that the present estimation is an improvement over the previous one.

Fig. 6 Numerical simulation and heteroclinic bifurcation for: H/k = 0.05, k/L = 1.0

Fig. 7 Numerical simulation and heteroclinic bifurcation for: H/k = 0.05, k/L = 1.5

Fig. 9 Numerical simulation and heteroclinic bifurcation for: H/k = 0.1, k/L = 1.637

Fig. 8 Numerical simulation and heteroclinic bifurcation for: H/k = 0.04, k/L = 1.5

6 Concluding remarks

This paper examined the applicability of a numerical method for identifying heteroclinic bifurcation to the surf-riding threshold prediction with the coupled surge-sway-yaw-roll mathematical model using a PD autopilot. Some modifications were made to obtain robustness and effi-ciency for convergence. The method is made sufficiently efficient by simplifying the iteration process, and its limi-tations due to unexpected physical phenomena are noted.

Acknowledgments This research was partly supported by a Grant-in-Aid for Scientific Research of the Japan Society for Promotion of Science (no. 18360415). The authors are grateful to Dr. M. Renilson from the Australian Maritime College for his technical and language advice for preparing this paper.

References

1. Grim O (1951) Das Schiff von Achtern Auflaufender See, Jahr-buch der Schiffbautechnischen Gesellschaft, 4 Band, pp 264–287 2. Makov Y (1969) Some results of theoretical analysis of surfriding

in following seas (in Russian). T krylov Soc 126

3. Ananiev DM (1966) On surf-riding in following seas (in Russian). Trans Krylov Soc 13:169–176

4. Umeda N, Renilson MR (1992) Broaching—a dynamic system analysis of yaw behaviour of a vessel in following and quartering seas. In: Wilson PA (ed) Manoeuvring and control of marine craft. Computational Mechanics, Southampton, pp 533–543 5. Spyrou KJ (1996) Dynamic instability in quartering seas: the

behaviour of a ship during broaching. J Ship Res 40(1):46–59 6. Umeda N (1999) Nonlinear dynamics on ship capsizing due to

broaching in following and quartering seas. J Mar Sci Technol 4:16–26

7. Umeda N, Hashimoto H (2006) Recent developments of capsiz-ing prediction techniques of intact ships runncapsiz-ing in waves. Proceedings of the 26th ONR symposium on naval hydrody-namics, Rome

8. Kan M (1989) Surging of large amplitude and surf-riding of ships in following seas (part 3 phase plane analysis) (in Japanese). J Naval Arch Jpn 166:267–276

9. Spyrou KJ (2001) Exact analytical solutions for asymmetric surging and surf-riding. In: Proceedings of 5th international workshop on stability and operational safety of ships, University of Trieste, Trieste, pp 4.51–4.53

10. Spyrou KJ (2006) Asymmetric surging of ships in following seas and its repercussions for safety. Nonlin Dyn 43:149–172 11. Kuznetsov YA (2004) Elements of applied bifurcation theory, 3rd

edn. Springer, Heidelberg, pp 581–585

12. Kawakami H, Yoshinaga T, Ueta T (1997) Methods of computer simulation on dynamical systems (in Japanese). Bull Jpn Soc Ind Appl Math 7(4):49–57

13. Kawakami H (1981) Calculating method of Separatrix loop of saddle-typed fixed points (in Japanese). J Inst Electron Inf Commun Eng J64-A(10):860–861

14. Doedel E, Friedman M (1989) Numerical computation of het-eroclinic orbits. J Comput Appl Math 26:159–170

15. Doedel E, Friedman M (1991) Numerical computation and continuation of invariant manifolds connecting fixed points. SIAM J Numer Anal 28(3):789–808

16. AUTO97 (1997)http://indy.cs.concordia.ca/auto/

17. Umeda N, Hori, M, Hashimoto H (2006) Theoretical prediction of broaching in the light of local and global bifurcation analysis. In: Proceedings of the 9th international conference on stability of ships and ocean vehicles. Federal University of Rio de Janeiro (Rio de Janeiro) 1:353–362

18. Umeda N, Renilson MR (2001) Benchmark testing of numerical prediction on capsizing of intact ships in following and quartering seas. In: Proceedings of the 5th international workshop on sta-bility and operational safety of ships, University of Trieste, pp 6.1.1–10

19. Umeda N (2000) Application of nonlinear dynamical system approach to ship capsize due to broaching in following and quartering seas. In: Contemporary ideas on ship stability. Else-vier, Amsterdam, pp 57–68 (This paper was firstly read at the first international ship stability workshop held at Ross Priory, UK in 1995)

20. Umeda N (1990) Probabilistic study on surf-riding of a ship in irregular following seas. In: Proceedings of the 4th international conference on stability of ships and ocean vehicle, Naples, pp 336–343

21. Spyrou KJ (1996) Dynamic instability in quartering waves part II: analysis of ship roll and capsize for broaching. J Ship Res 40(4):326–336

22. Umeda N, Hashimoto H (2002) Qualitative aspects of nonlinear ship motions in following and quartering seas with high forward velocity. J Mar Sci Technol 6:111–121

23. Guckenheimer P, Holmes P (1983) Nonlinear oscillations, dynamical systems and bifurcations of vector fields. Springer, New York, pp 12–16

24. Luenberger DG (1979) Introduction to dynamic systems. Wiley, New York

25. Maki A, Umeda N, Hori M (2007) Prediction of global bifurca-tion points as surf-riding threshold in following seas (in Japanese). J Soc Naval Arch Jpn 5:205–215

26. Kan M, Saruta T, Yasuno M, Yamakoshi Y, Suzuki S (1987) Surging of large amplitude and surf-riding of ships in following seas. J Soc Naval Arch Jpn 187:152–162

27. Spyrou KJ (1995) Surf-riding, yaw instability and large heeling of ships in following/quartering waves. Schiffstechnik 42:103– 112

28. Umeda N, Matsuda A (2000) Broaching in following and quar-tering seas—theoretical attempts and new prevention device. In: Proceedings of the 7th international conference on stability of ships and ocean vehicles, vol A. Australian Maritime Engineering CRC, Launceston, pp 460–470