Extreme response prediction for nonlinear floating offshore structures by Monte Carlo simulation

Applied Ocean Research 29 (2007) 221–230 www.elsevier.com/locate/apor

Extreme response prediction for nonlinear floating offshore structures by Monte Carlo simulation A. Naess a,b,∗ , O. Gaidai a , P.S. Teigen c a Centre for Ships and Ocean Structures, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway b Department of Mathematical Sciences, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway c StatoilHydro Research Centre, NO-7005 Trondheim, Norway

Received 13 April 2007; received in revised form 8 October 2007; accepted 24 December 2007 Available online 8 February 2008

Abstract The paper describes a method for the prediction of extreme response statistics of floating offshore structures subjected to random seas by Monte Carlo simulation. The particular case of the horizontal surge motions of a tension leg platform is considered, taking into account both the first order, wave frequency and the second order, slow-drift motions. The advantage of the Monte Carlo method is its simplicity and versatility, which allows us to account for the effect of time-variant wave-drift damping, as well as nonlinear mooring characteristics without noticeable increase in the computational complexity. It is demonstrated in this paper that the commonly assumed obstacle against using the Monte Carlo method for estimating extreme responses, i.e. excessive CPU time, can be circumvented, bringing the computation time down to quite acceptable levels. c 2008 Elsevier Ltd. All rights reserved.

Keywords: Monte Carlo simulation; Mean crossing rate; Extreme response prediction; Tension leg platform

1. Introduction The problem of predicting the extreme response statistics of a floating offshore structure subjected to random waves by Monte Carlo simulation is the challenge we address in this paper. Specifically, we shall discuss the wave forces and motion response of a tension leg platform (TLP), with particular focus on the estimation of the total extreme surge response in a random sea way. The state-of-the-art formulation of the surge response of a TLP in random seas is in terms of a second order stochastic Volterra series. Such models and their statistics have been studied over several decades [1–10]. In the case of a linear dynamic system there is a possibility of accurately estimating the extreme response statistics [11–13]. In these papers a numerical calculation technique based on characteristic functions and the so-called saddle point method is developed. If the dynamic system is nonlinear, the saddle point method is not applicable. In [14] a numerical path integration method was elaborated for a dynamic system with ∗ Corresponding address: CeSOS/NTNU, Otto Nielsens vei 10, NO-7491 Trondheim, Norway. E-mail address: [email protected] (A. Naess).

c 2008 Elsevier Ltd. All rights reserved. 0141-1187/$ - see front matter doi:10.1016/j.apor.2007.12.001

nonlinear damping and nonlinear stiffness, but only the slowdrift response could be calculated. In principle, such limitations do not exist for the Monte Carlo method. However, the standard Monte Carlo method has not been very efficient for estimating extreme response statistics. Accurate estimates of events with a very low probability of occurring, which is typical for extreme responses, require very long simulation times. We shall show in this paper that this obstacle can be effectively circumvented by adopting an appropriate extrapolation technique obtained by exploiting the properties of the mean upcrossing rate function, which, under nonrestrictive assumptions, completely determines the extreme value distribution. It will be shown by numerical examples that by combining empirical estimates of the mean upcrossing rate obtained from the standard Monte Carlo method with an appropriate linear extrapolation technique, a dramatic reduction of the CPU time needed to obtain accurate estimates of the extreme response distribution can be achieved. In fact, this reduction in the required CPU time makes the standard Monte Carlo method available in practice for the accurate estimation of extreme responses with a modern laptop computer.

222

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

Over the years, several simulation based approaches for the estimation of extreme response of nonlinear floating structures subjected to random waves have been proposed. For a discussion, the reader is referred to [15]. The work described in the present paper, which follows a different approach from the methods described in [15], extends some previous work by the authors reported in [16]. 2. The dynamic model The equations of motion for a floating, rigid-body structure subjected to environmental forces like wind, waves and current would generally be written as ¨ ˙ M Z(t) + H(Z(t), Z(t), t) = F(t),

(1)

Here M denotes a generalized 6 × 6 mass matrix, Z = (Z 1 , . . . , Z 6 )T = the structure’s response vector, while H denotes a nonlinear vector function. F(t) denotes a stochastic loading process, which in general also depends on the response of the structure. This point will be discussed at some length below. The methodology developed in this paper applies directly to a dynamic model as expressed in Eq. (1). The only requirement is that the response time histories can be simulated. Since the main purpose of this paper is to illustrate the versatility and accuracy of the proposed method, we have chosen to discuss a simplified single-degree-of-freedom (SDOF) model for the surge response of a tension leg platform (TLP) in random waves. Except for interaction effects between different motion modes, the SDOF model allows for the introduction of most of the relevant nonlinear effects that may influence the response of the TLP. Hence, in this paper the following SDOF equation of motion will be studied M Z¨ (t) + H (Z (t), Z˙ (t), t) = F(t),

(2)

