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{C}(\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 damping according to Equation (16) by using the Fast Fourier Transform (FFT). The FFT requires the frequency dependent damping to be sampled at equally spaced frequencies between 0 and the Nyquist frequency, which involves both interpolation and extrapolation:

  • Below the zero-frequency and the lowest specified frequency, the frequency dependent damping is assumed to grow with \(\mathrm {{\omega }^{2}}\) from zero at zero frequency to the specified value at the lowest specified frequency.

  • Above the highest specified frequency, the frequency dependent damping is assumed to decrease with \(\mathrm {{\omega }^{-3}}\) from the specified value at the highest specified frequency.

  • Linear interpolation is used between the specified frequencies.

Relevant background information can be found in (Ogilvie, 1964), Greenhow(1986), Taghipour et. al. (2008). Note that these references are general and does not describe how this is implemented in SIMO.

1.1.2. Method used to establish added mass at infinite frequency

Added mass at infinite frequency should preferably be computed directly by WADAM or other supported boundary element solver. However, if SIMA does not find the infinite frequency added mass when importing hydrodynamic data, SIMA will compute the added mass at infinite frequency from the retardation functions by combining eq. Equation (8) and Equation (17).

  1. Fourier transform the truncated and corrected retardation functions according 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}\]