## Abstract

The current intact stability criteria for ships do not account for their dynamic behaviour in irregular waves sufficiently. We provide analytical criteria for ship and sea state parameters that indicate large roll amplitudes or capsizing and determine the mean time for these events. These results can be used to formulate vulnerability criteria for the occurrence of parametric roll. We use reduced order models in order to perform an analytical analysis to this complicated problem. Therefore, we study the dynamics of a ship in random seas as a perturbation of a Hamiltonian system. The random process across the energy levels of the Hamiltonian system is approximated by a Markov process, which is obtained by stochastic averaging. The noise owing to wave excitation is a non-white stationary stochastic process. This is an extension to previous results, where white noise excitation was used.

## 1. Introduction

Today’s modern vessels are designed towards efficiency with respect to cargo load and energy consumption. Stability criteria should lead to the safety of crew, cargo and environment. Because of the complex ship dynamics resulting from fluid–structure interaction in random seas, the current criteria for the evaluation of intact stability of ships developed by the International Maritime Organization (IMO) are based on ship restoring capabilities in still water (IMO 2009). There are no dynamic criteria for vulnerability to parametric roll besides operational recommendations yet, which is not appropriate when considering the amount of intact vessel losses (Kreuzer & Wendt 2000). The IMO intact stability code notes, however, the necessity of improved intact stability criteria including the dynamic behaviour in a realistic seaway. Researchers worldwide are working on the improvement of knowledge on the intact stability (France *et al.* 2003; Shin *et al.* 2004; Spyrou 2006; Belenky & Sevastianov 2007). Important results are summarized in the technical regulatory document American Bureau of Shipping (2004). Currently, second generation intact stability criteria are under development (cf. Peters *et al.* 2011) using a multi-tiered approach in three levels of complexity for the stability failures: pure loss of stability, parametric roll and surf-riding/broaching-to. The current definitions of the three levels of complexity are given in Peters *et al.* (2011) and . Level 1 vulnerability criteria should be conservative and given by a simple procedure, algorithm, or formula. Level 2 vulnerability criteria should use models with reduced complexity and can involve numerical calculations such as numerical integration and hydrostatic calculations. The last and computationally most expensive level 3 direct stability assessment has to involve state-of-the-art numerical tools. Typically, ships are travelling in head or quartering or sometimes following seas, in order to reduce the amplitude of roll angle. Beam sea is avoided if possible because large roll oscillations occur if waves come at an angle of 90^{°} measured from the ship’s forward direction. However, beam sea conditions cannot be excluded if the ship has manoeuvering problems. For vessels such as modern container ships and fast roll-on/roll-off ferries (RoRo ferries), roll motions owing to waves from the front or rear can appear suddenly in open seas, leading to dangerous situations up to capsizing. For a ship in unidirectional head or following waves, no direct excitation of roll motion occurs. But it is well known that parametric resonance is possible owing to oscillations of roll-restoring moment. A nonlinear model has to be used when analysing large-amplitude ship rolling. Theoretical and experimental studies of this issue were presented in Oh *et al.* (2000). Parametric roll excitation in head or following seas was already recognized by Froude (1861) and later experimentally verified (Kempf 1938; Paulling & Rosenberg 1959). The mechanisms of parametrically excited roll motion were studied by Grim (1952) on the basis of the Mathieu equation. By investigating ship motion in regular waves, it is found that parametrically excited roll occurs at a wave frequency that is twice the roll motion eigenfrequency. A realistic description of the seaway has to include random wave elevation because of unpredictable forces due to wind and current. This yields that probabilistic methods have to be used for the evaluation of ship dynamics in realistic seas. Arnold *et al.* (1986) computed the largest Lyapunov exponent for a parametrically excited stochastic linear oscillator, which arises when the harmonic parametric forcing in the Mathieu equation is replaced by a Gaussian stochastic process. A negative maximal Lyapunov exponent yields the almost-sure (a.s.) stability of the trivial solution of the oscillators. This means that the stability of roll motion can also be determined if parametric excitation owing to incident waves is modelled by a Gaussian stochastic process. Bulian (2008) extended the concept of the Grim effective wave (Grim 1961) to a travelling effective wave described by stochastic processes for its amplitude and phase. The essential task is to determine under which sea environment a given ship design is endangered. Or more precisely, under which circumstances extreme vessel dynamics occur more probably and can be reduced by changing the design and operational parameters.

The main purpose of this article is to obtain analytical results for the energy of the ship roll mode and use these results to formulate a first passage time formula using the most simple analytical formulation if taking into account excitation due to random seas. This formula could be used for the development of a probabilistic level 2 criterion for parametric roll. In order to achieve such a result, we reduce the dimension of the original process to a one-dimensional process by stochastic averaging. The second important statement of this paper is a background for vulnerability criterion for parametric roll based on the top Lyapunov exponent.