where Z = Z (t) denotes the surge response of the TLP, M is the mass of the platform, including added mass, and H is a nonlinear function to be specified. The hydrodynamic loading process F(t) will be assumed to consist of two components: A linear, first order (wave frequency) term F1 (t), and a nonlinear, second order term F2 (t). The linear, first order wave excitation in short-crested random seas in deep water can be expressed as a stochastic process F1 (t) as follows Z ∞Z π k1 (τ ; β) dβ X (t − τ ) dτ (3) F1 (t) = 0

−π

where dβ X (t) denotes the incremental directional components adding up to give the short-crested random seas X (t), which represent the R πwave elevation at some reference point. Formally, X (t) = −π dβ X (t), cf. [6]. X (t) is assumed to be a stationary zero-mean Gaussian process describing a short-term sea state. k1 (τ ; β) denotes the impulse response function of the linear transfer mechanism between a long-crested wave field propagating in the direction β and the surge excitation forces.

Similarly, the nonlinear, second order wave excitation forces in surge can be expressed as Z ∞Z ∞Z π Z π k2 (τ, τ 0 ; β, β 0 ) dβ X (t − τ ) F2 (t) = 0

0

−π −π 0

× dβ 0 X (t − τ ) dτ dτ 0

(4)

where k2 (τ, τ 0 ; β, β 0 ) is called the bidirectional quadratic impulse response function for surge. The total surge excitation forces F(t) on the structure that are commonly accounted for is then F(t) = F1 (t) + F2 (t), which is a second order stochastic Volterra series model. This representation gives a description of the hydrodynamic forces accurate to second order on a large volume structure, that is, a case where inertia forces dominate and the potential theory is applicable. However, there is one aspect of the second order forces on a floating structure that has turned out to be of some significance. This is related to what is known as wave-drift damping, which can be attributed to a dependence of F2 (t) on the motion of the TLP. A discussion of this is given below. To be able to express the force process F(t) through the random directional sea way as expressed in Eqs. (3) and (4), it is required to provide the relevant hydrodynamic transfer functions for the surge mode. As already stated, for simplicity of exposition, coupling with other motion modes is neglected. However, the method we present in this paper is equally suitable for a fully coupled analysis. Transfer functions are calculated by using commercial software, e.g. WAMIT. This program was used to calculate the (directional) linear and the (bidirectional) quadratic (surge force) transfer functions for the TLP. These are denoted by Kˆ 1 (ω; β) and Kˆ 2 (ω, ω0 ; β, β 0 ), respectively. These two functions are calculated on a discrete grid in parameter space. Specifically, let 0 ≤ ω0 < · · · < ω L+1 be the chosen discretization of the (positive) part of the frequency axis, and −π ≤ β0 < · · · < β M+1 ≤ π the corresponding discretization of the direction angle range (−π, π ). The spectral representation of the short-crested random wave elevation X (t) at the chosen reference location is then expressed as follows X (t) =

L X M X j=−L k=1

r

1 S X (|ω j |; βk )∆ω j ∆βk B jk eiω j t . 2

(5)

Here S X (ω j ; βk ) denotes the directional one-sided spectral density of the short-crested random seas. ∆ω j = (ω j+1 − ω j−1 )/2 and ∆β j = (β j+1 − β j−1 )/2. ω− j = −ω j . {B jk }, j = 1, . . . , L, k = 1, . . . , M, is a set of independent, complex Gaussian N(0,1)-variables with independent and identically distributed real and imaginary parts. These variables also satisfy the relations B− j,k = B ∗jk , where * signifies complex conjugation. This provides the appropriate extension to negative frequency indices. Also note that throughout this paper, when the summation index runs from negative to positive values, it invariably omits zero. As the focus of the present paper is on the first order and second order, slowly-varying response, it is appropriate to

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

introduce the slow-drift approximation. This amounts to putting Kˆ 2 (ω, −ω0 ; β, β 0 ) = 0

for ω · ω0 < 0.

(6)

This condition clearly implies that all sum-frequency components of the second order excitation forces are neglected, retaining only the difference frequency components, which are the only terms that can give rise to the slowly-varying response. In order to set up the proper dynamic model for the surge response Z (t), it is necessary to take into account the fact that hydrodynamic loading on a floating body depends on the motions of the body. In case of the slow-drift motions of the TLP, it has turned out that it is of some importance to take into account the dependence of the slow-drift force F2 (t) on the slowly-varying surge velocity Z˙ 2 (t). Since we have a nonlinear dynamic model, we shall have to introduce a definition of the slow-drift response Z 2 (t). A suitable definition would seem to be the following: The slow-drift response Z 2 (t) is obtained from the total response Z (t) by a low-pass filter that removes all wave frequency components. In practice, this may be achieved by a running mean operator used iteratively. To account for the dependence of the slow-drift force on the slowly-varying velocity, it is appropriate to write F2 (t, Z˙ 2 (t)) rather than F2 (t). However, F2 (t, Z˙ 2 (t)) is not directly available to us, but only F2 (t) ≡ F2 (t, 0) given by Eq. (4). Since, in the context of slow-drift motions, Z˙ 2 (t) is small, the following approximation is adopted, F2 (t, Z˙ 2 (t)) ≈ F2 (t, 0) +

∂ F2 (t, 0) ˙ Z 2 (t). ∂ Z˙ 2

