Generation of time series

This chapter describes different ways of generating realizations (samples) of Gaussian stochastic processes for which the mean values and the variance spectra are known. These are as follows:

  • wind gust.

  • wave elevation, particle velocity and acceleration.

  • \(1^{\mathrm {st}}\) order wave responses such as \(1^{\mathrm {st}}\) order excitations, \(1^{\mathrm {st}}\) order wave-induced motions or \(1^{\mathrm {st}}\) order diffracted wave kinematics.

  • \(2^{\mathrm {nd}}\) order wave responses such as full or simplified second order wave forces.

The time series are generated by superposition of harmonic components, with uniformly distributed phases by means of pre-generation by the Fast Fourier Transform (FFT), by time domain summation of the harmonic components or by a state-space model (wind gust) driven by white noise.

1. Superposition of harmonic components

The contribution to the variance of each harmonic component is set equal to the spectrums contribution in the area

The time series are generated by discretizing the variance spectrum into a finite number of finite-valued harmonic components, and by sampling phases from a uniform distribution over \(\mathrm {[0,2\pi ]}\) , see Equation (1). The time series of the sample is written:

\[x(t)=\sum_{k=1}^{N}\,x_k\cos(\omega t+\phi _k)\]

where \(\mathrm {\phi _k}\) is the sampled phase. Processes that are linear transformations of \(\mathrm {x(t)}\) are written:

\[y(t)=\sum_{k=1}^{N}\,\frac{y_k}{x_k}x_k\cos(\omega t+\phi _k+\phi _y)\]

where:

\(\displaystyle \frac{y_{k}}{x_k}\)

response amplitude operator

\(\phi_y\)

phase angle

The definitions of phase angles are described in Environment.

Time series generated in this way will repeat themselves with a period \(T_p=2\pi /\Delta \omega _{k\mathrm {min}}\) , where \(\omega _{k\mathrm {min}}\) is the smallest frequency increment.

The user can currently choose between two different algorithms for generating pseudo-random numbers. The "mersenne twister" is the recommended method and should be used unless backwards compatibility with earlier SIMO versions are required. The "legacy" method is included for backwards compatibility. It has been identified that the legacy method can give non-gaussian and non-stationary wave elevation in short crested waves with more than about 30-50 discrete directions, depending on wave spectrum and simulation duration. By choosing the mersenne twister, these issues are avoided.

1.1. Fast Fourier transform (FFT)

Normally, the addition of harmonic components to obtain time series is done by the Fast Fourier transform (FFT) which requires equal frequency spacing and a number of frequencies (or number of time samples) to be \(\mathrm {N=2^r}\) where \(\mathrm {r}\) is an integer number.

A Cooley-Tukey Fourier transform algorithm is used to compute FFT in SIMO.

The relations between time increment, \(\mathrm {\Delta t}\) , the number of time steps, \(\mathrm {N_t}\) , number of positive frequency components, \(\mathrm {N_{\omega }^+}\) , and the frequency increment is:

\[\Delta \omega =\frac{2\pi }{N_t\Delta t}\]
\[N_{\omega }^+=\frac{N_t}{2}+1\]

Thus, the duration of the time series is limited to:

\[T_p=N_t\Delta t=\frac{2\pi }{\Delta \omega }\]

When using the FFT algorithm, the time series of excitations or linear responses must be calculated for predefined positions and vessel headings. It is not possible to utilise the very high efficiency of the FFT if positions must be modified by for example calculating for short periods and then updating the position. In this case, different strategies are possible:

  • Calculation with traditional addition of time-domain sine or cosine functions for every harmonic component at the wanted position and time increment. This is described in Section Summation of harmonic components in the time domain.

  • Introduction of a linear frequency-dependent phase filter (not implemented).

  • Calculation of time series by FFT in a grid of positions and directions and interpolation during the simulation to obtain values at the wanted position and time increment (not implemented).

1.2. Summation of harmonic components in the time domain

Calculation of wave responses for each time step during the simulation can be performed by summation of harmonic components (sine/cosine series). The calculation is made for the instantaneous locations of each body (and for slender elements: each section). This will normally give more correct wave kinematics and loads for each body, irrespective of changes in the (x,y,z)-position. This method will enable simulation of horizontal and vertical transport and also, by some stretching, enable changes in the sea state with time. The computing time will be longer than for an FFT pre-generation, especially for long time series.

Due to these time-consuming calculations, a combination of pre-generated time series and cosine series in the time domain is made possible. For such a combination it is also possible to ensure that the same realization of wave components is used for both options for wave response calculation.

