Nonlinear Container Ship Model for the Study of Parametric Roll Resonance

Parametric roll is a critical phenomenon for ships, whose onset may cause roll oscillations up to ±40◦, leading to very dangerous situations and possibly capsizing. Container ships have been shown to be particularly prone to parametric roll resonance when they are sailing in moderate to heavy head seas. A Matlab/Simulinkr parametric roll benchmark model for a large container ship has been implemented and validated against a wide set of experimental data. The model is a part of a Matlab/Simulink Toolbox (MSS, 2007). The benchmark implements a 3-order nonlinear model where the dynamics of roll is strongly coupled with the heave and pitch dynamics. The implemented model has shown good accuracy in predicting the container ship motions, both in the vertical plane and in the transversal one. Parametric roll has been reproduced for all the data sets in which it happened, and the model provides realistic results which are in good agreement with the model tank experiments.

1 Introduction Parametric roll is an autoparametric resonance phenomenon whose onset causes a sudden rise in roll oscillations.The resulting heavy roll motion, which can reach 30-40 degrees of roll angle, may bring the vessel into conditions dangerous for the ship, the cargo, and the crew.The origin of this unstable motion is the time-varying geometry of the submerged hull, which produces periodic variations of the transverse stability properties of the ship.
Parametric roll is known to occur when a ship sails in moderate to heavy longitudinal or oblique seas; the wave passage along the hull and the wave excited vertical motions result in variations of the intercepted waterplane area, and in turn, in relevant changes in the restoring characteristics.The onset and build-up of parametric roll is due to the occurrence of concomitant conditions: the wave length is close to the ship length (λ w ≈ L PP ), the ship approaches waves with encounter frequency almost twice the roll natural frequency (ω e ≈ 2ω 0 ), and the wave height is greater than a ship-dependent threshold (h w > h s ).
The risk of parametric roll has been known to the maritime community since the early fifties, but only for small vessels with marginal stability -e.g.fishing boats -sailing in following seas.However, the phenomenon has recently attracted significant interest by the scientific community after accidents occurred with container ships sailing in head seas, incidents that involved significant damage to cargo as well as structural damages for millions of dollars (France et al., 2001;Carmel, 2006).
Several different types of vessel have reported to experience parametric roll in head seas, e.g.destroyers (Francescutto, 2001), ro-ro paxes (Francescutto and Bulian, 2002) and PCTC (Palmquist and Nygren, 2004).Container carriers, however, are the most prone to parametric roll because of the current particular hull shape, i.e. large bow flare and stern overhang, and hence abrupt variation in the intercepted water-plane area when a wave crest or trough is amidships.
This has called for deep investigations into the nature of parametric roll in head/near head seas, and for the development of mathematical models able to capture and reproduce the physical aspects driving the resonant motion.In the last six years mathematical models of different complexity have been proposed by the scientific community, most of them relying on the Mathieu Equation to describe the dynamics of the ship subject to parametric resonance.
One-DOF models considering the uncoupled roll motion have been widely used to analyze the critical parameters of the phenomenon and derive stability conditions.Examples can be found in the papers by France et al. (2001) and Shin et al. (2004) where the authors employed the 1-DOF roll equation to show that, in regular waves, the Mathieu Equation can explain the onset of heavy roll motion in head seas.Bulian (2006) proposed a 1.5-DOF model where the dynamic interaction between the vertical motions and the roll oscillation was relaxed by the assumption of quasi-static heave and pitch.Moreover, that assumption allowed an analytical description of the GZ curve that was approximated as a surface varying with roll angle and wave crest position.This model is considered valid for moderate ship speed in head seas, and has lead to reasonable results in predicting parametric roll.
A 3-DOF nonlinear fully coupled model was first developed by Neves (2002).A first attempt was done by using Taylor series expansion up to 2 nd -order to describe the coupled restoring forces and moments in heave, pitch and roll.This model, although it provided a quite thorough description of the nonlinear interactions among the different modes, tended to overestimate the roll oscillation above the stability threshold.Neves and Rodríguez (2005) proposed a 3 rd -order analytical model where the couplings among the three modes are expressed as a 3 rd -order Taylor series expansion.In this new model the nonlinear coefficients are mathematically derived as a functions of the characteristics of the hull shape.This 3-DOF model has been applied for the prediction of parametric roll to a transom stern fishing vessel (Neves and Rodríguez, 2006a,b) providing outcomes which better match the experimental results than the 2 nd -order model.
It is noted that the above-mentioned literature have attempted to model parametric roll from an analytical points of view.Jensen (2007) takes a statistical approach instead, motivated by the difficulties inherent in describing the interaction between a 3-dimensional wave pattern and the motion of a ship hull.He shows how the statistical distribution of nonlinear ship responses can be estimated very accurately using a firstorder reliability method.A commercial implementation in a system to predict parametric roll (SeaSense ) was reported in Nielsen et al. (2006).
The direction of this paper is the analytical one, aiming at providing simulation tools that could e.g.be used in studies of active stabilization and control.The model proposed by Neves and Rodríguez (2005) is applied to describe the dynamics of a container vessel subject to parametric roll resonance conditions.The model parameters are identified based upon the ship line drawings and the loading conditions.A Matlab/ Simulink implementation of the above model is then presented.The reliability of the implemented model in simulating parametric resonance behavior is validated against experimental data.The validation has shown good agreement with the experimental results for roll both in the experiments where parametric roll resonance occurred, and in the experiments where it did not occur.
The main goal of this work is to provide a benchmark for simulating parametric roll of a container ship over a large range of ship speeds and sea states.This benchmark has been designed to be a fully integrated part of Matlab/Simulink Toolbox for marine systems (MSS, 2007).The availability of such a powerful tool opens up a great wealth of opportunities, notably the design and testing of novel model-based roll motion stabilizers.