(7)

≈ It is shown in [14] that for a TLP structure ∂ F∂2Z(t,0) ˙2 −c F2 (t, 0) ≡ −c F2 (t) for a suitable constant c > 0, may to some extent serve as a useful approximation to capture qualitatively the time-variant damping effect, which is the result of the expansion in Eq. (7). The first dynamic model adopted for Z (t) is now the following, M˜ Z¨ + C Z˙ + K (Z + ε Z 3 ) = F1 (t) + F2 (t, Z˙ 2 ) ≈ F(t) − cF2 (t) Z˙ 2 .

(8)

Here M˜ = M + m, ˜ where m˜ is an appropriately chosen (constant) added mass. C, K and ε are suitably chosen positive constants. This equation is rewritten in the equivalent form F(t) Z¨ + 2ωe ζ Z˙ + 2ωe cF ˜ 2 (t) Z˙ 2 + ωe2 (Z + ε Z 3 ) = , M˜

(9)

˜ ζ = C/(2ωe M), ˜ c˜ = c/(2ωe M). ˜ where ωe2 = K / M, Thus, the dynamic system is nonlinear with a Duffing-type hardening stiffness nonlinearity and time-varying damping. Since for the TLP the relative damping coefficient ζ is usually small, the contribution from the time-varying term is nonnegligible, especially for severe seas for which the slowdrift response is significant. The third order term in the restoring force, generally referred to as the “set down” effect, is caused by the fact that the tethers will induce the TLP to act like an inverted pendulum. To lowest order, the surge restoring force

223

is S sin θ = S ZL , where S is the total tension in the tethers, θ is the angular displacement corresponding to a surge displacement Z , and L is the effective length of the tethers. However, due to the pendulum effect, a surge motion will also involve a vertical displacement by an amount L(1 − cos θ ). Given that the pontoons are fully submerged, this will lead to an additional restoring force ∆S, so that the total surge restoring force equals (S + ∆S) sin θ = S sin θ + ρπ D 2 L(1 − cos θ ) sin θ ≈

S ρπ D 2 3 Z+ Z . L 2L 2

(10)

The above equation gives a reasonable representation of the surge restoring force to third order in Z provided tether elongation can be ignored and the tethers can be modelled as straight lines. This is, however, not quite correct and additional modifications of the restoring term would result by considering both tether elongation, tether dynamics and weight/buoyancy effects. Note that the set-down effect will have an influence also on the hydrodynamic loading process, which becomes dependent not only on Z˙ but on Z as well. Even if this dependence could have been taken into account, it has been neglected in this paper as it is of minor importance. For a TLP subjected to wave action only (i.e. in the absence of wind and current), Eq. (9) can be expanded further in order to account for other important nonlinear effects. The most important further extension of the equation of motion can be obtained by including viscous drag on the hull. The most basic way to accommodate a drag term is to include an expression of the type FD = 0.5ρ AC D | Z˙ | Z˙ ,

(11)

where A denotes the projected area of the TLP on a plane normal to the wave direction, ρ is the density of the sea water, and C D the nondimensional drag coefficient. This expression has the apparent weakness that it only accounts for the velocity of the hull relative to still water, and only indirectly (through the drag coefficient) accounts for shielding effects and diffraction. A more appropriate force term would result if one includes the relative velocity between the hull and the fluid and accounts for diffraction, that is, X FD = 0.5ρ Ai C D,i | Z˙ − Ui |( Z˙ − Ui ), (12) i

where the summation index runs over different component members of the submerged part of the TLP, and where the fluid velocities Ui in the surge direction (x-direction) are computed at the location of the structural component, and in the absence of that particular component. Such a procedure has been attempted with reasonable success in the past [17] for a gravity base structure. However, this approach is cumbersome and requires a lot of computations. Furthermore, for a compliant structure the accuracy may be questionable. The route taken here is therefore a more pragmatic one, which includes the undisturbed fluid velocity field only (i.e. the fluid velocity field in the absence of the structure) by employing a suitable global drag coefficient.

224

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

The drag coefficient can be estimated on the basis of model tests or obtained by superposition of the drag force on individual members. In this case shielding and diffraction can only be incorporated in an approximate sense. An argument for introducing this simplification is that for the heavy sea states considered in this paper, it is a reasonable approximation to assume that the wavelengths are “long” relative to the dimensions of the TLP. In such a case, the phase differences in the velocity field over the TLP structure may be effectively ignored. Furthermore, the cross flow principle can be assumed to be valid for the individual structural members of the TLP (i.e. the vertical TLP columns and the horizontal pontoons). For the one-dimensional case considered in the present paper the “drag-summation” therefore runs over the columns and the fore and aft pontoons, and include the x-velocities in the fluid domain only. Strictly speaking, the latter choice is not fully consistent with the cross flow principle as the x-component of the drag to (Z − q force on the pontoons would be proportional q

Ux ) (Z − Ux )2 + Uz2 , and to (Z − Ux ) (Z − Ux )2 + U y2 for the columns. However, the main uncertainty here lies with the actual magnitude of the drag coefficients, and discarding the x yz coupling may therefore seem like a small sacrifice, compared to the reduction achieved in the computational work. Allowing for a certain amount of linear damping the proposed equation of motion then takes the form: X M˜ Z¨ + C Z˙ + cF2 (t) Z˙ 2 + C˜ D, j | Z˙ − U j |( Z˙ − U j ) j