Time series from simulations with system changes including position changes may be in conflict with the assumption of stationary (ergodic) Gaussian processes. Hence, the analysis of these time series and their resulting responses should be handled with care.

As stated above, the time series generated in this way will repeat themselves with a period \(T_p=2\pi /\Delta \omega _{k\mathrm {min}}\) . For a user-given number of wave components for wind sea and swell, it is important that the number of components is sufficiently high to avoid a repeating wave response time series for the longest duration of the simulated operation that is close to stationary.

The usual methods of estimating extreme response are not applicable for processes with a moving average or which includes transients. Thus the standard deviation and estimates of extremes by use of the Rayleigh distribution should not be calculated for a body that i.e. is lowered through the splash zone or is lifted or landed during the simulation. Problems with repeating time series of wave kinematics, as described above, will thus not be a relevant subject in many of these cases.

Stochastic amplitudes may be specified to be used for time domain wave response generation. This is another conflict with the assumption of ergodic processes, and analyses will often require simulations with multiple wave realizations in order to be able to give reliable extreme value estimates.

1.3. State-space model driven by white noise

Fluctuating wind (gust) can be modelled as a state space model driven by white noise of unit density, (Kaasen, 1999). The first step is to approximate the given wind spectrum in non-dimensional form by a rational transfer function, so that

\[S'(f')=|H'(jf')|^2,\quad j=\sqrt{-1}\]

where the non-dimensional transfer function is on the form:

\[H'(p')=\frac{b'_0+b'_1p'+\cdots+b'_mp'^m}{a'_0+a'_1p'+\cdots+a'_np'^n}\]

where \(\mathrm {p'=jf'}\)

The two scaling factors, one for frequency \(\mathrm {(\lambda_f)}\) , and one for spectral density \(\mathrm {(\lambda_s)}\) are:

\[f=\lambda_ff', S=\lambda_sS'\]

The next step is to estimate the model parameters by iterative least squares fitting of the rational spectrum functions to the given spectrum and deriving transfer functions and hence state space models from these.

For frequencies approaching infinity, the spectral density functions follow \(\mathrm {const}\cdot f^{-5/3}\). By choosing \(\mathrm {n=m+1}\) in Equation (7) , the fitted spectra will decline as \(\mathrm {const}\cdot f^{-2}\) towards high frequencies, i.e. slightly higher than the given spectra.

The dimensionless coefficients were found by least squares fitting of each spectrum for 100 dimensionless frequencies over an exponential frequency range from \(\mathrm {f=0.1}\) to \(\mathrm {f=10}\) .

For the Davenport and Harris spectra it was found that satisfactory transfer functions could be obtained with \(2^{\mathrm {nd}}\) -order models. For the ISO 19901-1 (NPD) spectrum a \(3^{\mathrm {rd}}\) -order model was required. The coefficients for these spectra are given in Wind. To obtain coefficients for a given set of spectrum parameters, the following scaling factors for the model coefficients can be derived from the scaling factors \(\lambda_f\) and \(\lambda_s\):

\[\begin{array}{lll}a_i&=(2\pi \lambda_f)^{n-i}a'_i,&i=0,\cdots,n\\\\b_i&=(2\pi \lambda_f)^{m-i}\sqrt{\displaystyle \frac{\lambda_s}{2\pi }}b'_i,&i=0,\cdots,m\end{array}\]

The state space model is written:

\[\boldsymbol{\dot x}=\boldsymbol{Ax}+\boldsymbol{B}\zeta\]

Here, \(\mathrm {\boldsymbol{x}}\) is the \(\mathrm {n}\) -element state vector, and \(\mathrm {\boldsymbol{A}(n\times n)}\) and \(\mathrm {\boldsymbol{B}(n\times 1)}\) are constant matrices. \(\mathrm {{\zeta}}\) is white noise. For the \(2^{\mathrm {nd}}\) -order models the following state space representation may be chosen:

\[\begin{bmatrix}\dot x_1\\\dot x_2\end{bmatrix}=\begin{bmatrix}0&1\\-a_0&-a_1\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}+\begin{bmatrix}b_1\\b_0-a_1b_1\end{bmatrix}{\zeta}\]

where \(\mathrm {x_1}\) is the wind gust velocity and \(\mathrm {x_2}\) is an auxiliary variable.

For the \(3^{\mathrm {rd}}\) -order models a similar state space representation is chosen:

\[\begin{bmatrix}\dot x_1\\\dot x_2\\\dot x_3\\\end{bmatrix}=\begin{bmatrix}0&1&0\\0&0&1\\-a_0&-a_1&-a_2\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}+\begin{bmatrix}b_2\\b_1-a_2b_2\\b_0-a_2b_1-(a_1-a_2^2)b_2\end{bmatrix}\zeta\]