Mathematical Model for Parametric Roll
The proposal and the adoption of an analytical model for representing a specific phenomenon should be driven by a trade-off between complexity and agreement with physical laws governing that phenomenon and/or experimental results.Tondl et al. (2000) define an autoparametric system as follows: Definition 1 Autoparametric systems are vibrating systems that consist of at least two constituting subsystems.One is the primary system that will usually be in a vibrating state.This primary system can be externally forced, self-excited, parametrically excited, or a combination of these.The second constituting subsystem is called the secondary system.The secondary system is coupled to the primary system in a nonlinear way, but such that the secondary system can be at rest while the primary system is vibrating.
An autoparametric system is, hence, characterized by these main aspects: 1. two nonlinearly coupled subsystems; 2. a normal mode where the primary system is in a vibrating state and the secondary system is at rest; 3. the presence of instability regions where the normal mode becomes unstable; 4. in the region of instability of the normal mode the overall system is in autoparametric resonance: the secondary system is parametrically excited by the vibrations of the primary system and it will not be at rest anymore.
Considering Definition 1, 1 DOF models have too little complexity to describe an autoparametric system, since the roll motion for a ship sailing in longitudinal seas represents only the secondary system.They are useful to obtain insight in the parametric roll resonance phenomenon, but they will have difficulty predicting the real amplitude of the oscillations about the transverse plane.
The model proposed by Neves and Rodríguez (2005) is complex enough to capture the dynamics of a container vessel behaving as an autoparametric system; it includes both the primary system (heave and pitch dynamics) which is externally excited by the wave motion, and the secondary system (roll dynamics) which is parametrically excited by the primary.

Equations of Motion
The 3-DOF nonlinear mathematical model of the container vessel is presented in the following way (using the notation of Neves and Rodríguez (2005)): Let be the generalized coordinate vector, where z is the heave displacement, φ is the roll angle, and θ is the pitch angle, as shown in Figure 1.Then the nonlinear equations of motion can be expressed in matrix form as where • M ∈ R 3×3 is the diagonal rigid-body generalized mass matrix; • A ∈ R 3×3 is the generalized added mass matrix; • B ∈ R 3×3 is the hydrodynamic damping (nonlinear in roll); • c res ∈ R 3 is the nonlinear vector of restoring forces and moments expressed as functions of the relative motion between ship hull and wave elevation ζ(t); • c ext ∈ R 3 is the vector of the external wave excitation forces and moments which depends on wave heading, encounter frequency, wave amplitude and time.