+ K (Z + ε Z ) = F(t), 3

(13)

where C˜ D, j now are dimensional drag terms, proportional to the total projected surface area of the structural members lumped together into group no. j, and where U j = Ux, j denotes the water particle velocity in the surge direction at a suitable reference point in group no. j. 3. Extreme value prediction There is in general no simple procedure for estimating the exact statistical distribution of the extreme response values. However, for the applications we have in mind in this paper, it is possible to obtain a very good approximation, which can be expressed in terms of the so-called mean upcrossing rate. To explain how this is achieved, let N + (ξ ; t1 , t2 ) denote the random number of times that the process Z (t) upcrosses the level ξ during the time interval (t1 , t2 ). By introducing the mean rate of ξ -upcrossings of Z (t) at time t, denoted by ν + (ξ ; t), then under suitable regularity conditions on the response process, which can be adopted here, the following formula is obtained E[N + (ξ ; t − ∆t/2, t + ∆t/2)] ∆t ∆t→0 Z ∞ = s f Z (t) Z˙ (t) (ξ, s) ds

ν + (ξ ; t) = lim

(14)

0

where f Z (t) Z˙ (t) (·, ·) denotes the joint PDF of Z (t) and Z˙ (t) = d Z (t)/dt, and E[·] denotes expectation. Eq. (14) is often

referred to as the Rice formula [18]. ν + (ξ ; t) is assumed throughout to be finite. Let M(T ) = max{Z (t) : 0 ≤ t ≤ T } denote the extreme value of the surge response process Z (t) over the long-term time interval of length T . The cumulative distribution function of M(T ) under the assumption that the stream of upcrossing events constitute a Poisson process, is then given in terms of the mean upcrossing rate as follows [19]. Prob(M(T ) ≤ ξ ) = Prob(N + (ξ ; 0, T ) = 0)  = exp −E[N + (ξ ; 0, T )]   Z T + ν (ξ ; t) dt , = exp −

(15)

0

where ν + (ξ ; t) denotes the mean level upcrossing rate at time t. For practical purposes, this is rewritten as [19,20]   Z + Prob(M(T ) ≤ ξ ) = exp −T ν (ξ ; w) f W (w) dw , (16) w

where f W (w) denotes the long-term (ergodic) PDF of relevant sea state parameters w, and ν + (ξ ; w) equals the mean upcrossing rate for the short-term (stationary) condition specified by w. In our case, W = (Hs , T p ), where Hs denotes the significant wave height and T p is the spectral peak period of a short-term sea state. This implies that the long-term PDF can be estimated from the scatter diagram of the sea states at the specified location. The practical long-term analysis would typically be carried out by weighting the crossing rate function over a selected group of sea states that is believed to provide a sufficient resolution over the scatter diagram. For the extreme responses, small to moderate seas do not contribute noticeably, despite their dominating rate of occurrence. Therefore, it is the large to severe sea states that are of particular importance, reducing substantially the number of sea states that need to be included. Eq. (16) can then be written as ! X + Prob(M(T ) ≤ ξ ) = exp − ck ν (ξ, k) T , (17) k∈K

where K denotes the set of included sea states parameterized by the index k. ν + (ξ, k) denotes the mean rate of upcrossings of the level ξ by the response process X k (t) for sea state k, and ck is its relative frequency of occurrence given by the scatter diagram. Time T may be e.g. a 20 year period, expressed in seconds T ≈ 6.3 × 108 s. 4. Empirical estimation of the mean crossing rate It is realized from the previous section that a key to providing estimates of the extreme values of the response process Z (t) on the basis of simulated response time histories, lies in the estimation of the mean upcrossing rate function for each shortterm sea state. For the simplicity of notation, in the following we suppress the sea state parameter k. By assuming the requisite ergodic properties of the short-term response process,

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

the mean upcrossing rate is conveniently estimated from the ergodic mean value. That is, it can be assumed that n + (ξ ; 0, T ) , T →∞ T

ν + (ξ ) = lim

(18)

where n + (ξ ; 0, T ) denotes a realization of N + (ξ ; 0, T ), that is, n + (ξ ; 0, T ) denotes the counted number of upcrossings during time T from a particular simulated time history. It is assumed that initial transient effects have been removed from these time series, that is, only the stationary parts of the time series are included. In practice, m time histories of a suitable, specified length, T0 say, are provided by simulation. The appropriate ergodic mean value estimate of ν + (ξ ) is then νˆ + (ξ ) =

m 1 X n + (ξ ; 0, T0 ), m T0 j=1 j

(19)