An outline of the paper is as follows: in §2, we describe the modelling of random ocean wave amplitude and phase and we provide additional background information and motivation for the problem that are needed to determine the random fluid–ship interaction forces. For the reader’s convenience, a summary of fluid–ship interaction forces is presented in §3. In §3, we also show under which circumstances the governing equations for ship roll dynamics can be modelled by a weakly perturbed softening spring-type Duffing oscillator with additional cubic damping (cf. Arnold *et al.* 2004). We formulate the ship roll dynamics problem as a stochastically perturbed weakly Hamiltonian system. Analysis for the hardening spring-type noisy Duffing–van der Pol oscillator was done in Sri Namachchivaya (1991) and Sri Namachchivaya *et al.* (2001), where the excitation was a standard Wiener process. A method of finding the top Lyapunov exponent and moment Lyapunov exponent, featuring an expansion of an invariant measure, is presented in §4. The results of Sri Namachchivaya (1991) and Sri Namachchivaya *et al.* (2001) were extended for a more realistic non-white stationary stochastic process in §5. In §7, the analytical results from §§5 and 6 are verified by numerical simulations. By this, it is shown that the stochastic averaging theory is applicable for the case of real ship roll dynamics. Applications of our results to related systems are discussed in §8 and the paper concludes in §9 with a summary of our results and some additional remarks.

## 2. Description of random ocean waves

A well-known model of random long-crested sea waves is given by superposition of infinitely many harmonic waves using linear wave theory. The amplitude of single waves with a given frequency depends on the underlying sea state given by the corresponding spectral density. Such an irregular long-crested wave surface can be written as
2.1with wavenumber *κ*(*ω*) and frequencies *ω* corresponding to a one-sided spectral density *S*(*ω*). We obtain an initially irregular wave surface, if a random phase shift *ζ*(*ω*) is added, which is uniformly distributed in the interval [0,2*π*). The above integral is not a Riemann integral but a summation rule over the frequencies *ω*. Note that *Z*(*x*,*t*) is evolving in space *x* and time *t*. Various spectral densities were defined to describe different sea states. Common sea spectral densities are the Pierson–Moskowitz (PM) spectrum for deep-water waves and the Joint North Sea Wave Project (JONSWAP) spectrum for shallow-water waves. Formulations of these spectra, depending on significant wave height and modal frequency, can be found in the electronic supplementary material. The stochastic process defined by *Z*(*x*,*t*) is too complicated to be used in second-level vulnerability criteria because it involves infinitely many components, and thus too many possibilities how waves can surround a ship hull in the proposed approach. Therefore, the random wave surface *Z*(*x*,*t*) is approximated by an effective wave surface
2.2where *L* is the length of the effective wave. The random processes *η*_{c}(*t*), *η*_{s}(*t*), *η*(*t*) and *ψ*(*t*) of *Z*_{eff} are related by
2.3
2.4The error between the random wave surface *Z* and the effective wave surface *Z*_{eff} is minimal, if *η*_{c} and *η*_{s} are solutions of the unconstrained quadratic optimization problem
2.5The optimal solution of (2.5) is obtained by solving for ∂*I*(*η*_{s},*η*_{c})/∂*η*_{s}=∂*I*(*η*_{s},*η*_{c})/∂*η*_{c}=0, which gives
2.6and
2.7Here, the transfer functions *f*_{s} and *f*_{c} are given by
2.8where *r*=(*L*/2)*κ*(*ω*). Details of the derivation of equations (2.6) and (2.7) are given in the electronic supplementary material. Comparison of the processes (2.6) and (2.7) with the irregular wave surface (2.1) shows, that the spectral densities *S*_{ηc} and *S*_{ηs} of the stochastic processes *η*_{c} and *η*_{s}, respectively, are related to the spectral density *S*(*ω*) of the random sea elevation process by
2.9The major advantage of the effective wave *Z*_{eff}(*x*,*t*) in comparison with the irregular wave surface *Z*(*x*,*t*) is that we have a model for the excitation owing to irregular seas depending only on two parameters *η*_{c} and *η*_{s} or *η* and *ψ*, respectively, which actually are stochastic processes with known spectral densities. If we set *ψ*=0 or equivalently *η*_{s}=0 in equation (2.2), we obtain the original Grim effective wave (Grim 1952). For most ship hulls, a wave of the same length as the ship length will cause the largest variations of the righting lever. Therefore, the effective wavelength is chosen equal to the ship length *L*. If the ship is moving with mean speed *U* in waves that approach from an angle *χ* with respect to the ship’s forward direction, it is convenient to describe the sea state in a ship fixed coordinate system. Then, the frequency of wave encounter *ω*_{e} in the ship fixed frame is given by
2.10and the spectrum with respect to the encounter frequency is
2.11where |d*ω*/d*ω*_{e}| is the Jacobian. In this case, the encounter spectrum *S*(*ω*_{e}) has to be used instead of *S*(*ω*) in equation (2.9). In deep water, *κ*(*ω*)=*ω*^{2}/*g* and thus |d*ω*/d*ω*_{e}| can be calculated as
2.12The effect of speed *U* and effective wave on a PM spectrum is shown in figure 2 in the electronic supplementary material.

## 3. Ship motion model

The motion of a ship and wave propagation are described with respect to the right-handed Cartesian coordinate system *O*(*x*,*y*,*z*), where *z* is pointing upwards out of the fluid domain. The corresponding ship fixed coordinate system is *O*′(*x*′,*y*′,*z*′). The (*x*,*y*) plane defines the free water surface at rest. The centre of gravity of the ship is located at the origin *O*′ of the ship fixed coordinate system, and the *x*′-axis coincides with its forward direction. The respective ship motions in *x*-, *y*- and *z*-directions are surge, sway and heave. Roll *Φ*, pitch *θ* and yaw *ζ* are the rotations around the *x*′-, *y*′- and *z*′-axes, respectively. The forces ** F**=(

*F*

_{x},

*F*

_{y},

*F*

_{z})

