1. Force models

For sinusoidal motion, the equation of motion may be written

\[\boldsymbol{M\ddot {x}}+\boldsymbol{C\dot {x}}+\boldsymbol{D_1\dot {x}}+\boldsymbol{D_2f}(\boldsymbol{\dot {x}})+\boldsymbol{K}(\boldsymbol{x})\boldsymbol{x}=\boldsymbol{q}(t,\boldsymbol{x},\boldsymbol{\dot {x}})\]
\[\begin{array}{lllll}\boldsymbol{M}&=&\boldsymbol{m}+\boldsymbol{A}(\omega )\\\\\boldsymbol{A}(\omega )&=&\boldsymbol{A}_\infty+\boldsymbol{a}(\omega )\\\\\boldsymbol{A}_\infty&=&\boldsymbol{A}(\omega =\infty)\\\\\boldsymbol{C}(\omega )&=&\boldsymbol{C}_\infty+\boldsymbol{c}(\omega )\\\\\boldsymbol{C}_\infty&=&\boldsymbol{C}(\omega =\infty)&\equiv&\boldsymbol{0}\end{array}\]

where:

\(\boldsymbol{M}\)

frequency-dependent mass matrix

\(\boldsymbol{m}\)

body mass matrix

\(\boldsymbol{A}\)

frequency-dependent added-mass

\(\boldsymbol{C}\)

frequency-dependent 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 right-hand side of Equation (1) are given by

\[\boldsymbol{q}(t,\boldsymbol{x},\boldsymbol{\dot {x}})=\boldsymbol{q}_{WI}+\boldsymbol{q}_{WA}^{(1)}+\boldsymbol{q}_{WA}^{(2)}+\boldsymbol{q}_{CU}+\boldsymbol{q}_{ext}\]

where:

station-keeping 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 station-keeping 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

\[(\boldsymbol{m}+\boldsymbol{A}(\omega ))\boldsymbol{\ddot {x}}+\boldsymbol{C}(\omega )\boldsymbol{\dot {x}}+\boldsymbol{Kx}=\boldsymbol{f}'(t)=\boldsymbol{q}-\boldsymbol{D}_2\boldsymbol{f}(\boldsymbol{\dot {x}})-\boldsymbol{D}_1\boldsymbol{\dot {x}}\]

With regard to frequency-dependent coefficients only, the equation of dynamic equilibrium can be written

\[\boldsymbol{A}(\omega )\boldsymbol{\ddot {x}}+\boldsymbol{C}(\omega )\boldsymbol{\dot {x}}=\boldsymbol{f}(t)=\boldsymbol{f}'(t)-\boldsymbol{Kx}-\boldsymbol{m\ddot {x}}\]

provided that the right hand side force varies sinusoidally at one single frequency, \(\mathrm {\omega }\).

In the frequency domain, the equation is written

\[\big(-\omega ^2\boldsymbol{A}(\omega )+i\omega \boldsymbol{C}(\omega )\big)\boldsymbol{X}(\omega )=\boldsymbol{F}(\omega )\]

or

\[\big(i\omega \boldsymbol{A}(\omega )+\boldsymbol{C}(\omega )\big)i\omega \boldsymbol{X}(\omega )=\boldsymbol{F}(\omega )\]

Using the following relations,

\[\begin{array}{l}\boldsymbol{A}(\omega )=\boldsymbol{A}_\infty+\boldsymbol{a}(\omega )\\\\\boldsymbol{B}(\omega )=\boldsymbol{C}_\infty+\boldsymbol{c}(\omega )\end{array}\]

where \(\mathrm {\boldsymbol{A}_\infty=\boldsymbol{A}(\omega =\infty)}\) and \(\mathrm {\boldsymbol{C}_\infty=\boldsymbol{C}(\omega =\infty)\equiv0}\), Equation (7) is written

\[-\omega ^2\boldsymbol{A}_\infty\boldsymbol{X}(\omega )+\big(i\omega \boldsymbol{a}(\omega )+\boldsymbol{c}(\omega )\big)i\omega \boldsymbol{X}(\omega )=\boldsymbol{F}(\omega )\]

Using the inverse Fourier transform gives

\[\boldsymbol{A}_\infty\boldsymbol{\ddot {x}}(t)+\int^\infty_{-\infty}\boldsymbol{h}(t-\tau)\boldsymbol{\dot {x}}(\tau)\mathrm {d}\tau=\boldsymbol{f}(t)\]

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}\).

\[\boldsymbol{A}_\infty\boldsymbol{\ddot {x}}(t)+\int_0^t\boldsymbol{h}(t-\tau)\dot {x}(\tau)\mathrm {d}\tau=\boldsymbol{f}(t)\]