where n +j (ξ ; 0, T0 ) denotes the counted number of upcrossings of the level ξ by time history no. j. This will be the approach for the estimation of the mean upcrossing rate adopted in this paper. For a suitable number m, e.g. m ≥ 20 − 30 [21], and provided that T0 is sufficiently large, a fair approximation of the 95 % confidence interval CI0.95 (ξ ) for the value ν + (ξ ) can be obtained as √ CI0.95 (ξ ) = νˆ + (ξ ) ± 1.96 sˆ (ξ )/ m, (20) where the empirical standard deviation sˆ (ξ ) is given as !2 m n +j (ξ ; 0, T0 ) 1 X + 2 − νˆ (ξ ) . sˆ (ξ ) = m − 1 j=1 T0

(21)

It is pointed out that the choice of m and T0 for a specific case is mainly determined by an optimization of the estimate given by Eq. (21). Hence, the initially simulated time series may be divided into shorter blocks for the sake of optimization of this estimate. Note that the point estimate νˆ + (ξ ) is not affected by this operation. The consistency of the estimates obtained by Eq. (21) can be checked by the observation that Var[N + (ξ ; 0, T )] = ν + (ξ ) T since N + (ξ ; 0, T ) is a Poisson random variable by assumption; here Var[·] denotes variance. This leads to the equation s(ξ )2 /m = ν + (ξ )/(m T0 ). The empirical counterpart of Eq. (21) if the assumption that N + (ξ ; 0, T ) is a Poisson random variable is explicitly used, is then sˆ (ξ )2 /m ≈ νˆ + (ξ )/(m T0 ). It is worth pointing out that this last estimate is invariant with respect to blocking of the data, hence the optimization step is avoided. The approach to extreme value prediction presented in this paper derives from the realization that for the dynamic models of the type considered, the short-term as well as the long-term mean ξ -upcrossing rate as a function of the level ξ is in general highly regular in a specific way. In fact, it has been observed that for a wide range of dynamical systems the mean upcrossing rate tail behaves with good approximation like exp{−a(ξ −b)c } (ξ ≥ ξ0 ) where a, b and c are suitable constants. Hence, as discussed in detail in [22], it may be assumed that ν + (ξ ) ≈ q(ξ ) exp{−a(ξ − b)c } , ξ ≥ ξ0 ,

(22)

225

where the function q(ξ ) is slowly varying compared with the exponential function exp{−a(ξ − b)c } for tail values of ξ . Therefore, under the assumptions made, log | log(ν + (ξ )/ q(ξ ))| plotted versus log(ξ − b), exhibits an almost perfectly linear tail behaviour. If the function q(ξ ) could be replaced by a constant value, q say, which could be determined, we would immediately be in a position to apply a linear extrapolation strategy to find the values of ν + (ξ ) at extreme levels of ξ . In general, q(ξ ) is not constant, but its variation in the tail region is usually sufficiently slow to allow for its replacement by a constant, cf. [22]. In practice, an initial choice of the parameter q is extracted from an average of the ratio ν + (ξ )/ f Z (ξ ) for large values of ξ [22]; here f Z (ξ ) denotes the response probability density function, which would then also be estimated empirically. Having chosen a value for q, the value of the parameter b must be based on optimizing the linear fit, while at the same time keeping ξ0 as small as possible. The last proviso comes from the desire to obtain accurate extrapolation, which is related to the accuracy in the slope of the fitted straight line. Hence the region over which the fitting is done, should be as long as the linear fit allows. As mentioned above, and discussed in [22], the ratio ν + (ξ )/ f Z (ξ ) may vary with ξ . The initial choice of a value for the parameter q would consequently be somewhat arbitrary. It is therefore suggested to allow for optimization also with respect to the value of q, if that is deemed necessary. However, in this paper q is taken as an average of the ratio ν + (ξ )/ f Z (ξ ), empirically determined from the data tail, which seems to work well in most cases. The specific behaviour of ν + (ξ ) expressed by Eq. (22) is based on the underlying assumption that the appropriate asymptotic extreme value distribution for the stationary response data under study is the Gumbel distribution. For all relevant dynamic models for offshore structures as given by Eq. (1), with F(t) modelled as a stationary process, this would be the case. Therefore, the following asymptotic relation for the extreme value distribution is obtained [23], Prob(M(T ) ≤ ξ ) ∼ exp {− exp (−aT (ξ − bT ))} ,

T → ∞, (23)

for suitable parameters aT > 0 and bT . In practice, a commonly used approach is to assume that the observed extreme value response data do follow a Gumbel distribution. The problem with this approach is that classical extreme value theory cannot be used to decide to what extent the asymptotic distribution is actually valid for a given set of extreme value data. Note that the asymptotic Gumbel distribution corresponds to an asymptotic upcrossing rate that is purely exponential, that is, with c = 1.0 in Eq. (22). Hence, by adopting a much more general class of functions, with the purely exponential functions as a sub-class, to represent the upcrossing rate, as we have done here, the ability to capture subasymptotic behaviour is greatly enhanced. By this, the necessity to adopt a strictly asymptotic extreme value distribution, of questionable validity, is avoided. The extrapolation technique described above allows the empirical estimation of the ξ -upcrossing rate needed for a

226

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

Table 1 Particulars of the TLP Column diameter D (m) Eigenperiod surge Te (s) Relative damping ζ Total mass (incl. added mass) M˜ (kg)

