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 waveinduced 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 pregeneration by the Fast Fourier Transform (FFT), by time domain summation of the harmonic components or by a statespace 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 finitevalued 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:
where \(\mathrm {\phi _k}\) is the sampled phase. Processes that are linear transformations of \(\mathrm {x(t)}\) are written:
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 pseudorandom 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 nongaussian and nonstationary wave elevation in short crested waves with more than about 3050 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 CooleyTukey 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:
Thus, the duration of the time series is limited to:
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 timedomain 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 frequencydependent 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 pregeneration, especially for long time series.
Due to these timeconsuming calculations, a combination of pregenerated 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 usergiven 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. Statespace 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 nondimensional form by a rational transfer function, so that
where the nondimensional transfer function is on the form:
where \(\mathrm {p'=jf'}\)
The two scaling factors, one for frequency \(\mathrm {(\lambda_f)}\) , and one for spectral density \(\mathrm {(\lambda_s)}\) are:
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 199011 (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\):
The state space model is written:
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:
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:
With excitation by a whitenoise source of unit power density, these matrices now define a continuoustime statespace 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 statespace 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 firstorder 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:
then the surface elevation can be expressed as:
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.
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. Secondorder time series
This section is limited to describing generation of time series of secondorder forces. The quadratic term of the wave force is written (compare with Chapter Force Models):
The secondorder force can thus be generated by use of a:

quadratic filter

double sum of harmonic components
4.2. Double sum of harmonic components
The double sum of harmonic components can be written:
For pregeneration 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:
where
\(\mathrm {L^{\pm}}\) is found either by ordinary summation of the harmonic components or by use of FFT.
Alternatively,
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