^{T}and moments

**=(**

*M**M*

_{Φ},

*M*

_{θ},

*M*

_{ζ})

^{T}on the hull are determined by integrating the fluid pressure

*p*over the wetted surface

*Ω*of the ship hull, 3.1Here,

**is the unit normal vector in the outward direction on the ship hull**

*n**Ω*and

**is the corresponding position vector. The pressure distribution due to incident waves at a depth**

*r**z*below the calm water surface is 3.2where

*ζ*

_{a}is the wave amplitude and

*d*is the water depth. Higher-order effects owing to ship generated waves are neglected. The restoring forces of a ship in waves can be determined by the righting lever function

*GZ*. First the centre of buoyancy

*B*has to be computed for relevant roll angles by integration of pressure over the wetted surface, which is obtained by computation of (3.1). The righting lever is the horizontal distance from

*Z*to the centre of gravity

*G*, and the metacentric height

*GM*is given by the distance from

*G*to the metacentrum

*M*, cf. figure 1. As can be seen in figure 2, the slope of righting lever

*GZ*with respect to the roll angle is higher for a wave trough than for a wave crest. This leads to increased initial stability of the ship’s upright position in the wave trough and reduced initial stability in the wave crest. Coupling between pitch and roll is negligible in beam sea conditions. The strongest parametric resonance can be found at an encounter frequency to damped roll eigenfrequency ratio of 1:2, i.e.

*ω*=2

*ω*

_{d}. If heave and pitch motions are large for parametric roll resonance conditions, one has to study the coupled heave–pitch–roll motion, which can be stated as 3.3

The functions *Z*(*η*(*t*),*ψ*(*t*)), *Θ*(*η*(*t*),*ψ*(*t*)) and *M*(*η*(*t*),*ψ*(*t*)) are force and moments due to wave excitation, respectively. *F*_{z} and *M*_{θ} are components of (3.1). Note that the added hydrodynamic masses *A*_{z}(*ω*_{e}), *A*_{θ}(*ω*_{e}) and *A*_{Φ}(*ω*_{e}) depend on the encounter frequency *ω*_{e}. For the heave–pitch–roll system, the righting lever is a nonlinear function of heave, pitch and roll with a large variety of nonlinear dynamics in resonant frequency regions. For the development of a closed-form formula appropriate for level 2 stability criteria, the case of small heave and pitch has to be considered. The analysis of coupled heave–pitch–roll motion in random seas can only be done by extended numerical computations and should be included in level 3 stability criteria. However, to account for small heave and pitch motion, we assume a quasi-static behaviour. This means that for incident waves, the ship is in static equilibrium with respect to heave and pitch. The righting lever function is computed for various roll angles *Φ* and for long-crested waves of the same length as the ship length with an elevation of *η*, cf. Dostal & Kreuzer (2011). The wave crest is located at in the ship coordinate system. Therefore, the rolling behaviour of a ship can be represented by the following equation if heave and pitch motions are small:
3.4Here, *I*_{xx} is the roll moment of inertia, *A*_{xx} is the hydrodynamic added mass evaluated at the natural frequency *ω*_{n}, *b*_{1} and *b*_{3} are linear and cubic damping coefficients, *g* is acceleration due to gravity, *Δ* is the displacement and *M* is a roll excitation moment, which is small, if the ship travels in about the same direction as the incident waves. The roll equation (3.4) is valid for a ship travelling in an arbitrary direction in long-crested random seas as described in §2. For conventional ships, the *GZ* function approaches a global maximum as the roll angle *Φ* increases from zero. As the roll angle increases further, the *GZ* function crosses zero from positive to negative at the angle of vanishing stability. In calm water, we denote this angle as *Φ*_{c}. For a ship in a fixed wave, the angle of vanishing stability will in general be lower than *Φ*_{c} in the case of a wave crest amidship and usually higher than *Φ*_{c} for a wave trough amidship, cf. figure 2. Considering only hydrostatic forces due to long-crested incident waves, a ship can be regarded as capsized if the maximal value of the angle of vanishing stability in waves is exceeded. Then, the ship will approach its capsized equilibrium. In the region between the maximal and minimal values of the angle of vanishing stability in waves, the ship can capsize or survive, which depends on the instantaneous wave crest position and roll angle. The leading nonlinearity of righting lever *GZ* is cubic in *Φ*, and the dependence of righting lever variations with respect to the wave crest position is harmonic. Hence, we approximate the righting lever *GZ* by
3.5Here, is substituted by *η*_{c} according to equation (2.4). The use of higher-order nonlinearities for the approximation of *GZ* does not decrease the approximation error significantly, as can be seen in Dostal & Kreuzer (2011). For the following, we set *ξ*_{t}:=*η*_{c}(*t*) and approximate *M* by
3.6After some obvious transformations, equation (3.4) with *GZ* and *M* given by (3.5) and (3.6), respectively, can be written as the system of equations
3.7From equations (2.9) and (2.11), the spectral density *S*_{ξt} of the stationary stochastic process *ξ*_{t} is given by
Using the rescalings
3.8we obtain
3.9We are interested in the energy of roll dynamics. Hence, we use the Hamiltonian formalism and rewrite equation (3.9) in terms of the total energy defined by the Hamilton function
3.10where *α*_{1},*α*_{3}>0. Then, (3.9) can be written as
3.11Making use of the fact that the roll damping and excitation forces due to fluid–structure interaction are small, we can assume *ε*≪1 to be also a small parameter. Thus, we have transformed the original equation for roll motion (3.4) to a weakly perturbed Hamiltonian system (3.11) excited by the stationary random process *ξ*_{tε}. The coefficients *ν*_{1} and *ν*_{2} are additive and parametric noise intensities, respectively. Contour lines of the Hamiltonian are shown in figure 3. The fixed points of system (3.11) without dissipation and random perturbation, i.e. *ε*=0, are
3.12