With excitation by a white-noise source of unit power density, these matrices now define a continuous-time state-space model that generates an output that has a wind gust spectrum equal to that of the specified spectrum.

2. Wind sampling

The wind gusts in mean wind direction and normal to mean wind direction are generated by either FFT or by a state-space model driven by white noise.

The wind gust in different positions can be fully correlated by using the same series or uncorrelated by generating stochastic independent time series at each point.

3. Sampling of waves and first-order responses

The harmonic components differ slightly from what was said in previous sections due to the directionality of waves.

Based on the assumption that the directionality of the wave elevation spectra can be written:

\[S_\zeta(\beta ,\omega )=\theta (\beta )S_\zeta(\omega )\]

then the surface elevation can be expressed as:

\[\zeta(t)=\sum_{j=1}^{N_\beta }\,\sum_{k=1}^{N_\omega }\,\mathrm {Re}\,(\tilde\zeta_{jk}e^{i\omega _kt})\]
\[|\zeta_{jk}|=\sqrt{2S_\eta (\omega _k)\Delta \omega _k\theta _j(\beta _j)\Delta \beta _j}\]
\[\arg\tilde\zeta_{jk}=\phi _{P_{jk}}+\phi _{\zeta_{jk}}-\frac{\pi }{2}\]

where:

\(\phi_{P_{jk}}\) position dependent phase angle

\(\phi_{\zeta_{jk}}\)

random phase angle

The subtraction of \(\mathrm {\pi /2}\) gives the wanted shape of the surface according to Waves.

\[\zeta(t)=\zeta_a\sin{(\omega t+\phi _\zeta)}\]

The harmonic components of velocity and acceleration of water particles and first order responses are found simply by multiplying their complex transfer functions with the corresponding harmonic component of the wave elevation.

The order of summation in Equation (7) is not important. Thus, it is favourable first to sum over the directions and then perform one FFT or sum over the frequencies.

4. Second-order time series

This section is limited to describing generation of time series of second-order forces. The quadratic term of the wave force is written (compare with Chapter Force Models):

\[q_{W\!A2}(t)=\int_{-\infty}^{\infty}\!{\int_{-\infty}^{\infty}\!{h^{(2)}(\tau_1,\tau_2)\zeta(t-\tau_2)}\textrm{d}{\tau_1}}\textrm{d}{\tau_2}\]

The second-order force can thus be generated by use of a:

  • quadratic filter

  • double sum of harmonic components

4.1. Quadratic filter

To be completed.

4.2. Double sum of harmonic components

The double sum of harmonic components can be written:

\[q_{W\!A}(t)=\mathrm {Re}\sum_{j}\,\sum_{k}\,\sqrt{S_{\eta }(\omega _j)S_{\eta }(\omega _k)\Delta \omega _j\Delta \omega _k}H_2(\omega _j,\omega _k)e^{i(\epsilon_j+\epsilon_k)}e^{i(\omega _j+\omega _k)t}\]

For pre-generation of waves this is obtained by a \(2^{\mathrm {nd}}\) order FFT.

4.3. Simplified model (Newman)

According to Newman, Chapter Force Models the wave force can be written:

\[q(t)=\overline{(L^+(t))^2}-\overline{(L^-(t))^2}\]

where

\[L^{\pm}=\mathrm {Re}\sum_{m}\,\zeta_m\sqrt{\pm2H_{mm}^{(2-)}e^{i\omega _mt}}\]

\(\mathrm {L^{\pm}}\) is found either by ordinary summation of the harmonic components or by use of FFT.

Alternatively,

\[u(t)=\sum_{m}\,\sqrt{H_{mm}^{(2-)}}\mathrm {Re}{(\tilde\zeta_me^{i\omega _mt})}\qquad v(t)=\sum_{m}\,\sqrt{H_{mm}^{(2-)}}\mathrm {Im}{(\tilde\zeta_me^{i\omega _mt})}\]

where \(\mathrm {H_{mm}^{(2-)}}\) can be positive and negative, and \(\mathrm {v(t)}\) is the Hilbert transform of \(\mathrm {u(t)}\) .

Then a result for the slowly varying force can then be seen from

\[q(t)=\mathrm {Re}(\sum_{m}\,\sum_{n}\,\tilde\zeta_m\tilde\zeta_n^*H_{mm}^{(2-)}e^{i(\omega _m-\omega _n)t})=\mathrm {Re}(u(t)^2+v(t)^2)\]

4.4. Short-crested waves

The effect of short-crested waves can be found by determining the contribution to the force from unidirectional waves and superimposing contributions from wave components to find the total force.