10.0 128.8 0.05 1.5×107

specific extreme value estimation to be carried out with sufficient accuracy for most practical prediction purposes with much less numerical efforts than a direct estimation of the pertinent extreme values. Note that on the premise of a linear tail behaviour relative to a specified plotting procedure, the estimation uncertainty is related to the estimation of a slope. Hence, if the confidence band is sufficiently narrow over a suitable range of linear behaviour, the uncertainty in the slope will tend to be correspondingly small, and consequently, so will the uncertainty of the estimated upcrossing rate at the extreme response levels. This is a very important feature of the proposed extreme value prediction method, since the assumption that linear tail behaviour can be achieved by the proposed method seems to have wide applicability [22,24]. However, in practice, for short simulations the available data may not always allow for a very narrow confidence band, as will be illustrated in the next section. In [25] a general approach to the estimation of a confidence interval for the extrapolated estimate of the upcrossing rate at the extreme response level is proposed. 5. Numerical results In this section empirical estimates of the mean upcrossing rate obtained by Monte Carlo simulations combined with an optimized extrapolation procedure will be presented. A particular model of an offshore tension leg platform (TLP) is considered and the corresponding linear transfer function (LTF) and quadratic transfer function (QTF) are computed using a second order diffraction program (WAMIT). For simplicity, unidirectional seas are used, meaning that the directional argument β is skipped. This simplification should have no effect on the conclusions based on the comparison of accuracy. The combined first order and second order slowlyvarying surge deck motion is studied applying a single-degreeof-freedom model. The TLP particulars are detailed in Table 1 and the sub-surface part of the structure is shown in Fig. 1. The values in Table 1 are used to obtain the second order response. This means that for the second order response we have used a simplified version of Eq. (8), where mass ˜ stiffness K and damping coefficient C are frequency M, independent, which is a good approximation for the slow-drift motion. The time-invariant damping part ζ is considered to be 5%. In this paper, two versions of Eq. (9) will be used. The first version is a linear, time-invariant model obtained by putting c˜ = ε = 0. The second version is the fully nonlinear model where the parameter c˜ in (9) is chosen such that Var[2ωe cF ˜ 2 (t) Z˙ 2 (t)] ˙ is about 10% of Var[2ωe ζ Z (t)]. The parameter ε is estimated from the condition that 0.2Z (t) ≥ ε Z 3 (t) when Z (t) ≤ 6σ Z , i.e. even in the extreme response region, stiffness hardening

Fig. 1. Sketch of submerged part of TLP. Units in meters.