These saddle points *P*_{1} and *P*_{2} are connected by the heteroclinic orbit
3.13The time derivative of the Hamiltonian for system (3.11) is given by
3.14Combining the first equation of (3.11) and equation (3.14), and using
3.15obtained from (3.10), we obtain the system
3.16An important property of equation (3.16) is that the energy level *H* changes slowly compared to the oscillations of the variable *x*. This enables the application of stochastic averaging to this system in §5. But first, we state results for stability of the trivial solution of equation (3.4).

From now on, we omit the index *ε* in *t*_{ε} and use *t*:=*t*_{ε} as the scaled time.

## 4. Small-noise stability analysis

For the assessment of dangerous ship dynamics or even capsizing, it is necessary to determine for which sea state statistics and ship speeds large roll amplitudes occur. A first indication is the bifurcation behaviour of stability of the trivial solution in equation (3.4). By this means, response amplitudes can be excluded for specific significant wave heights and mean wave periods. For a ship in following or head sea conditions, no direct excitation of roll motion owing to incident waves is possible. But the roll motion can be excited by waves travelling along the ship, leading to a change of the submerged part of the hull, if roll velocity or roll angle is initially not zero. This results in variation of the righting lever in the time leading to parametric excitation of roll motion. The mechanisms of parametric excited roll motion were studied by Grim (1952) on the basis of the Mathieu equation
4.1It is common to plot the stability boundary of the trivial solution of (4.1) in the space of parameters *γ* and *ν*. Such a plot is known as an Ince–Strutt diagram. For the damped Mathieu equation
4.2the Ince–Strutt diagram is plotted in figure 4 for various values of damping *β* in the *γ*–*ν* plane. Here, *γ* is the squared damped natural frequency, which is given by *γ*=*α*−*β*^{2}. We are interested in the stability of the upright position for the case of parametric excited roll motion in random seas. As was also pointed out in Spyrou (2005), we have to investigate the a.s. stability of the trivial solution of equation (3.7) with *ν*_{1}=0. Thus, we have to linearize (3.7) and study the resulting problem
4.3where *ξ*_{t} is a stationary Gaussian process with spectral density *S*_{ξtξt}, as above.

### Remark 4.1

For determination of the linearized coefficients, the fitting of the righting lever curve should be applied for roll angles reaching from zero to 1^{°}. The reason is that here we are interested in the slope of the righting lever function *GZ* in the *Φ* direction, which corresponds to metacentric height (*GM*) variations with incident waves instead of *GZ* variations.

Equation (4.3) is the parametrically excited stochastic linear oscillator, which arises when the harmonic parametric forcing in the Mathieu equation (4.1) is replaced by a Gaussian stochastic process *ξ*_{t}. For this case, the top Lyapunov exponent was computed by Arnold *et al.* (1986).

### (a) Top Lyapunov exponent

The *Furstenberg–Khashminskii* formula yields the top Lyapunov exponent *λ* of system (4.3). A negative top Lyapunov exponent ensures the a.s. stability of the trivial solution of (4.3) and consequently the a.s. stability of the original nonlinear system. For the undamped stochastic oscillator given by
4.4we can use the following result.

### Theorem 4.2

(Arnold *et al.* 1986) *For the case* *and γ>0, the top Lyapunov exponent* *for the oscillator in (4.4) is given by
*4.5

Now, we consider the under-damped case *β*^{2}<*α*. Inserting
4.6in (4.3) yields
4.7Because of transformation (4.6), the Lyapunov exponents *λ* for (4.3) and for (4.4) are related by . Therefore, we can use the results stated in theorem 4.2 for equation (4.7) with *γ*=*α*−*β*^{2} and obtain
4.8With formula (4.8), the a.s. stability of the trivial solution can be computed for the stochastic process *ξ*_{t} with arbitrary spectral density *S*_{ξtξt}. From (4.8), the stability boundary is given by
4.9We can now formulate a simple criterion for the stability of a ship travelling in random head or following waves.

### Vulnerability criterion for parametric roll based on top Lyapunov exponent

*Let a ship travel in random longitudinal head or following waves and let the excitation process* *ξ*_{t} *due to random fluid forces be defined by the one-dimensional spectrum* *S*_{ξtξt}. *Then, the upright position of this ship is a.s. stable, if the Lyapunov exponent* *λ*, *as stated in equation (4.8 ), is smaller than zero, where* *α*, *β* *and* *ν* *are coefficients from the linearized roll equation (4.3)*.

