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.000510 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.210 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 ik  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