contributes not more than 20% relative to the linear part for severe sea, which leads to ε = 1.36 × 10−4 . Finally, the following approximate values were found, c˜ = 30/( M˜ g), g = 9.81 m/s2 for moderate seas, and c˜ = 90/( M˜ g) for severe seas. The adopted parameter values are largely arbitrary, but the choices made seem to provide a reasonable model for the chosen TLP structure. To get an accurate representation of the response process, there is a specific requirement that must be observed. Since the damping ratio is only 5%, the frequency resolution ∆ω must secure a sufficient number of frequency values over the resonance peak. This will ensure that the second order difference frequency response component captures the TLP surge dynamics with sufficient accuracy. It is commonly required that there are at least 5 discrete frequencies over the 2 is equal to or higher than half ˆ frequency range where | L(ω)| 2 ), where ˆ of the resonance peak height max(| L(ω)| −1  ˆ L(ω) = −ω2 + 2iζ ωe ω + ωe2 .

(24)

For the surge force QTF Kˆ 2 (ω, ω0 ) a suitable initial frequency grid must be chosen for which the values of the force QTF are calculated. The calculation of the force QTF is generally by far the most time consuming part of the numerical analysis. Therefore the initial grid is usually rather coarse to avoid excessive computer time. For the calculations at hand the discrete frequency range was the following: ω1 = 2π/30.0, . . . , ωn = 2π/4.0 (rad/s), n = 30. This necessitates the use of an interpolation procedure to be able to provide values of the QTF on a much finer grid than the initial one to comply with the requirement of sufficient frequency resolution to capture the dynamics of slow-drift motion. In this paper cubic spline interpolation has been used. In the particular case considered here, the resolution requirement led to the choice ∆ω = 0.0018 rad/s and L = 760 interpolated discrete frequencies.

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

227

Table 2 Representative sea states along with response standard deviations, nonlinear and linear TLP Hs (m)

T p (s)

σ Z (m), lin.

σ Z (m), nonlin.

10.0 15.0

11 17

9.3 6.4

8.2 6.0

Fig. 3. Wave exciting force amplitude, surge QTF.

Fig. 2. Wave exciting force amplitude, surge LTF.

Each random stationary sea state is specified by a JONSWAP spectrum, which is given as follows ( αg 2 5  ω p 4 Sη (ω) = 5 exp − + ln γ 4 ω ω "  2 #) 1 ω (25) −1 × exp − 2 ωp 2σ where g = 9.81 m s−2 , ω p denotes the peak frequency in rad/s and where γ and σ are parameters affecting the spectral shape. σ = 0.07 when ω ≤ ω p , and σ = 0.09 when ω > ω p . The parameter γ is chosen to be equal to 3.3. The parameter α is determined from the following empirical relationship !2 Hs α = 5.06 (26) (1 − 0.287 ln γ ) T p2 where Hs denotes the significant wave height and T p = 2π/ω p is the spectral peak wave period. For the present analysis, two short-term JONSWAP sea states have been used — moderate and severe. Table 2 presents the sea state parameters along with the corresponding response standard deviations. Fig. 2 shows the LTF for the wave exciting force amplitude, while Fig. 3 depicts the spline interpolated QTF. Since the QTF is complex valued, only its absolute value is plotted. For the two chosen sea states, Figs. 4–7 present the corresponding response tail crossing rates obtained by Monte Carlo simulation for the linear system given by putting ε = 0 and c˜ = 0 in Eq. (9). Figs. 4 and 6 show the results obtained from 1000 realizations, each of duration 1 h.

Fig. 4. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 1000 h of response time histories for the case of linear dynamics (c˜ = ε = 0). Saddle point integration results (—) coincide with the optimized linear fit with b = 0.75σ Z , q = 0.0205, where σ Z = 9.3 m (see Table 2). Sea state with Hs = 10 m, T p = 11 s.

The crossing rate plots are done on the transformed scale, i.e. log | log(ν + (ξ )/q)|versus log((ξ − b)/σ Z ), see Eq. (22). Extreme response prediction based on Eq. (17) will typically involve crossing rates of the orders 10−6 − 10−7 , but to illustrate the achieved accuracy the extrapolated results at the crossing rate level 10−10 are highlighted. Figs. 4 and 6 also show the highly accurate results obtained by using a saddle point integration technique [13]. These results cannot be distinguished from those obtained by linear extrapolation of the mean upcrossing rate function provided by Monte Carlo simulations. Hence, linear extrapolation can be done over several orders of magnitude with high accuracy. To illustrate the fact that good accuracy can be obtained with much shorter time simulation records, Figs. 5 and 7 show the results of using only 25 h of simulated response time histories, which required about 10 min. of CPU time on a standard desktop PC. Figs. 8–11 present response tail crossing rates for the two sea states obtained by Monte Carlo simulation for the fully

228

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

Fig. 5. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 25 h of response time histories, and by optimized linear fit (—) with b = 0.75σ Z , q = 0.020, where σ Z = 9.3 m (see Table 2), for the case of linear dynamics (c˜ = ε = 0). Sea state with Hs = 10 m, T p = 11 s.

Fig. 6. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 1000 h of response time histories for the case of linear dynamics (c˜ = ε = 0). Saddle point integration results (—) coincide with the optimized linear fit with b = 0.39σ Z , q = 0.042, where σ Z = 6.4 m (see Table 2). Sea state with Hs = 15 m, T p = 17 s.

nonlinear model given by Eq. (9). In this case the Monte Carlo simulation results are the only results available for verification of the extrapolation method for the nonlinear model, but the experience from the linear case indicates that the results obtained from Figs. 8 and 10 are very accurate. Again, the crossing rate plots are done on the transformed scale, see Eq. (22), and for illustration purposes the extrapolated results at the crossing rate level 10−10 are highlighted. These results can be compared with the corresponding results obtained from only 25 h of simulated response time histories shown in Figs. 9 and 11. The agreement is again good.

Fig. 7. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 25 h of response time histories for the case of linear dynamics (c˜ = ε = 0); optimized linear fit (—) with b = 0.39σ Z , q = 0.042, where σ Z = 6.4 m (see Table 2). Sea state with Hs = 15 m, T p = 17 s.

Fig. 8. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 1000 h of response time histories for the fully nonlinear model; optimized linear fit (—) with b = 0.54σ Z , q = 0.0248, where σ Z = 8.2 m (see Table 2). Sea state with Hs = 10 m, T p = 11 s.

6. Concluding remarks In this paper we have described a method for predicting the extreme response of floating marine structures in random waves based on Monte Carlo simulations. The prerequisite for application of the method is thus the availability of the equations of motion of the structure which makes possible the simulation of its response time histories caused by stochastic environmental loads. The method consists of a combination of the basic Monte Carlo method with an optimization based extrapolation technique. The method has been successfully implemented, and

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

Fig. 9. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 25 h of response time histories for the fully nonlinear model; optimized linear fit (—) with b = 0.54σ Z , q = 0.0234, where σ Z = 8.2 m (see Table 2). Sea state with Hs = 10 m, T p = 11 s.

229

Fig. 11. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 25 h of response time histories for the fully nonlinear model; optimized linear fit (—) with b = 0.39σ Z , q = 0.045, where σ Z = 6.0 m (see Table 2). Sea state with Hs = 15 m, T p = 17 s.

Technology. The authors would also like to acknowledge the support of the StatoilHydro Research Centre. References

Fig. 10. Empirical crossing rates by Monte Carlo simulation (∗) with 95% confidence bands (– – –) based on 1000 h of response time histories for the fully nonlinear model; optimized linear fit (—) with b = 0.39σ Z , q = 0.045, where σ Z = 6.0 m (see Table 2). Sea state with Hs = 15 m, T p = 17 s.

numerical results for the surge response of a TLP model have been presented. For two chosen sea states, complete agreement has been demonstrated between the predictions provided by the proposed method and an exact method for the case of linear dynamics. To verify the accuracy of the proposed method for the general case of nonlinear dynamics also, extensive Monte Carlo simulations were carried out for the same two sea states, and both cases lend further support to the applicability and accuracy of the method. Acknowledgements This work was supported by the Research Council of Norway (NFR) through the Centre for Ships and Ocean Structures (CeSOS) at the Norwegian University of Science and

[1] Neal E. Second order hydrodynamic forces due to stochastic excitation. In: Proc. 10th ONR symposium. 1974. [2] Vinje T. On the statistical distribution of second-order forces and motions. International Shipbuilding Progress 1983;30:58–68. [3] Naess A. Statistical analysis of second-order response of marine structures. Journal of Ship Research 1985;29(4):270–84. [4] Langley RS. A statistical analysis of low frequency second-order forces and motions. Applied Ocean Research 1987;9(3). [5] Kato S, Ando S, Kinoshita T. On the statistical theory of total second order response of moored floating structures. In: Proc. 19th annual offshore technology conference. Paper No. 5579. 1987. [6] Naess A. Statistical analysis of nonlinear, second-order forces and motions of offshore structures in short-crested random seas. Probabilistic Engineering Mechanics 1990;5(4):192–203. [7] Naess A, Johnsen JM. An efficient numerical method for calculating the statistical distribution of combined first-order and wave-drift response. Journal of Offshore Mechanics and Arctic Engineering 1992;114(3): 195–204. [8] Langley RS, McWilliam S. A statistical analysis of first and second order vessel motions induced by waves and wind gusts. Applied Ocean Research 1993;15(1):13–23. [9] McWilliam S, Langley RS. Extreme values of first and second order wave induced vessel motions. Applied Ocean Research 1993;15(3):169–81. [10] Grime AJ, Langley RS. On the efficiency of crossing rate prediction methods used to determine extreme motions of moored offshore structures. Applied Ocean Research 2003;25(3):127–35. [11] Naess A. Crossing rate statistics of quadratic transformations of gaussian processes. Probabilistic Engineering Mechanics 2001;16(3):209–17. [12] Naess A, Karlsen HC. Numerical calculation of the level crossing rate of second order stochastic Volterra systems. Probabilistic Engineering Mechanics 2004;19(2):155–60. [13] Naess A, Karlsen HC, Teigen PS. Numerical methods for calculating the crossing rate of high and extreme response levels of compliant offshore structures subjected to random waves. Applied Ocean Research 2006; 28(1):1–8. [14] Naess A, Johnsen JM. Response statistics of nonlinear, compliant offshore structures by the path integral solution method. Probabilistic Engineering Mechanics 1993;8(2):91–106.

230

A. Naess et al. / Applied Ocean Research 29 (2007) 221–230

[15] Adegeest L, Braathen A, Vada T. Evaluation of methods for estimation of extreme nonlinear ship responses based on numerical simulations and model tests. In: Proc. 22nd symposium on naval hydrodynamics (1998), vol. 1. Washington: National Academy Press; 1999. p. 70–84. [16] Naess A, Gaidai O, Teigen PS. Simulation based estimation of extreme response of floating offshore structures. In: Proceedings 26th international conference on offshore mechanics and arctic engineering, New York: ASME. 2007. p. OMAE–2007–29079. [17] Teigen PS. Nonlinear wave loads on the shafts of a GBS. In: Proceedings of ISOPE’2004. 2004 ISOPE Paper IL-01. [18] Rice SO. Mathematical analysis of random noise. In: Wax N, editor. Selected papers on noise and stochastic processes. New York: Dover Publications, Inc.; 1954. p. 133–294. [19] Naess A. On the long-term statistics of extremes. Applied Ocean Research 1984;6(4):227–8.

[20] Naess A, Moan T. Probabilistic design of offshore structures. In: Chakrabarti SK, editor. Handbook of offshore engineering, vol. 1. Amsterdam: Elsevier; 2005. p. 197–277 [chapter 5]. [21] Walpole RE, Myers RH, Myers SL, Ye K. Probability and statistics for engineers and scientists. Upper Saddle River (NJ): Prentice Hall; 2002. [22] Naess A, Gaidai O. Monte Carlo methods for estimating the extreme response of dynamical systems. Journal of Engineering Mechanics, ASCE. 2008 [in press]. [23] Leadbetter MR, Lindgren G, Rootzen H. Extremes and related properties of random sequences and processes. New York: Springer-Verlag; 1983. [24] Naess A, Gaidai O, Haver S. Efficient estimation of extreme response of drag dominated offshore structures by Monte Carlo simulation. Ocean Engineering 2007;34(16):2188–97. [25] Naess A, Gaidai O. Estimation of extreme values from sampled time series. Structural Safety 2008 [submitted for publication].