Substituting \(\mathrm {\boldsymbol{f}(t)}\) from Equation (5) and \(\mathrm {\boldsymbol{f}'(t)}\) from Equation (4), the equation of motion becomes

\[(\boldsymbol{m}+\boldsymbol{A}_\infty)\ddot {x}+\boldsymbol{D}_1\boldsymbol{\dot {x}}+\boldsymbol{D}_2\boldsymbol{f}(\boldsymbol{\dot {x}})+\boldsymbol{Kx}+\int_0^t\boldsymbol{h}(t-\tau)\boldsymbol{\dot {x}}(\tau)\mathrm {d}\tau=\boldsymbol{q}(t,\boldsymbol{x},\boldsymbol{\dot {x}})\]

\(\mathrm {\boldsymbol{h}(\tau)}\), the retardation function, is computed by a transform of the frequency-dependent added-mass and damping

\[\boldsymbol{h}(\tau)=\frac{1}{2\pi }\int^\infty_{-\infty}\boldsymbol{c}(\omega )+i\omega \boldsymbol{a}(\omega )e^{i\omega t}\mathrm {d}\omega =\frac{1}{2\pi }\int^\infty_{-\infty}\boldsymbol{H}(\omega )e^{i\omega t}\mathrm {d}\omega\]

or similarly

\[\boldsymbol{H}(\omega )=\int^\infty_{-\infty}\boldsymbol{h}(\tau)e^{-i\omega \tau}\mathrm {d}\tau=\boldsymbol{c}(\omega )+i\omega \boldsymbol{a}(\omega )\]

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

\[\boldsymbol{h}(\tau)=\frac{1}{\pi }\int_0^\infty\big(\boldsymbol{c}(\omega )\cos(\omega \tau)-\omega \boldsymbol{a}(\omega )\sin(\omega \tau)\big)\mathrm {d}\omega\]

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:

\[\boldsymbol{h}(\tau)=\frac{2}{\pi }\int^\infty_0\boldsymbol{c}(\omega )\cos(\omega \tau)\mathrm {d}\omega =-\frac{2}{\pi }\int^\infty_0\omega \boldsymbol{a}(\omega )\sin(\omega \tau)\mathrm {d}\omega\]

for \(\mathrm {\tau>0}\).

This means that the frequency-dependent mass and damping can be found from the retardation function

\[\boldsymbol{a}(\omega )=-\frac{1}{\omega }\int^\infty_0\boldsymbol{h}(\tau)\sin(\omega \tau)\mathrm {d}\tau\]
\[\boldsymbol{c}(\omega )=-\int^\infty_0\boldsymbol{h}(\tau)\cos(\omega \tau)\mathrm {d}\tau\]

For \(\mathrm {\tau=0}\) we have from Equation (15)

\[\boldsymbol{h}(0)=\frac{1}{\pi }\int^\infty_0\boldsymbol{c}(\omega )\mathrm {d}\omega\]

or, if \(\mathrm {\boldsymbol{c}}\) is not known, we make use of the fact that \(\mathrm {\boldsymbol{c}(\omega =0)=0}\), which gives

\[\int^\infty_0\boldsymbol{h}(\tau)\mathrm {d}\tau=0\]

The relations in Equation (17) and frequency-dependent added-mass and damping are known as the Kramers-Krönig relation.

Either frequency-dependent added-mass or frequency-dependent damping, and one value of the added mass, are required to calculate the retardation function. In the program, frequency-dependent 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:

  1. 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}}\)

  2. 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

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

  4. 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.

  5. 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:

  1. 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).

  2. 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.

  3. 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)}}\).

\[\begin{array}{lll}\boldsymbol{q}(t,\boldsymbol{x},\boldsymbol{\dot {x}})&=&\boldsymbol{q}^{(1)}+\boldsymbol{q}^{(2)}\\\\\boldsymbol{q}^{(1)}&=&\boldsymbol{q}_{WA}^{(1)}\\\\\boldsymbol{q}^{(2)}&=&\boldsymbol{q}_{WI}+\boldsymbol{q}_{WA}^{(2)}+\boldsymbol{q}_{CU}+\boldsymbol{q}_{ext}\end{array}\]

The position vector can then be separated,

\[\boldsymbol{x}=\boldsymbol{x}_{LF}+\boldsymbol{x}_{HF}\]

The high frequency motions to be solved in frequency domain are expressed by

\[\boldsymbol{m}+\boldsymbol{A}(\omega )\boldsymbol{\ddot {x}}_{HF}+\boldsymbol{D}_1+\boldsymbol{C}(\omega )\boldsymbol{\dot {x}}_{HF}+\boldsymbol{Kx}_{HF}=\boldsymbol{q}_{WA}^{(1)}(\omega )\]

In the frequency domain, using the equation for first order wave excitation force, we have

\[\boldsymbol{X}_{HF}(\omega )=\Big(-\omega ^2\big(\boldsymbol{m}+\boldsymbol{A}(\omega )\big)+i\omega \boldsymbol{D}_1+\boldsymbol{C}(\omega )+\boldsymbol{K}\Big)^{-1}\boldsymbol{H}_1(\omega )\tilde{\zeta}(\omega )\]

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

\[\boldsymbol{m}+\boldsymbol{A}(\omega =0)\boldsymbol{\ddot {x}}_{LF}+\boldsymbol{D}_1\boldsymbol{\dot {x}}_{LF}+\boldsymbol{D}_2\boldsymbol{f}(\boldsymbol{\dot {x}})+\boldsymbol{Kx}_{LF}=\boldsymbol{q}^{(2)}=\boldsymbol{q}_{WI}+\boldsymbol{q}_{WA}^{(2)}+\boldsymbol{q}_{CU}+\boldsymbol{q}_{ext}\]