Generalized Mass, Added Mass and Damping
The generalized mass matrix can be written as where m is the ship mass, I x is inertia in roll and I y is inertia in pitch.
The hydrodynamic added mass and damping matrices are expressed as where all entries except K φ( φ) can be evaluated by means of potential theory (Salvesen et al., 1970).
The hydrodynamic damping in roll may be expressed as where the linear term represents the potential and linear skin friction, whereas the nonlinear term takes into account viscous effects.The coefficients K φ and K φ| φ| can be calculated by the formulae given in Himeno (1981).The roll damping characteristics may also be derived from data of roll decaying tests at appropriate forward speeds of the vessel.

Waves
In regular seas, the incident wave elevation according to the Airy linear theory, see Newman (1977), is defined as where A w is the wave amplitude, k is the wave number, χ is the wave heading, and ω e is the encounter wave frequency.For head seas (χ = 180 • ), the wave elevation reads as

Nonlinear Restoring Forces and Moments
The nonlinear restoring actions are given by the combination of the effects of the vessel motion in calm water and the effect of the wave elevation along the hull.Therefore, the vector of restoring forces and moments can be written, up to 3 rd -order terms, as where c pos,s i ζ j = ∂ i+j cpos ∂s i ∂ζ j s i ζ j .The 1 st , 2 nd and 3 rd -order components in (9), which are independent of the displacement vector s, must be included in the external forces and moments acting on the vessel.These terms describe the linear and nonlinear Froude-Krylov forces/moments.
The 2 nd and 3 rd -order nonlinear effects due to hullwave interactions must, instead, be included in the restoring vector c res because of their affinity, from the mathematical point of view, with the hydrostatic actions.Then the restoring force and moments due to body motion are given by where Therefore the restoring force/moments in each degree of freedom are given by the following terms: • 1 st -order body motions (c pos,s ) • 2 nd -order hull-wave interactions (c pos,sζ ) The time varying terms depend explicitly on the wave elevation ζ(t) and thus implicitly on the time t.
Looking at the 1 st , 2 nd and 3 rd -order coefficients, a strong cross-coupling between all three degrees of freedom becomes evident.

External Forcing
The interaction between ship motion and wave passage is modeled as a variation of the geometry of the submerged hull defined by the instantaneous wave position.The external forcing vector c ext (ζ, ζ, ζ) includes only contributions independent of ship motions, such that τ 1w represents the 1 st -order wave excitation forces generated by the wave motion.These forces are characterized by two contributions: the first one is due to Froude-Krylov forces, which are caused by incident waves considering the hull restrained from moving and that the presence of the hull does not influence the wave field.The second contribution gives the diffraction forces, which provide the corrections necessary for the variation of the flow field produced by the hull.
τ 2w are the 2 nd -order wave excitation forces which include three important components.The first contribution is given by the mean wave drift forces caused by nonlinear wave potential effects; the second one is due to low-frequency wave drift forces caused by nonlinear elements in the wave loads; and the third component is given by high-frequency wave drift forces.
In the present analysis the external force and moments are defined as being proportional to the first order wave motion, whereas higher order terms are neglected.Therefore the external force/moments vector c ext reads as The wave excitation forces are defined by the waveforce response amplitude operator (force RAO) for each degree of freedom.The Force RAO is computed (Perez, 2005) as where τ1wi is the complex 1 st -order wave excitation forces, and ζ is the complex wave elevation.Since the model only considers head seas, (18) simplifies to With these force RAOs, it is possible to obtain the wave excitation loads in each degree of freedom as for i = 3, 4, 5, where α i = arg[ Fi (ω e )].For example, the external force acting on heave is given by 3 Identification of Model Parameters from Hull Form and Wave Characteristics The identification of model parameters is completely based upon the hull shape of the container vessel and upon the wave characteristics.In this section the formulas are presented.The numerical values of those parameters, computed for the considered container ship, can be found in Appendix A.
In Table 1 the main characteristics of the containership are reported.