As an example, we define the spectrum of *ξ*_{t} by
4.10with *a*_{1}=0.2, *a*_{2}=1 and *b*_{0}=1. These parameters lead to a modal frequency of *ω*_{m}=1. In analogy to the Ince–Strutt diagram for the deterministic equations (4.1) and (4.2), we are interested in the stability boundaries in the plane of stiffness *α* and noise intensity *ν* for different values of damping *β* in the under-damped case *β*^{2}<*α*. The under-damped case is the relevant case for the ship roll behaviour, since the roll damping is usually small in comparison with the restoring of upright position *α*. The noise intensity *ν* depends on the sea state as well as on the change of righting lever for different wave crest positions. Figure 5 shows the first instability region, where the main resonance is located at . Contrary to the deterministic case (4.2), the whole *α*–*ν* plane would represent instability if the damping *β* goes to zero. This result is obtained from equation (4.9), by neglecting higher-order terms in *ν*. For vanishing damping, it follows that level 1 vulnerability criteria based on the stability boundary (figure 4) of the damped Mathieu equation (4.2) are not conservative enough, if taking into account random excitation. Therefore, using the above vulnerability criterion for parametric roll based on the top Lyapunov exponent would improve the level of safety.

### (b) Moment Lyapunov exponent

We are also interested in the stability of the *p*th moment of the solution *x*(*t*;*x*_{0}) of equation (4.3), since the top Lyapunov exponent does not ensure the stability of moments. The exponential growth rate of the *p*th moment is given by the moment Lyapunov exponent, which is defined as
4.11Expansion in *p* yields
4.12It is difficult to obtain a closed-form solution for the moment Lyapunov exponent of equation (4.3). However, if the damping and noise intensity in (4.3) are small, a closed-form solution can be obtained. Therefore, consider
4.13with a small parameter *ε*≪1. The system (4.13) has eigenvalues
Then, we can apply the following theorem, which is proved in Arnold *et al.* (1997).

### Theorem 4.3

*For small noise, the pth moment Lyapunov exponent of system (4.13) is given by
*4.14

A condition ensuring the negative moment Lyapunov exponent *g*(*p*) is more stringent than ensuring a negative top Lyapunov exponent *λ*.

## 5. Stochastic averaging

We extend the present theory for averaging of equation (3.11). Up to now, asymptotic methods in *ε* were used for determining the averaged equations of system (3.16) subjected to real noise excitation (Roberts 1982*a*; Moshchuk *et al.* 1995; Roberts & Vasta 2000). For the case of beam sea, Roberts (1982*b*) considered Jacobian elliptic functions in order to obtain an averaged solution for the roll energy. In that approach, the Jacobian elliptic functions were approximated by Fourier series, and the final results were obtained by numerical integration. Stochastic averaging of (3.16) was done for the case of white noise excitation (Roberts 1978; Zhu 1982; Sri Namachchivaya 1991) and a white noise excited Duffing oscillator with a two-well potential (Sri Namachchivaya *et al.* 2001). In the following, we extend the results from Sri Namachchivaya *et al.* (2001) and Sri Namachchivaya (1991) for equation (3.11) under non-white noise excitation. It should be noted that equation (3.16) is the softening spring-type Duffing oscillator with linear and cubic damping, by which many physical systems can be approximated. Therefore, the following results are widely applicable.

### (a) White noise case

First, we state the results for averaging system (3.11) subjected to white noise excitation. Therefore, let , where *W*_{t} is a standard Wiener process and d*W*_{t} its increment. Then, using the Itô lemma, we obtain from system (3.11)
5.1Averaging (5.1) according to Khasminskii (1968), we obtain the one-dimensional Itô equation
5.2where
5.3and
5.4Here, *T*(*H*) is the period of one oscillation of the fast variable *x* in the absence of noise and damping, i.e. *ε*=0, starting at the energy level *H*. Additionally, we use
5.5For , the period is given by
5.6where
5.7and
5.8The limits of integration ±*b*(*H*) are the points where and the periodic orbit intersects the *x*-axis, i.e. *b*(*H*) is the maximum value of *x* for each energy level *H* and is given by
5.9The function *K*(*k*) is the complete elliptic integral of the first kind, cf. Byrd & Friedman (1954). The elliptic modulus *k* is given by
5.10The integrals appearing in the equations for drift (5.3) and diffusion (5.4) exist for . They can be computed in terms of complete elliptic integrals of the first and second kind, *K*(*k*) and *E*(*k*), respectively. Then, we obtain
5.11
5.12
5.13
5.14
5.15
5.16
and
5.17Thus, we have obtained a one-dimensional Itô equation (5.2) for the process of total energy *H*(*t*) of system (3.11) subjected to white noise excitation.

### (b) Real noise case

