1. Force models
For sinusoidal motion, the equation of motion may be written
where:
\(\boldsymbol{M}\) 
frequencydependent mass matrix 
\(\boldsymbol{m}\) 
body mass matrix 
\(\boldsymbol{A}\) 
frequencydependent addedmass 
\(\boldsymbol{C}\) 
frequencydependent potential damping matrix 
\(\boldsymbol{D_1}\) 
linear damping matrix 
\(\boldsymbol{D_2}\) 
quadratic damping matrix 
\(\boldsymbol{f}\) 
vector function where each element is given by \(f_i = \dot{x}_i \dot{x}_ i\) 
\(\boldsymbol{k}\) 
hydrostatic stiffness matrix 
\(\boldsymbol{x}\) 
position vector 
\(\boldsymbol{q}\) 
exciting force vector 
The exciting forces on the righthand side of Equation (1) are given by
where:
stationkeeping and coupling elements, etc.)
\(\boldsymbol{q}_{WI}\) 
wind drag force 
\(\boldsymbol{q}_{WA}^{(1)}\) 
first order wave excitation force 
\(\boldsymbol{q}_{WA}^{(2)}\) 
second order wave excitation force 
\(\boldsymbol{q}_{CU}\) 
current drag force 
\(\boldsymbol{q}_{ext}\) 
any other forces (wave drift damping, specified forces and forces from stationkeeping and coupling elements, etc.) 
Two different solution methods described in the following two subsections are available in SIMO.
1.1. Solution by convolution integral
Assume that the equations of motion can be written
With regard to frequencydependent coefficients only, the equation of dynamic equilibrium can be written
provided that the right hand side force varies sinusoidally at one single frequency, \(\mathrm {\omega }\).
In the frequency domain, the equation is written
or
Using the following relations,
where \(\mathrm {\boldsymbol{A}_\infty=\boldsymbol{A}(\omega =\infty)}\) and \(\mathrm {\boldsymbol{C}_\infty=\boldsymbol{C}(\omega =\infty)\equiv0}\), Equation (7) is written
Using the inverse Fourier transform gives
Physically, values of \(\mathrm {\boldsymbol{h}(t\tau)}\) for, \(\mathrm {t<0}\), i.e. before the "experiment" started, are zero. Also, causality implies that \(\mathrm {\boldsymbol{h}(t\tau)\equiv0}\) for \(\mathrm {\tau>t}\).
Substituting \(\mathrm {\boldsymbol{f}(t)}\) from Equation (5) and \(\mathrm {\boldsymbol{f}'(t)}\) from Equation (4), the equation of motion becomes
\(\mathrm {\boldsymbol{h}(\tau)}\), the retardation function, is computed by a transform of the frequencydependent addedmass and damping
or similarly
Note that this transform pair is identical to the Fourier transform defined in General Assumptions and Notation.
Using the fact that \(\mathrm {\boldsymbol{c}(\omega )=\boldsymbol{c}(\omega )}\) and \(\mathrm {\boldsymbol{a}(\omega )=\boldsymbol{a}(\omega )}\) gives
From causality, \(\mathrm {\boldsymbol{h}(\tau)=0}\) for \(\mathrm {\tau<0}\); the process can not have any memory effect of the future. This means that the two parts in the integral, Equation (15) must be opposite for \(\mathrm {\tau<0}\) and identical for \(\mathrm {\tau>0}\), or mathematically:
for \(\mathrm {\tau>0}\).
This means that the frequencydependent mass and damping can be found from the retardation function
For \(\mathrm {\tau=0}\) we have from Equation (15)
or, if \(\mathrm {\boldsymbol{c}}\) is not known, we make use of the fact that \(\mathrm {\boldsymbol{c}(\omega =0)=0}\), which gives
The relations in Equation (17) and frequencydependent addedmass and damping are known as the KramersKrĂ¶nig relation.
Either frequencydependent addedmass or frequencydependent damping, and one value of the added mass, are required to calculate the retardation function. In the program, frequencydependent damping is used for calculating retardation functions.
1.1.1. Method used for retardation function calculation
The retardation function is calculated on the basis of frequency dependent demping according to Equation (16). A Fast Fourier Transform (FFT) that requires equidistant frequency resolution is applied for the inverse Fourier transform.
The specified (truncated) frequency depended damping is given for the frequency range \(\mathrm {[{\omega }_{min}^{T},{\omega }_{max}^{T}]}\). The values of damping at \(\mathrm {{\omega }_{min}^{T}}\) and \(\mathrm {{\omega }_{max}^{T}}\) are given by \(\mathrm {C({\omega }_{min}^{T})}\) and \(\mathrm {C({\omega }_{max}^{T})}\) .
The following steps are applied:

Add damping values on both sides of the truncated frequency range of damping to ensure that the values decreases to zero both in the low frequency and the high frequency range:

in the low frequency range by setting \(\mathrm {C(\omega =0)=0}\) and damping values between \(\mathrm {(\omega =0)}\) and \(\mathrm {C({\omega }_{min}^{T})}\) assuming that damping value increases with \(\mathrm {{\omega }^{2}}\)

in the high frequency range assuming that the damping decrease from \(\mathrm {C({\omega }_{max}^{T})}\) with \(\mathrm {{\omega }^{3}}\)


Resample damping values with \(\mathrm {\Delta \omega =\frac{2\pi }{n\cdot \Delta T}}\) where \(\mathrm {n}\) is number of time steps (initially) and \(\mathrm {\Delta T}\) is time step length of the retardation function

Compute the retardation function by inverse Fourier transfom of the resampled frequency dependent damping according to Equation (16)

Truncate the length of the retardation function. The default limit value is set to 0.5% of the maximum absolute retardation function value. The retardation function is truncated by searching backwards until this limit value is exceeded.

Correct the truncated retardation function so that the sum of all values in the retardation function is zero. All values are shifted with the same amount.
1.1.2. Method used to establish added mass at infinte frequency
The following steps are applied:

Fourier transform the truncated and corrected retardation functions accordng to Equation (17). This gives the shifted added mass at equidistant FFT frequencies, \(\mathrm {\boldsymbol{a}(\omega )}\), Equation (8).

Interpolate the shifted added mass to find values at the original truncated frequency range \(\mathrm {[{\omega }_{min}^{T},{\omega }_{max}^{T}]}\). These frequencies are the same for specified frequency dependent added mass and damping.

Accumulate the difference between the original and the interpolated added mass from the previous step for all original frequencies. The added mass at infinite frequency is found by dividing this value by the number of original frequencies.
1.2. Separation of motions
As an alternative to solving the whole differential equation Equation (1) in time domain by use of the retardation function, the motions can be separated in a high frequency part and a low frequency part.
The high frequency motions are solved in the frequency domain, which require the motions to be linear responses to waves. This means that \(\mathrm {\boldsymbol{D}_2}\), the quadratic damping, is set to be zero and \(\mathrm {\boldsymbol{K}}\) constant.
The exciting force is separated in a high frequency part, \(\mathrm {\boldsymbol{q}^{(1)}}\) and a low frequency part, \(\mathrm {\boldsymbol{q}^{(2)}}\).
The position vector can then be separated,
The high frequency motions to be solved in frequency domain are expressed by
In the frequency domain, using the equation for first order wave excitation force, we have
where \(\mathrm {\boldsymbol{H}_1}\) is the first order transfer function between excitation force and wave elevation and \(\mathrm {\boldsymbol{X}_{HF}}\) is the first order transfer function between motion and wave elevation.
The low frequency motions are solved in time domain for which the dynamic equilibrium equation is written