Body Motion Coefficients
The 1 st -order body motion coefficients refer to calm water hydrostatics and are given by where ρ is the water density, g is the acceleration of gravity, A 0 is the waterplane area, x f0 is the longitudinal coordinate of the centre of floatation, and GM l is the longitudinal metacentric height.The 2 nd and 3 rd -order body motion coefficients correspond to the variations in the restoring characteristics of the ship due to the changes in pressure related to the vessel motions.In order to compute them numerically, it is necessary to express the nonlinear hydrostatic actions as function of the three modes heave, pitch, and roll.In particular, it is possible to demonstrate that where ∇ 0 is the mean displacement, is the instantaneous displacement, z G is the vertical coordinate of the centre of gravity, x B1 , y B1 , and z B1 are the coordinates of the instantaneous centre of buoyancy.
Tables 2-3 show the 2 nd and 3 rd -order coefficients for each degree of freedom.

Hull-Wave Interaction Coefficients
Under the assumption of regular waves, the periodic wave passage along the hull produces cyclic variation in the restoring characteristics of the vessel.These changes are taken into account by the 2 nd and 3 rdorder coefficients included in the nonlinear interactions c pos,sζ and c pos,s 2 ζ + c pos,sζ 2 .In order to determine the hull-wave interaction coefficients, the Froude-Krylov forces must be defined.The velocity potential for the undisturbed wave, as defined in ( 7), is given by A w g ω e e kz sin(kx cos χ − ky sin χ − ω e t).(24) Therefore, the 1 st and 2 nd -order Froude-Krylov forces are: where n is the normal to the hull surface and j addresses the specific mode for which the force is computed.The coefficients are then given by the formulas is Tables 4-5.
Due to the assumption of regular waves, the coefficients can be described as a sum of a sine and a cosine term.For instance, the 2 nd -order term K ζφ (t), which is proportional to wave amplitude, can be written as where K ζφc and K ζφs are constants.
Analogously, the 3 rd -order terms K ζzφ (t) and K ζφθ (t) are given by These functions play an important role since they parametrically excite the coupled system, being multiplied with, respectively, z(t)φ(t) and φ(t)θ(t).
The 3 rd -order term K ζζφ (t), which is proportional to the wave amplitude squared, is given by Heave-roll-pitch coupling where it can be noticed the presence of a constant term plus a super-harmonic term of double the encounter frequency.

Matlab Implementation of the Model
A Matlab/Simulink model for the container ship model was developed for the purposes of simulating parametric roll resonance, based on the model of Section 2.
For each time instant and system state, a function generates the instantaneous value of [ṡ T sT ] T .Numerically integrating with an explicit Runge-Kutta method of order 4, with the fixed time step h = 1 s, the state [s T ṡT ] T is calculated for any given time instant.
The parameters used in the calculations are listed in Appendix A. While this was not done for the results presented in this paper, for other encounter frequencies than the ones used in the experiments, interpolation can be applied to calculate approximate parameter values.
The code is part of the Marine Systems Simulator (MSS, 2007).

Validation of the Model Against Experimental Data
A comparison of the simulation and the experimental results can be seen in Figures 2-24.
The experiments were conducted with a 1:45 scale model ship in a towing tank.The experiments were done with varying forward speed, and wave frequency and height.This is summarized in Table 6.All data in the table and in the figures are in full scale.
The simulations were done with the code described in Section 4. All simulations were made ballistically.Initial conditions can be found in Table 7.Initial conditions not listed in the table (θ, ż, φ and θ) were all zero.The experiments were all assumed to start at t = 0.
A comparison of the simulation results with the experimental results can be seen in Table 8.The first column is the experiment number.The second column is wave amplitude A w .The third column is wave frequency ω.The fourth column is the ratio of encounter frequency to natural roll frequency (ω e /ω 0 ).The fifth     In Figure 24, we can see the maximum roll angle achieved in the simulations and experiments for certain conditions, plotted against the ratio of encounter frequency to natural roll frequency (ω e /ω 0 ).The data in the figure is all for A w = 2.5 m, and ω = 0.4640 rad/s.