For practical problems, it is crucial to consider real noise excitation. For the non-white noise case, we need to determine the functions *x*(*t*+*τ*) and *y*(*t*+*τ*) in terms of *x*(*t*) and *y*(*t*) in order to apply stochastic averaging and obtain a closed-form solution. Therefore, an analytical solution of the differential equation (3.11) is needed, which can be obtained for *ε*=0 in terms of Jacobian elliptic functions. Then, addition formulae for Jacobian elliptic functions can be used to eliminate the time shift and obtain the states *x* and *y* as functions of time *t* only. In the case (non-capsized case), we obtain for *ε*=0 a solution of (3.11) by
5.185.19with . The expressions for *x*(⋅) and *y*(⋅) contain the Jacobian elliptic functions sn(⋅,*k*), cn(⋅,*k*) and *dn*(⋅,*k*), see Byrd & Friedman (1954). The elliptic modulus is the same as in the white noise case. From now on, we omit the elliptic modulus *k* and use the abbreviations
5.20If the subscript *τ* or *t*+*τ* is used, we refer to the argument *qτ* or *q*(*t*+*τ*), respectively. Applying stochastic averaging to system (3.16), we obtain a one-dimensional Itô stochastic differential equation for the energy level *H*. For this, we shall assume that the stochastic process *ξ*_{t} is stationary with mean zero and has autocorrelation function *R*_{ξtξt}(*τ*)=*E*{*ξ*_{τ+t}*ξ*_{t}}, which approaches zero sufficiently fast as *τ* increases, implying *ξ*_{t} is strong mixing as defined in the electronic supplementary material, cf. Papanicolaou & Kohler (1975).

### Proposition 5.1

*Let ξ*_{t} *in (3.16) be a zero-mean stationary stochastic process that is strong mixing and satisfies the mixing rate condition. Let H be the solution of system (3.16). Then, for fixed* , *the integrals*
5.21
and
5.22*exist, and the process H converges, as* *weakly on the time interval of order* *to a diffusion Markov process* *H* *satisfying the Itô stochastic differential equation*
5.23

Note that we use the variable *H* for the solution of the original system (3.16) and for the solution of the Itô equation (5.23), since in the following it is clear, where we refer to the averaged solution of *H* (i.e. solution of (5.23)).

### Proof.