Analysis of the Model Based Upon the Validation Results
The 3 rd -order model developed for the 281m long container ship shows high capabilities in reproducing the vertical and transversal dynamics of the vessel under parametric resonance conditions, as shown by the comparison of the experimental results (Figures 2-24).
Considering the 13 experiments where parametric resonance did occur, the implemented model performs well: starting from similar initial conditions and being subjected to the same excitation forces used during the experiments, the model develops parametric resonance within the same time frame as the 1:45 scale model ship in most of the cases.
The most obvious differences between the simulation and the experimental results consists of the amplitude of the oscillations.In all the experiments where parametric resonance occurred, the peak value of the roll oscillations is higher than the saturation level at which the model settles.Although the model has a general tendency to underestimate the peak value of the roll motion, the gap is relatively small in most cases.
Considering the 9 experiments where parametric roll did not occur, the model produced 5 false positive cases developing resonant motion.In order to understand this disagreement between model behavior and experimental results, the tuning factor ω e /ω 0 must be taken into consideration.In fact all the 5 false-positive cases occur with a tuning factor close to the limits of the first instability region of the Hill-Mathieu Equation (ω e ≈ 2ω 0 ), as shown in Figure 24.Looking at the peak value of the roll oscillations (Figure 24 and Table 8), it is seen that the largest differences are in the region of high tunings ( ωe ω0 2.1) for which the model predicts large roll motion whereas the experiments showed no amplification.It seems obvious that when the experimental conditions are close to the limits of stability the model does not match exactly the frequency at which the abrupt variation in roll motion take place.
The errors indicated in tests 1173, 1179 and 1181 have no real physical meaning, since the initial condition of 2 degrees was chosen arbitrarily high in order to indicate a decaying motion.
For all 22 experiments, heave and pitch dynamics have shown relatively good agreement with the experiments.In all the test runs the two modes oscillates at the excitation frequency, matching the experimental records.The amplitude of the oscillations is close to that of the experimental values.

Conclusions
A Matlab/Simulink benchmark for the simulation of parametric roll resonance for a large container ship has been implemented and validated against experimental results.The implementation reflects the coupled 3 rdorder nonlinear model for parametric roll first introduced by Neves and Rodríguez (2005).
The mathematical model for the container ship at hand, already presented in Rodríguez et al. (2007), has been reviewed, illustrating in details the ship dynamics that has been taken into account.1 st , 2 nd and 3 rd -order contributions have been described and analytical formulas of all the couplings coefficients due to heave, roll, pitch, and wave motion have been given.Furthermore, numerical values of all coefficients are computed based upon the hull geometry and wave characteristics.
The benchmark has been tested on a set of 22 different conditions, which have been chosen to match the experimental conditions.Each test run differs from the previous for at least the value of one parameter among ship speed, wave frequency, and wave height.The heave and pitch dynamics described by the model are in good accordance with the experimental results.The model seems to reproduce the pitch motion slightly better, catching the right amplitude in most of the cases.
The results obtained in roll have shown good agreement with the records of the experiments.In particular, the model agrees with the experiments in all the cases where parametric roll occurred, although the amplitude of the roll oscillations does not quite reach the experimental peak value in most cases.
For the experiments where there was no parametric roll, then the model produces false positives in about 50% of the cases.This disagreement between the simulation and the experimental results is believed to be related to the specific values of the tuning factor ω e /ω 0 which were too close to the limits of the first instability region of the Hill-Mathieu Equation.In these cases, it is very difficult to get the correct response with ballistic simulations.
The availability of this benchmark offers a wide range of opportunities for development of new model-based control strategies for counteracting or preventing parametric roll resonance.

A Tables of Coefficients
The parameters can be found in Tables 9-16.All numbers are given in the kg-m-s (SI) system.Only non-zero numbers are listed.
Table 9 contains the rigid body inertia matrix, and

Table 1 :
Main characteristics of the container ship

Table 4 :
2 nd -order hydrostatic restoring coefficients due to wave passage

Table 5 :
3 rd -order hydrostatic restoring coefficients due to wave passage

Table 7 :
Simulation initial conditions column is maximum roll amplitude for the simulations (max |φ sim |).The sixth is maximum roll amplitude for the experiments (max |φ exp |).The seventh and final column is the percentage error given by 100 max |φ sim | − max |φ exp | max |φ exp | ,

Table 9 :
Table 10 the added mass.Table 11 contains the hydro-dynamic damping parameters, while Table 12 contains the body motion parameters.Table 13 contains the wave motion parameters for heave.Table 14 contains the wave motion parameters for roll.Table 15 contains the wave motion parameters for pitch.Table 16 contains the external wave excitation parameters.Note that α 3 and α 5 are given in radians.Rigid body inertia

Table 10 :
Added mass