The system (3.16) exhibits slow and fast oscillations, which are the energy and roll oscillations, respectively. Additionally, the driving stochastic process *ξ*_{t} fulfills sufficient mixing conditions. Therefore, we can apply the limit theorem as given in Papanicolaou & Kohler(1975) to system (3.16). In our case, the fast roll oscillation *x*(*t*) is averaged by using the averaging operator for periodic functions , which is given by
Substituting (5.18), (5.19) and d*t*=d*u*/*q*, we obtain for the drift *m*(*H*)
5.24In analogy, the diffusion *σ*(*H*) is given by
5.25Note that there is no additional term in equation (5.24), because the diffusive part of the first equation of system (3.16) is zero (cf. Papanicolaou & Kohler 1975; Sri Namachchivaya 1990). We obtain the expressions for *m*(*H*) and *σ*^{2}(*H*) as stated in proposition 5.1 after cancelling the terms which integrates to zero.

### Remark 5.2

The stochastic process *H* defined by the averaged equations (5.23) can be interpreted as the total energy variation of roll oscillation owing to wave forces and roll damping acting on the vessel. Each energy level *H* corresponds to specific oscillatory roll motion. The energy level *H* defines the roll angle and velocity by means of equation (3.10). Therefore, the roll oscillations become larger with increasing energy Level *H*.

Numerical evaluation of the integral in (5.21) is expensive because there is an integrable singularity at , which corresponds to *y*=0 in the original system (3.11). Fortunately, the integrals in equations (5.21) and (5.22) can all be solved in closed form in terms of complete elliptic integrals of the first, second and third kind. The explicit evaluation of the drift and diffusion terms is presented in appendix A*a*. This is achieved by using the addition formulas of Jacobian elliptic functions, which eliminates the time shift of the functions sn_{t+τ}, cn_{t+τ}, and dn_{t+τ} in equations (5.21) and (5.22).

It is not difficult to show that for *R*_{ξtξt}(*τ*)=*δ*(*τ*), where *δ*(⋅) is the Dirac delta, the drift *m*(*H*) and diffusion *σ*^{2}(*H*) from proposition 5.1 are given by (5.3) and (5.4), respectively. This is shown in the electronic supplementary material.

## 6. Probability density

A common probability measure is the probability density function (PDF). The differential generator of the averaged process *H* defined by the Itô equation (5.23) is given by
6.1We introduce the speed and scale measures, which transform the diffusion process (5.23) on its natural scale, where the drift is identically zero. This task is achieved via the transformation *S*(*x*) given by
is called the scale density. The reduced generator takes the form
6.2After defining the speed density *μ*(*y*)=1/(*σ*^{2}(*y*)*s*(*y*)), we can reduce to
6.3where *M* is the speed measure defined by
6.4A discussion on the meaning of scale density *s* and speed density *μ* can be found in Sri Namachchivaya (1989). The stationary PDF *p*_{st} associated with (5.23) is the solution of the Fokker–Planck equation
6.5and can be given in terms of
6.6The coefficients *c*_{1} and *c*_{2} are determined by the boundary and normality conditions. Following Karlin & Taylor (1981), the (left) boundaries can be classified as: entrance, if and ; reflecting, if , and *M*(*x*_{l})=0; exit, if . Here, *M*(⋅) is the speed measure, *Σ*_{l}(*x*_{l}) is the time to reach the left boundary *x*_{l} starting in *x*∈[*x*_{l},*x*_{r}], whereas *N*_{l}(*x*_{l}) is the time to reach *x*∈[*x*_{l},*x*_{r}] starting in *x*_{l}. These measures are given by
6.7If the righting lever curve *GZ* is of the softening type, then the angle of vanishing stability *Φ*_{c} exists and is an exit boundary (), since *Φ*_{c} can be reached in finite time. Moreover, due to negative roll restoring for |*Φ*|>|*Φ*_{c}|, sample trajectories cannot return to |*Φ*|<|*Φ*_{c}| if the initial energy *H*_{0} is smaller than *H*_{c}. This means that a capsized ship remains capsized and cannot return to an upright position any more. The energy corresponding to the angle *Φ*_{c} can be computed from (3.10) and (3.12). Since the system (3.16) has an exit boundary at *H*=*H*_{c}, a stationary density for (3.16) does not exist. If we assume a reflecting boundary at *H*_{c} (i.e. no capsizing), which is approximately the case for low-to-moderate noise intensities, then *c*_{1}=0. In this case, a stationary solution of the Fokker–Planck equation (6.5) exists and is given by
6.8From the stationary PDF of energy level *p*_{st}(*H*), we can compute the stationary PDF *p*_{st}(*b*) of maximal roll amplitudes *b*(*H*) corresponding to the respective energy levels *H*. Using the probability density transformation *p*_{st}(*b*)=*p*_{st}(*H*)((d/d*H*)*b*)^{−1}, we obtain
6.9

## 7. Simulation

In this section, we show, for two example cases, that the stationary PDF of maximal amplitudes *b* (6.9) is a good approximation to the time variant PDF of the original system (3.8) in the case of very small probability flow out of the region bounded by the heteroclinic orbit *γ* (see figure 3). Of course, a good approximation can be obtained only if additionally the assumptions of proposition 5.1 are satisfied. However, in the context of ship roll dynamics, it is important to validate that these assumptions are also valid for typical parameter sets of real ships and practically relevant sea states. For this, we use realistic data of a RoRo ferry (Kreuzer & Sichermann 2005) as stated in table 1, and compare the PDF of the original system (3.11) with the stationary PDF (6.9) of the averaged system (5.23). After computation of the righting lever approximation (3.5) for *χ*=0 rad, we obtain the coefficients for equations (3.7) and (3.9) as stated in table 2. In this section, we refer to *t* as the real time, and to *t*_{ε}:=*εt* as the scaled time. We use a continuous time autoregressive moving average process (CARMA processes) in order to generate the stationary zero-mean random process *ξ*_{t}=:*χ*_{1}(*t*). A sea process can be approximated accurately by a second-order CARMA processes, (cf. Dostal & Kreuzer 2011). The equations for the second-order CARMA process *χ*_{1} in continuous time are
7.1

where d*W*_{t} is the increment of a standard Wiener process. To generate the random process on the scaled time, we use the rescalings
leading to the system of filter equations
7.2Thereby we recall the fact that . The autocorrelation function of the second-order CARMA process *ξ*_{t} is given by
7.3

and the corresponding spectral density of *ξ*_{t} is
7.4

Because the theory of stochastic averaging yields an asymptotic statement for *ε*↦0, we can expect good results only for small *ε*. If we set *ε*=0.1, the parameters of the RoRo ferry are of the same order, as shown in table 2. This means that damping and noise intensity are small compared with the stiffness coefficients and , and stochastic averaging is applicable. Additionally, it is important that the random excitation process *ξ*_{t} satisfies the mixing conditions. For the PM spectra from table 3, the mixing conditions are fulfilled. We use the relation (5.9) between the energy level *H* and *b* and display our results in terms of *b*, which is the maximal value of *x* for each energy level *H*. It can be seen in figures 6 and 7, that the agreement between *p*_{st} and simulation is very good for the PM sea states from table 3. In particular, the accuracy of the solution for the PM sea state with *H*_{s}=14 m is remarkable. For the case of a JONSWAP sea state, with *H*_{s}=10 m given in table 3, there is also a good agreement between the probability density *p*_{st} of the averaged system and the probability density of the original system (3.9) obtained by simulation, as shown in figure 8. If the random excitation intensity is too high or if the spectral bandwidth of the excitation process is too short, then the assumption of small probability flow out of the region bounded by the heteroclinic orbit (3.13) is not valid anymore, and the approximation of probability density becomes less accurate. This can be seen in figure 9, where the random excitation is due to a JONSWAP sea state with *H*_{s}=13 m. As a consequence, we need additional theory at hand, to obtain reasonable results for cases where the probability density of energy level *H* cannot be approximately seen as stationary. Therefore, we compute mean first passage times for capsizing in the next section.

## 8. Mean first passage time

It is not economically feasible, or even unavoidable under storm conditions, to develop a ship design such that extreme motion can be ruled out completely. In rough open seas, there is always a positive probability that a vessel encounters very large waves. As a consequence, the occurrence of rogue waves results in a positive probability for capsize, if a ship is travelling a sufficient amount of time in rough open seas. But how long does it take on average for a specific ship to capsize or reach critical roll amplitudes that lead to shift of load or to danger for passengers and crew? To answer this question, it is necessary to determine the mean time for reaching certain limit roll amplitudes starting from lower initial roll amplitudes, i.e. solve a first passage time problem. The results are used as a measure for the security level of the ship design. The maximal roll amplitude for a specific energy level *H* is given by *b*, see equation (5.9). Therefore, we are now looking for the mean time it takes for the roll motion to reach certain energy levels *H*_{c} starting at an initial energy level *H*_{0}. Hence, we use the following result.

### Theorem 8.1

(Sri Namachchivaya 1989) *Let x*_{e} *be the entrance boundary (i.e. * *and* *and let x*_{c} *be the regular or exit boundary (i.e.* *then the mean time to reach x*_{c} *starting at any x∈[x*_{e}*,x*_{c}*] is given by
*8.1

In figures 10 and 11, the mean time (8.1) for capsizing is plotted as a function of the initial roll angle *Φ*_{0}:=*b*(*H*_{0}). Here, the roll motion can reach the critical energy at roll angle , where the ship capsizes within a few wave cycles. The case of a sea state defined by the PM spectrum is less dangerous for a vessel because this spectrum has a broader bandwidth than a JONSWAP spectrum with similar significant wave height and modal frequency. Therefore, large roll motion due to parametric resonance needs a longer time to develop. This yields that the mean capsizing time for the case of a PM sea state is much higher than for a corresponding JONSWAP sea state, which can be verified in figures 10 and 11. By looking at these figures, it is also interesting to see the dramatic decrease of capsizing time if changing the significant wave height from 10 to 12 for a JONSWAP sea state. Note that, since we have used *ε*=0.1 to obtain the first passage times *t*_{ε}, the real time in seconds is *t*=10*t*_{ε}. So the mean capsizing time for a PM sea state with *H*_{s}=10 m is about 13.3 years (see figure 10), and the mean capsizing time for a JONSWAP sea state with *H*_{s}=10 m is about 3 days (see figure 11), if starting with an initial roll angle *Φ*_{0}=0.05^{°}. For the studied RoRo ferry, the modal frequency of *ω*_{m}=0.64 rad s^{−1} and the ship speed were chosen such that the modal wave encounter frequency is close to the resonant frequency 2*ω*_{d}≈0.95 rad s^{−1} for parametric roll.

### First passage time for parametric roll

*The total mean first passage time* *T*_{mfp}(*H*_{c}) *until the occurrence of a critical energy level* *H*_{c} *under the desired operational conditions has to be at least* *Y* *years, to achieve a sufficient level of safety. This time is given by*
8.2*The mean first passage time* *u*_{e} *for reaching a critical energy level* *H*_{c} *starting at* *H*=0, *is defined in theorem 8.1. For a given ship design, the time* *u*_{e} *depends not only on* *H*_{c}, *but also on the excitation due to random seas, on ship speed* *U* *and on the wave encounter angle* *χ*. *This set of operational conditions is denoted as* . *The spectral density of random sea state* *S*(*ω*) *is a function of sea state parameters such as significant wave height* *H*_{s} *and modal frequency* *ω*_{m}. *From sea statistics of the desired operational region and from operational guidelines*, *is the probability of the set* .

The above first passage time formula for parametric roll could be useful for the development of a level 2 vulnerability criterion for parametric roll. The critical energy level *H*_{c} and *Y* have to be specified depending on the purpose of the ship design. For example, a container ship should not exceed *H*_{c} corresponding to a critical maximal roll angle *Φ*_{c} at which loaded containers will shift or turn over. For a passenger ship, an *H*_{c} corresponding to a critical maximal roll velocity could be the limiting design property.

## 9. Summary and conclusions

In this paper, we have obtained analytical results for the behaviour of the roll motion of ships in random seas. These results were used to state a *vulnerability criterion for parametric roll based on a top Lyapunov exponent* and a *first passage time formula for parametric roll*. It was shown that the developed theory yields accurate results for the analysis of large real ships. An approximation for the probability density in a finite-time interval was obtained, assuming reflecting boundary conditions for the corresponding Fokker–Planck equation. The approximate probability density can be used in the design process of a ship to optimize its stability. As the model for roll motion is a softening spring-type Duffing oscillator with additional cubic damping, which is excited by a non-white stochastic process, our results are also applicable to various engineering problems.

## Acknowledgements

The authors are indebted to the DFG (Deutsche Forschungsgemeinschaft/German Research Foundation) for funding the project under contract Kr 752/29-1. Part of this research was carried out while L.D. was visiting the Department of Aerospace Engineering of University of Illinois at Urbana-Champaign, which was funded by NSF grant no. CMMI 07-58569. Any opinions, findings and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation.

## Appendix A

(*a*) *Real noise*

The integrals appearing in the equations for drift (5.21) and diffusion (5.22) can be computed in terms of complete elliptic integrals *K*(*k*),*E*(*k*) and *Π*(⋅,*k*) of the first, second and third kind, respectively. Here, and in the following, we use the abbreviations *s*:=sn(*qτ*), *c*:=cn(*qτ*) and *d*:=dn(*qτ*). The last integral in equation (5.21) was already computed in §5*a* by
A1The integral *I*_{m} in (5.21) is given by
A2Using the transformation , a solution of the above integral can be obtained for as
A3Evaluating the above integral and simplifying the resulting expression, we obtain
A4where
A5
A6and
A7If or if , then and . Then, we obtain the limit
A8In analogy, a solution of the integral in (5.22) can be obtained by
A9with
A10and
A11If , then and we obtain the limits
A12and
A13

- Received April 28, 2012.
- Accepted August 31, 2012.

- This journal is © 2012 The Royal Society