Force models

1. General

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 Equation (34), 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}\]

2. Body forces

2.1. Mass forces

The mass forces are described by the product \(\mathrm {\boldsymbol{m\boldsymbol{\ddot {x}}}}\) where \(\mathrm {\boldsymbol{m}}\) is the (constant) mass matrix.

2.2. Time-dependent mass

The data group Time-Dependent Mass (TDM) allows the user to model the effects of a mass that can vary both in magnitude and in position.

A TDM can be used to model:

  1. Time variation of a mass located at a defined point at the body. In that case, the mass has no volume.

  2. Filling liquid into one or a set of tanks with specified geometry. Free surface effects can then be accounted for. In this case, one TDM may be modelled: - in terms of PORTIONS. Up to 5 PORTIONS are allowed. Each PORTION has a simple geometry (CONE or BOX). The TDM’s location and orientation relative to the body (vessel) needs to be specified. - in terms of a set of tanks. Each tank are represented by a mesh of triangular panels described in a STL (STereoLithography) file. The maximum number of ballast tanks is unconstrained.

  3. A mass that moves relative to the body (future extension)

In all cases, the effect of the force terms \(\frac{\mathrm{d}m}{\mathrm{d}t}\cdot v\) and \(\frac{\mathrm{d}I}{\mathrm{d}t}\cdot \omega\) are neglected (mass \(\mathrm {m}\) and inertia \(\mathrm {I}\)). This means that forces that are due to relative velocity between the TDM and the body are assumed to be negligible. In each time instant in a simulation, the mass matrix of the body will be updated and the gravity force due to the TDM will be calculated. In cases with tank filling, the centre of gravity of the liquid will also be updated. No sloshing forces are accounted for in the model and the surface of the liquid is always horizontal.

2.2.1. Operational relevance

A model of time-varying mass is useful for simulation of operations involving ballasting or de-ballasting, such as heavy lift operations in general and in particular mating and platform removal operations. Examples are a ballast tank that is emptied, or water that is moved from one tank to another, within a defined time period. An example of time dependent mass is presented in Example of time-dependent mass. If a TDM without geometry is used, the TDM mass may be negative.

2.2.2. Point mass model

Model characteristics and limitations

The model enables specification of one or more point masses on a body. Each point mass may vary with time. The location of the point is fixed in the local body coordinate system. The mass point can be a mass point or a tank with specified geometry. In both cases the mass is taken as an addition to the system mass specification.

Note that the mass will always contribute to gravity force.

Model description

The instantaneous point mass, \(\mathrm {m(t)}\), is found by using the prescribed mass rate. During static analysis the specified value at time zero \(\mathrm {m(t=0)}\), is used.

A point mass, \(\mathrm {m(t)}\), located at coordinates \(\mathrm {(x,y,z)}\) in local (body fixed) system, will give the following change in the body mass matrix.

\[\Delta \boldsymbol{M}(t)=\begin{bmatrix}1&0&0&0&z&-y\\\\0&1&0&-z&0&x\\\\0&0&1&y&-x&0\\\\0&-z&y&y^2+z^2&-xy&-zx\\\\z&0&-x&-xy&z^2+x^2&-yz\\\\-y&x&0&-zx&-yz&x^2+y^2\end{bmatrix}\cdot m(t)=\begin{bmatrix}\boldsymbol{I}&\boldsymbol{A}^T\\\\\boldsymbol{A}&\boldsymbol{A}\boldsymbol{A}^T\end{bmatrix}\cdot m(t)\]

where \(\mathrm {\boldsymbol{I}}\) is a 3x3 unit matrix and \(\mathrm {\displaystyle \boldsymbol{A}=\begin{bmatrix}0&-z&y\\\\z&0&-x\\\\-y&x&0\end{bmatrix}}\)

The updated mass matrix will then be \(\mathrm {\boldsymbol{M}(t=t_i)=\boldsymbol{M}_0+\Delta \boldsymbol{M}(t=t_i)}\)

where \(\mathrm {\boldsymbol{M}_0}\) is the mass matrix originally specified in the system description file.

A general expression of the inertia term in the equation of motion is

\[\frac{\partial (\boldsymbol{M}\boldsymbol{\dot {x}})}{\partial t}=\boldsymbol{M}\boldsymbol{\ddot {x}}+\boldsymbol{\dot {x}}\frac{\partial \boldsymbol{M}}{\partial t}\]

In the present implementation, the last term is not included, because the change in mass does not necessarily give an impulse to the body.

2.2.3. Tank with simple geometry

When a TDM with geometry is specified and the volume of one portion (portion B) shall be subtracted from portion A (IADD= -1), then no part of portion B may be outside portion A. There is no check in SIMO that this event does not occur. If the geometry of a TDM is defined by adding portions, the portions must not be overlapping (the intersection must be zero).

For a TDM with geometry, trim of the TDM in the range \(\mathrm {(89\deg<}\) trim \(\mathrm {<91\deg)}\) is not allowed. If this becomes a problem, the orientation of the tank x-axis relative to the vessel should be modified.

image243
Figure 1. Example of time-dependent mass.

2.2.4. Set of tanks with general geometry (ballast tank system)

When a TDM is defined as a ballast tank system, each tank geometry is represented by a triangular mesh described in a STL file. The number of tanks is unconstrained.

Computing the mass properties of a tank

A tank containing an amount of fluid contributes to the gravity force and to the mass matrix of the body. In order to compute this contribution, we assume that we know where the free surface of the fluid is relative to the tank. SIMO first proceeds to sort each panel of the mesh as follows:

  • if all the three vertices of the panel are above the free surface line, the panel is disregarded

  • if one vertex is above the free surface line, the panel is cut at the free surface line level. The intersection of the free surface line with the edges of the panel defines two vertices which lay on the free surface line. Two new triangular panels are then defined based on the two vertices under the free surface line and the two vertices on the free surface line. They are included in the "contributing panels".

  • if two vertices are above the free surface line, the panel is cut at the free surface line level. The intersection of the free surface line with the edges of the panel defines two vertices which lay on the free surface line. One new triangular panel is then defined based on the vertex under the free surface line and the two vertices on the free surface line. It is included in the "contributing panels".

  • if all the three vertices are under the free surface line, the panel is included in the "contributing panels".

Thanks to this cutting procedure, the mesh can be coarse as long as it describes the geometry correctly.

For each "contributing" panel, the following integrals are computed analytically:

  • \(\mathrm {V_{m0}=\int_{\Gamma}d\Gamma}\)

  • \(\mathrm {V_{m1}=\int_{\Gamma}\boldsymbol{X}d\Gamma}\)

  • \(\mathrm {V_{m2}=\int_{\Gamma}\boldsymbol{X}\cdot \boldsymbol{X^T}d\Gamma}\)

where:

  • \(\mathrm {\Gamma}\) is the integration domain defined by the volume "under" the triangular panel

  • \(\mathrm {\boldsymbol{X}}\) is the position vector in the Earth-fixed coordinate system

  • \(\mathrm {V_{m0}}\) is the contribution of the panel to the total volume

  • \(\mathrm {V_{m1}}\) is the contribution of the panel to the total volume center

  • \(\mathrm {V_{m2}}\) is the contribution of the panel to the total second volume moment

By summing up the contribution of each "contributing panel", the total fluid volume is computed, as well as the position of the center of volume. The two latter quantities are used to compute the total gravity force and moment, using the density of the fluid contained in the tank.

The second volume moments are used to build the mass matrix of the tank, which is transposed to the body-fixed coordinate system, using the generalized Huygens theorem.

Computing the position of the fluid free surface

One challenge consists in finding the position of the fluid free surface in the tank, given the attitude of the tank and a fluid volume. This is done in SIMO by using an iterative process based on an improved version of the secant method.

Permeability factor

The permeability factor is used to reduce the maximum amount of fluid a tank can contain. The total volume of a tank is computed by SIMO based on the geometry mesh of the tank. The maximum volume of fluid that can be contained in the tank is then equal to the total volume of the tank times the permeability factor.

Damaged tank

When the state of tank is "damaged", it is assumed that the tank is severely damaged. The fluid density in the tank is replaced by the sea water density. The amount of fluid in the tank is such that the fluid level corresponds to the calm sea water surface level. For example, if the tank is completely submerged, it will be completely filled with sea water.

Limitations

In the present SIMO version, it is not possible to change the amount of fluid in the a tank during a time domain simulation. The amount of fluid present in the tank is defined at the beginning of the simulation.

2.3. Hydrodynamic reaction forces

Based on first order potential wave theory, the hydrodynamic reaction forces are described by added-mass coefficients, \(\mathrm {\boldsymbol{A}}\), and damping coefficients, \(\mathrm {\boldsymbol{C}}\). Both added-mass and damping are frequency-dependent. They are included in the retardation functions, as described in Solution by convolution integral or the first-order motion transfer functions as described in Separation of motions.

2.4. Damping forces

Linear damping is expressed by \(\mathrm {\boldsymbol{D}_1\boldsymbol{\dot {x}}}\), where \(\mathrm {\boldsymbol{D}_1}\) is the linear damping matrix.

Quadratic damping is expressed by \(\mathrm {\boldsymbol{D}_2\boldsymbol{f}(\boldsymbol{\dot {x}})}\) where \(\mathrm {\boldsymbol{D}_2}\) is the quadratic damping matrix and \(\mathrm {\boldsymbol{f}}\) is a vector function where each element is given by \(\mathrm {f_i=\dot {x}_i|\dot {x}_i|}\).

2.5. Hydrostatic stiffness

2.5.1. Linear stiffness matrix model

The force is expressed by \(\mathrm {\boldsymbol{K(x-x_0)}}\) where \(\mathrm {\boldsymbol{K}}\) is the hydrostatic stiffness matrix and \(\mathrm {\boldsymbol{x_0}}\) the stiffness reference position.

2.5.2. Non-linear stiffness model

When using this formulation, the hydrostatic force model will return the complete buoyancy force (and not the restoring forces as it is the case when using linear stiffness matrix). The buoyancy force is computed based on the part of the hull geometry volume which is under the calm water line (waves are ignored).

The geometry of the hull is described by a STL (STereoLithography) file. The mesh of the geometry is then based on triangular panels. At each time step, SIMO will update the position of the body and proceed to a cut of the mesh as follows:

  • if all the three vertices of the panel are above the calm water line, the panel is disregarded

  • if one vertex is above the calm water line, the panel is cut at the calm water line level. The intersection of the calm water line with the edges of the panel defines two vertices which lay on the calm water line. Two new triangular panels are then defined based on the two vertices under the calm water line and the two vertices on the calm water line. They are included in the "contributing panels".

  • if two vertices are above the calm water line, the panel is cut at the calm water line level. The intersection of the calm water line with the edges of the panel defines two vertices which lay on the calm water line. One new triangular panel is then defined based on the vertex under the calm water line and the two vertices on the calm water line. It is included in the "contributing panels".

  • if all the three vertices are under the calm water line, the panel is included in the "contributing panels".

Thanks to this cutting procedure, it is not necessary to refine the mesh about the calm water line. The mesh can be coarse as long as it describes the geometry correctly.

For each "contributing" panel, the following integrals are computed analytically:

  • \(\mathrm {V_{m0}=\int_{\Gamma}d\Gamma}\)

  • \(\mathrm {V_{m1}=\int_{\Gamma}\boldsymbol{X}d\Gamma}\)

where:

  • \(\mathrm {\Gamma}\) is the integration domain defined by the volume "under" the triangular panel

  • \(\mathrm {\boldsymbol{X}}\) is the position vector in the Earth-fixed coordinate system

  • \(\mathrm {V_{m0}}\) is the contribution of the panel to the total volume

  • \(\mathrm {V_{m1}}\) is the contribution of the panel to the total volume center

By summing up the contribution of each "contributing panel", the total submerged volume is computed, as well as the position of the center of buoyancy. The two latter quantities are used to compute the total buoyancy forces and moments which are transformed to the body-fixed coordinate system.

2.6. Aerodynamic forces

The wind force is calculated based on the instantaneous wind and body velocities.

2.6.1. Classic wind coefficients

The force is calculated according to the following formula

\[q_j=C_j(\alpha )v^2\]

where:

\(j\)

degree of freedom

\(C\)

wind force coefficient for the instantaneous relative direction

\(v\)

relative wind speed between body and wind

\(\alpha\)

relative velocity direction in local coordinate system

The low frequency body motion is used to calculate the relative wind speed. Optionally, the velocity of the body may be omitted when calculating relative velocity.

2.6.2. Position dependent wind coefficients

This model is adapted for situations where large changes of positions and attitudes are expected during the simulation (a damaged condition for example).

It is possible to define the wind coefficients as a function of:

  • \(\mathrm {\alpha }\) : relative velocity direction in local coordinate system

  • \(\mathrm {z}\) : vertical position of the body in Earth-fixed coordinate system

  • \(\mathrm {\phi }\) : roll angle of the body

  • \(\mathrm {\theta }\) : pitch angle of the body

In that case, the coefficients are interpolated for a given body position and a relative wind direction. The force is computed as in equation Equation (28).

2.7. Wave excitation forces

Wave forces may be split into three contributions:

  • first-order forces that oscillate with wave frequencies

  • second-order mean, rapidly varying, and slowly varying wave drift forces.

  • higher-order ringing forces.

It is assumed that the wave forces can be formally expressed by the following equation. The two first terms are a Volterra series truncated after the quadratic term, see (Hasselmann, 1966) and (Dalzell, 1976), while the last term is the time-dependent ringing load.

\[\begin{array}{l}q(t)=q^{(1)}(t)+q^{(2)}(t)+q^{(R)}(t)\\\\\displaystyle =\frac{1}{2\pi }\int^\infty_{-\infty}h^{(1)}(\tau_1)\zeta(t-\tau_1)\mathrm {d}\tau_1+\frac{1}{4\pi ^2}\int^\infty_{-\infty}\int^\infty_{-\infty}h^{(2)}(\tau_1,\tau_2)\zeta(t-\tau_2)\zeta(t-\tau_1)\mathrm {d}\tau_1\mathrm {d}\tau_2+q^{(R)}(t)\end{array}\]

where:

ringing moment

\(q^{(1)}(t)\)

time-dependent first-order wave force

\(q^{(2)}(t)\)

time-dependent second-order wave force

\(q^{(R)}(t)\)

time-dependent third-order ringing force or third- and fourth-order ringing moment

\(h^{(1)}\)

linear impulse response function

\(h^{(2)}\)

second-order impulse response function

It is assumed that the kernels \(\mathrm {h^{(1)}}\) and \(\mathrm {h^{(2)}}\) are smooth, absolutely integrable and possess Fourier transforms as follows

\[h^{(1)}(\tau)=\frac{1}{2\pi }\int^\infty_{-\infty}H^{(1)}(\omega )e^{i\omega \tau}\mathrm {d}\omega\]
\[H^{(1)}(\omega )=\int^\infty_{-\infty}h^{(1)}(\tau)e^{-i\omega \tau}\]
\[h^{(2)}(\tau_1,\tau_2)=\int^\infty_{-\infty}\int^\infty_{-\infty}H^{(2)}(\omega _1,\omega _2)e^{i(\omega _1\tau_1+\omega _2\tau_2)}\mathrm {d}\omega _1\mathrm {d}\omega _2\]
\[H^{(2)}(\omega _1\omega _2)=\frac{1}{4\pi ^2}\int^\infty_{-\infty}\int^\infty_{-\infty}h^{(2)}(\tau_1,\tau_2)e^{-i(\omega _1\tau_1+\omega _2\tau_2)}\mathrm {d}\tau_1\mathrm {d}\tau_2\]

where:

\(H^{(1)}\) first-order transfer function

\(H^{(2)}\)

second-order transfer function

The physical interpretation of first-order transfer functions is well established in linear response analysis.

The interpretation of the second-order transfer function is not obvious. A simple example taken from (Dalzell, 1976) shows that the second-order response to a dual harmonic excitation with frequencies \(\mathrm {\omega _i}\) and \(\mathrm {\omega _j}\) consists of a mean value in addition to the contribution from sum frequencies \(\mathrm {(\omega _i+\omega _j)}\) and difference frequencies \(\mathrm {(\omega _i-\omega _j)}\).

The magnitude of the second-order forces is relatively small compared to first-order forces, but may be important for structures with natural frequencies below the wave period range. A typical example is resonant heave/roll/pitch motions of TLPs. The second-order forces are usually of importance only if any of the difference or sum frequencies coincide with natural frequencies of the system which are outside the frequency range of the first-order forces.

The ringing response of TLPs with high natural heave and roll/pitch periods \((3-5~\mathrm {s})\) or large GBSs in deep water with high fundamental modal period \((3.5-6~\mathrm {s})\) has been shown to introduce extreme dynamic loads causing major concerns for the platform designs. From several studies one has concluded that a sufficiently accurate ringing load model is the most critical part of the modelling of the ringing phenomenon. Ringing is related to response in irregular seas.

2.8. First-order wave excitation forces

First order wave forces are described in the frequency domain as a transfer function between wave elevation and force

\[\boldsymbol{q}_{WA}^{(1)}(\omega )=\boldsymbol{H}^{(1)}(\omega )\zeta(\omega )\]

where:

\(\boldsymbol{H}^{(1)}\)

complex first order transfer function

\(\tilde{\zeta} (\omega)\)

complex harmonic wave component

2.9. Simplified second-order wave force model

The Newman method can be used to calculate second-order wave forces. This method is outlined in (Newman, 1974). (Faltinsen and Løken, 1979) and (Pinkster, 1981) compared Newman’s simplified approach with the different contributions from the second-order transfer function.

The second-order force can be expressed in the form

\[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_nH_{mn}^{(2+)}e^{i(\omega _m+\omega _n)t}\}+\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}e^{i(\omega _m-\omega _n)t}\}\]

where:

\(\tilde{\zeta}_m\)

complex Fourier component of sea surface elevation with frequency \(\omega_m\)

\(H_{mn}^{(2+)} \equiv H^{(2+)} (\omega_m, \omega_n )\)

second-order transfer function for the sum-frequency force

\(H_{mn}^{(2-)} \equiv H^{(2-)} (\omega_m, \omega_n )\)

second-order transfer function for the difference-frequency force

Neglecting the sum-frequency force in equation Equation (36), the slowly-varying second-order force is expressed by

\[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}e^{i(\omega _m-\omega _n)t}\}\]

where the time average is

\[\bar{q}=\frac{1}{2}\mathrm {Re}\{\sum_m\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}\}\]

The main contribution to the slowly-varying force is assumed to be associated with frequency combinations for which, \(\mathrm {\omega _m\approx \omega _n}\) or

\[|\omega _m-\omega _n|<<\frac{1}{2}(\omega _m+\omega _n)\]

Assuming that

\[H_{mn}^{(2-)}\approx H_{mm}^{(2-)}\]

the slowly-varying force can be approximated by

\[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}e^{i(\omega _m-\omega _n)t}\}\]

The double sum is circumvented using the relations

\[\begin{array}{lll}\displaystyle (\mathrm {Re}\{\sum_m\tilde{\zeta}_m\sqrt{|H_{mn}^{(2-)}|}e^{i\omega _mt}\})^2&=&\displaystyle \frac{1}{2}\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n\sqrt{|H_{mm}^{(2-)}H_{nn}^{(2-)}|}e^{i(\omega _m+\omega _n)t}\}\\\\\displaystyle &+&\displaystyle \frac{1}{2}\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*\sqrt{|H_{mm}^{(2-)}H_{nn}^{(2-)}|}e^{i(\omega _m-\omega _n)t}\}\end{array}\]

and

\[\sqrt{|H_{mm}^{(2-)}H_{nn}^{(2-)}|}\approx |H_{mm}^{(2-)}|\]

Thus, the slowly varying force may be written

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

where

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

The overline symbolizes low-pass filtering. These sums are to be calculated only for those terms where the argument of the square root is positive in each case.

(Kaasen, 1999) has shown that a modification to Newman’s method can be done in order to eliminate the high frequency noise generated by this method. By not using pre-generation of force time series, filtering out this noise must be done simultaneously with the simulation where there is a significant challenge to avoid a phase lag that may affect the low frequency components. Hence, avoiding the high frequency noise is a valuable improvement.

By using the difference-frequency force transfer functions in both parts of Equation (35), the force is

\[q(t)=\mathrm {Re}\{\sum_m\sum_nH_{mn}^{(2-)}[\tilde{\zeta}_m\tilde{\zeta}_n^*e^{i(\omega _m-\omega _n)t}+\tilde{\zeta}_m\tilde{\zeta}_ne^{i(\omega _m+\omega _n)t}]\}\]

As in Equation (38) it is assumed that \(\mathrm {H_{mn}^{(2-)}\approx H_{mm}^{(2-)}}\), and two new complex variables are introduced

\[\begin{array}{lll}\displaystyle u(t)=\sum_m\sqrt{H_{mm}^{(2-)}}\mathrm {Re}\{\tilde{\zeta}_me^{i\omega _mt}\}&&\displaystyle v(t)=\sum_m\sqrt{H_{mm}^{(2-)}}\mathrm {Im}\{\tilde{\zeta}_me^{i\omega _mt}\}\end{array}\]

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 as in Equation (40), can then be seen from

\[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{nm}^{(2-)}e^{i(\omega _m-\omega _n)t}\}=\mathrm {Re}\{u(t)^2+v(t)^2\}\]

Generally, the total velocity potential in the fluid is described by the perturbation scheme

\[\Phi=\varepsilon \Phi^{(1)}+\varepsilon ^2\Phi^{(2)}+\dot s\]

\(\mathrm {\varepsilon }\) is a small quantity. \(\mathrm {\Phi^{(1)}}\) is the first-order potential, including the undisturbed incoming wave, wave diffraction and influence from first-order motions in heave, sway and roll.

(Faltinsen and Løken, 1979) and (Pinkster, 1981) have shown that the second-order potential is only important for the off-diagonal terms of \(\mathrm {c_{ij}}\). At low wave frequencies, where first-order diffraction effects are small, the contribution due to second-order velocity potential is of greater importance.

The discrepancy between the forces based on Equation (36) based on Newman’s approximation, Equation (43), is, according to Faltinsen and Løken, within \(\mathrm {20~\%}\) in force amplitudes and almost negligible in the characteristic period. A better frequency resolution would probably reduce the discrepancy even more.

It is possible to pre-calculate wave drift forces for a number of headings before time domain simulation starts, and interpolate between these time series in time domain.

An alternative is to calculate the derivative of the average wave drift force components with respect to wave heading, which is estimated as the difference

\[\frac{\delta \bar{q}}{\delta \alpha }\cong\frac{\big(\bar{q}(\alpha +\Delta \alpha )-\bar{q}(\alpha -\Delta \alpha )\big)}{2\Delta \alpha }\]

In the time domain simulation, correction of forces with respect to heading changes is implemented as

\[q'=q+\frac{\delta \bar{q}}{\delta \alpha }\Delta \alpha\]

This should be used only for small yaw angles, e.g. within \(\mathrm {\pm10\deg}\).

2.9.1. Second order wave forces in short-crested seas

The following method is implemented in order to include the effect of short-crested seas on wave drift forces.

The average drift force in irregular waves is

\[\bar{q}=\int^\infty_02\bar{H}^{(2)}(\omega )S_\zeta(\omega )\mathrm {d}\omega\]

in which the average, or uni-directional drift force coefficient for the direction \(\mathrm {\alpha _1}\) is

\[\bar{H}^{(2)}(\omega )=\int_{\alpha _1-\frac{\pi }{2}}^{\alpha _1+\frac{\pi }{2}}H^{(2)}(\alpha ,\omega )\theta (\alpha -\alpha _1)\mathrm {d}\alpha\]

where \(\mathrm {\theta (\alpha )}\) is the spreading function.

Time series of slowly varying forces are derived from a sampled wave pattern represented by an array of harmonic components.

\[\begin{array}{lllll}\tilde{\zeta}_{kl}(\alpha _k,\omega _l)&&k=1,2,\dot sc,n_k&&l=1,2,\dot sc,n_l\end{array}\]

The surface elevation at the origin is represented by the components \(\mathrm {\tilde{\zeta}_l(\omega _l)}\)

\[\tilde{\zeta}_l(\omega _l)=\sum_{k=1}^{n_k}\tilde{\zeta}_{kl}\]

The drift force due to the component \(\mathrm {\tilde{\zeta}_{kl}}\) is

\[\bar{q}_l=\sum_{k=l}^{n_k}|\tilde{\zeta}_{kl}|^2H^{(2)}(\alpha _k,\omega _l)\]

An average drift force coefficient for the frequency is thus computed as

\[\hat{H}_l^{(2)}=\frac{\displaystyle \sum_{k=l}^{n_k}|\tilde{\zeta}_{kl}|^2H^{(2)}(\alpha _k,\omega _l)}{\displaystyle \sum_{k=l}^{n_k}|\tilde{\zeta}_{kl}|^2}\]

The expected value of this coefficient is equal to \(\mathrm {\bar{H}^{(2)}(\omega _l)}\), but the discrete sampling over directions gives a scatter similar to the scatter of \(\mathrm {|\tilde{\zeta}_l|}\).

Representing the sea by its surface elevation and using the average coefficient \(\mathrm {\hat{H}_l^{(2)}}\), the methods of Hsu and Blenkarn and Newman can be used directly for generating time series of slowly varying wave force. The validity for short-crested seas is unknown and should be subjected to further investigation. The validity also becomes questionable when the frequency dependence of \(\mathrm {H^{(2)}(\omega )}\) is strong, i.e. when \(\mathrm {H^{(2)}(\omega )}\) varies rapidly compared to the wave spectrum \(\mathrm {S_\zeta(\omega )}\).

2.9.2. Second order wave forces calculated with FFT and COS method in short-crested seas

Long crested waves.

The calculation is based on the complex wave Fourier components \(\mathrm {X(\omega )}\) and the real drift force coefficients \(\mathrm {C(\omega )}\) for the actual direction.

The drift force is calculated from the complex coefficients

\[Z(\omega )=X(\omega )\sqrt{C(\omega )}\]

FFT and COSINE simulation method will give the same result.

FFT method short-crested seas.

For each frequency, SIMO will calculate \(\mathrm {Z(\omega )}\) from equivalent coefficients \(\mathrm {X(\omega )}\) and \(\mathrm {C(\omega )}\)

\[\begin{array}{l}C(\omega )=\frac{\displaystyle \sum_{i=1}^{n_{dir}}X_i(\omega )X_i^*(\omega )c_i(\omega )}{\displaystyle \sum_{i=1}^{n_{dir}}X_i(\omega )X_i^*(\omega )}\\\\\displaystyle X(\omega )=\sum_{i=1}^{n_{dir}}X_i(\omega )\end{array}\]

COS-series short-crested seas.

For each frequency, SIMO will calculate an equivalent coefficient \(\mathrm {Z(\omega )}\) as

\[Z(\omega )=\sum_{i=1}^{n_{dir}}X_i(\omega )\sqrt{c_i(\omega )}\]

Depending on the degree of spreading, these two methods will give different results.

2.10. Full QTF second-order wave excitation forces

Generally, the second-order force for multi-directional or short-crested seas may be written

\[q_2=\int_{-\pi }^\pi \int_{-\pi }^\pi \int^\infty_{-\infty}\int^\infty_{-\infty}Z(\omega _1,\beta _1)Z^*(\omega _2,\beta _2)H^{(2)}(\omega _1,\omega _2,\beta _1,\beta _2)\mathrm {d}\omega _1\mathrm {d}\omega _2\mathrm {d}\beta _1\mathrm {d}\beta _2\]

where \(\mathrm {H^{(2)}}\) is the full quadratic transfer function. The equation can be solved directly by substituting all integrals by summations, for example according to procedures by (Stansberg, 1994). Various simplifications are available depending on the amount of data available, for example the Newman approximation instead of a full second order formulation. Directional superposition is not allowed due to bad convergence for this simplification.

2.11. Wave drift damping

Two different methods are possible:

  • a simplified model based on the derivatives of non-dimensional average wave drift forces with respect to relative velocity

  • Newman’s method

2.11.1. Simplified method

The wave drift damping model is based on the derivatives of non-dimensional average wave drift forces with respect to relative velocity.

\[c_{wd}=\frac{\partial \bar{F}_i}{\partial v_{ir}\bar{F}_i}\]

where \(\mathrm {\bar{F}_i}\) is the average force for a wave of unit wave amplitude. The coefficients \(\mathrm {c_{wd}}\) normally depend on wave frequency, and can be estimated from model tests or 3D sink-source analyses.

Assuming that

\[\frac{\partial F(t)}{\partial vF(t)}\approx \frac{\partial \overline{F(t)}}{\partial v\overline{F(t)}}\]

and

\[\frac{\partial \overline{F(\alpha )}}{\partial v\overline{F(\alpha )}}\approx \frac{\partial \overline{F(\alpha +\Delta \alpha )}}{\partial v\overline{F(\alpha +\Delta \alpha )}}\]

the wave drift forces in the time domain simulation are modified according to

\[\hat{F}_i(t)=F_i(t)(1+c_{wd}v_{wr}), i=1,2,6\]

where \(\mathrm {v_{wr}=v_{wc}-v_{wb}}\) is the difference between current and body velocity in the wave direction.

\(\mathrm {c_{wd}}\) is implemented in the program as a function of estimated peak wave period, and is fixed during the time domain simulation.

This is a simplified model that reflects both wave drift damping and wave drift - current interaction in the time domain.

2.11.2. Newman’s method

The wave drift damping force in the three horizontal degrees of freedom is calculated by:

\[\boldsymbol{F}_B=\begin{bmatrix}F_1\\\\F_2\\\\F_6\end{bmatrix}=\begin{bmatrix}B_{11}&B_{12}&B_{16}\\\\B_{21}&B_{22}&B_{26}\\\\B_{61}&B_{62}&B_{66}\end{bmatrix}\begin{bmatrix}v_{1r}\\\\v_{2r}\\\\v_{6r}\end{bmatrix}\]

Where \(\mathrm {B_{ij}}\) is the wave drift damping force coefficients and \(\mathrm {v_{ir}=v_{ic}-v_{ib}}\) where \(\mathrm {v_1}\), \(\mathrm {v_2}\) and \(\mathrm {v_6}\) are the absolute or relative velocities in direction surge, sway and yaw, respectively. Using velocities relative to current is a simplified way to account for the wave-current interaction effect. If the coefficients already contained corrections for wave-current interaction effects, absolute velocities must be used. This is for example the case if coefficients are computed by using the method described in Wave drift damping coefficients calculation.

The wave drift damping force coefficients in the matrix above can be expressed by:

\[B_{ij}(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*B_{mn,ij}e^{i(\omega _m-\omega _n)t}\}\]

where direction (i) due to velocity in direction (j) in bi-chromatic waves with frequencies ( _m ), ( _n ) and directions ( _m ), ( _n )

\(\tilde{\zeta}_m\)

complex Fourier components of sea surface elevation with frequency \(\omega_m\)

\(B_{mn,ij}\)

\(= B_{ij}(\omega_m, \omega_n, \alpha_m, \alpha_n)\) wave drift damping force coefficient in direction \(i\) due to velocity in direction \(j\) in bi-chromatic waves with frequencies \(\omega_m\), \(\omega_n\) and directions \(\alpha_m\), \(\alpha_n\)

The expression above is similar to Equation (36) for the slowly varying second-order force. The calculation method with assumptions are further similar to the Newman model outlined in First-order wave excitation forces for the wave drift force, with the additional assumption \(\mathrm {\alpha _m=\alpha _n}\).

It is possible to calculate wave drift damping forces for a number of headings before time domain simulation starts, and interpolate between these time series in time domain to correct the forces with respect to heading changes.

In SIMA, it is possible to compute the coefficients \(B_{ij}\) from Wave Drift Coefficients or QTF, using Aranha’s formulation described in Wave drift damping coefficients calculation.

2.12. Wave-current correction of second order wave drift force and wave drift damping calculation

The following subsections presents the theory backgound for the calculation of:

  • Second order wave drift force coefficients for a specified current condition (velocity and direction)

  • Wave drift damping coefficients ffor a specified current condition (velocity and direction)

2.12.1. Wave-current interaction correction of full QTFs using extended Aranha method

Consider a stationary floating body subjected to a current speed \(U\) with a going-to direction of \(\theta\) relative to the body x-axis and waves propagating in direction \(\beta\) relative to the body x-axis. The current components in x and y direction then becomes:

\[\begin{equation} \begin{split} u & =U\cos\theta\\ v & =U\sin\theta \end{split} \end{equation}\]

The full QTF with wave-current interaction effects are then taken as,

\[Q_{U}\left(\omega_{1},\omega_{2},\beta_{1},\beta_{2}\right)=B\left(\omega_{1},\omega_{2},\beta_{1},\beta_{2}\right)Q_{U}\left(\omega_{e,1},\omega_{e,2},\beta_{e,1},\beta_{e,2}\right)\]

Here, the frequency of encounter can be written as,

\[\omega_{e,i}=\left(1+\frac{u\cos\left(\beta_{i}\right)+v\sin\left(\beta_{i}\right)}{C_{p,i}}\right)\omega_{i}\]

and the direction of encounter can be written as:

\[\beta_{e,i}=\beta_{i}+\frac{1}{C_{g,i}}\left(-u\sin\left(\beta_{i}\right)+v\cos\left(\beta_{i}\right)\right)\]

Here, \(C_p\) and \(C_g\) is the phase and group velocity the waves (finite or infinite water depth), respectively.

Two alternative scaling factors can be used. The first is based on the aritmetic average of the Aranha scaling factors for the two interacting wave components:

\[B\left(\omega_{1},\omega_{2},\beta_{1},\beta_{2}\right)=\frac{1}{2}\left(A\left(\omega_{1},\beta_{1}\right)+A\left(\omega_{2},\beta_{2}\right)\right)\]

The other is based on the geometric average:

\[B\left(\omega_{1},\omega_{2},\beta_{1},\beta_{2}\right)=\sqrt{A\left(\omega_{1},\beta_{1}\right)A\left(\omega_{2},\beta_{2}\right)}\]

The Aranha scaling factor can be written as,

\[A\left(\omega,\beta\right)=1-\frac{\left(\omega\frac{\partial\alpha}{\partial\omega}-2\right)}{C_{g}}\left(u\cos\left(\beta\right)+v\sin\left(\beta\right)\right)\]

where,

\[\alpha=\frac{C_{g}}{C_{p}}\]

The dispersion relation in the finite water depth case reads,

\[\omega=\sqrt{kg\tanh\left(kh\right)}\]

which gives a phase velocity and group velocity:

\[C_{p}=\frac{\omega}{k}=\sqrt{\frac{g\tanh\left(kh\right)}{k}}\]
\[C_{g}=\frac{d\omega}{dk}=\alpha C_{p}\]
\[\alpha =\frac{C_{g}}{C_{p}}=\frac{1}{2}\left(1+\frac{2kh}{\sinh\left(2kh\right)}\right)=\frac{1}{2}+\frac{kh}{\sinh\left(2kh\right)}\]
\[\begin{equation} \begin{split} \frac{d\alpha}{d\omega} & =\frac{d\alpha}{dk}\frac{dk}{d\omega}\\ & =\frac{1}{C_{g}}\frac{d\alpha}{dk}\\ & =\frac{h}{C_{g}}\frac{\left(1-2kh\coth\left(2hk\right)\right)}{\sinh\left(2kh\right)}\\ & =\frac{h}{C_{g}\sinh\left(2kh\right)}\left(1-\frac{2kh}{\tanh\left(2kh\right)}\right) \end{split} \end{equation}\]

2.12.2. Wave drift damping coefficients calculation

Using Equation (73) and introducing the notation \(\kappa=1-\frac{\omega}{2}\frac{\partial(\frac{C_{g}}{C_{p}})}{\partial\omega}\) , we express Aranha’s formula for second order wave drift force as:

\[\mathbf{F}(\omega,\beta,\mathbf{U}_{r})=(1+2\frac{κ}{C_{g}}(u_{r}\cos\left(\beta\right)+v_{r}\sin\left(\beta\right)))\mathbf{F}(\omega_{e},\beta_{r},\mathbf{0})\]

where \(\mathbf{U}_{r}=\begin{bmatrix} u_{r} \\ v_{r} \end{bmatrix} = \begin{bmatrix} u_{0}-u_{b} \\ v_{0}-v_{b} \end{bmatrix}\) is the relative velocity between the current and the body low frequency velocity \(\mathbf{U}_{b}=\begin{bmatrix} u_{b}\\ v_{b} \end{bmatrix}\)

Using Taylor expansion at first order, about the current velocity \(\mathbf{U}_{r}=\begin{bmatrix} u_{0}\\ v_{0} \end{bmatrix}\) :

\[\mathbf{F}(\omega,\beta,\mathbf{U}_{r})=\mathbf{F}(\omega,\beta,\mathbf{U}_{0})-\frac{\partial\mathbf{F}}{\partial\mathbf{U}_{r}}(\omega,\beta,\mathbf{U}_{0})\mathbf{U}_{b}\]

The term \(-\frac{\partial\mathbf{F}}{\partial\mathbf{U}_r}(\omega,\beta,\mathbf{U}_{0})\mathbf{U}_{b}\) represents the wave drift damping that can be interpreted as the change of the drift force due to the low frequency velocity of the vessel. It is possible to calculate the coefficient \(\mathbf{B}(\omega,\beta,\mathbf{U}_0)=\frac{\partial\mathbf{F}(\omega,\beta,\mathbf{U}_r)}{\partial\mathbf{U}_r}_{|\mathbf{U}_r=\mathbf{U}_0}\) by using Equation (80) only for surge, sway and yaw (respectively noted as degrees of freedom no.1,2 and 6):

\[\begin{equation} \begin{aligned} \mathbf{B}(\omega,\beta,\mathbf{U}_0) & = \begin{bmatrix} B_{11} & B_{12}\\ B_{21} & B_{22}\\ B_{61} & B_{62} \end{bmatrix} \\ & = \frac{\partial\mathbf{F}(\omega,\beta,\mathbf{U}_r)}{\partial\mathbf{U}_r}_{|\mathbf{U}_r=\mathbf{U}_0} \\ & =P(\mathbf{U}_r)\mathbf{F}(\omega_{e},\beta_{r},\mathbf{0})\left[\begin{array}{cc} \frac{\partial}{\partial u_r} & \frac{\partial}{\partial v_r}\end{array}\right]_{\begin{array}{|c} u_r=u_{0}\\ v_r=v_{0} \end{array}} \end{aligned} \end{equation}\]

where \(P(\mathbf{U}_r)=1+2\frac{\kappa}{C_{g}}(u_r\cos\left(\beta\right)+v_r\sin\left(\beta\right))\)

After some calculation, we obtain:

\[\frac{\partial\mathbf{F}(\omega,\beta,u_r,v_r)}{\partial u_r}_{\begin{array}{|c} u_r=u_{0}\\ v_r=v_{0} \end{array}}=W_{c0}+P(\mathbf{U}_0)W_{u_r}\]

and

\[\frac{\partial\mathbf{F}(\omega,\beta,u_r,v_r)}{\partial v_r}_{\begin{array}{|c} u_r=u_{0}\\ v_r=v_{0} \end{array}}=W_{s0}+P(\mathbf{U}_0)W_{v_r}\]

where:

\[W_{c0}=2\frac{\kappa}{C_{g}}\cos\left(\beta\right)\mathbf{F}(\omega_{e_{0}},\beta_{r_{0}},0,0)\]
\[W_{s0}=2\frac{\kappa}{C_{g}}\sin\left(\beta\right)\mathbf{F}(\omega_{e_{0}},\beta_{r_{0}},0,0)\]
\[P(\mathbf{U}_0)=1+2\frac{\kappa}{C_{g}}(u_{0}\cos\left(\beta\right)+v_{0}\sin\left(\beta\right))\]
\[W_{u_r}=\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\omega}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}k\cos\left(\beta\right)-\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\beta}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}\frac{1}{C_{g}}\sin\left(\beta\right)\]
\[W_{v_r}=\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\omega}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}k\sin\left(\beta\right)+\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\beta}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}\frac{1}{C_{g}}\cos\left(\beta\right)\]

where \(\omega_{e_{0}}\) and \(\beta_{r_{0}}\) are respectively the encounter frequency and the encounter heading for the body at rest (i.e. for \(\mathbf{U}_b=\mathbf{0}\)).

The wave drift damping coefficients \(\mathbf{B}(\omega,\beta,\mathbf{U}_0)\) can thus be computed from the following inputs:

  • wave drift force coefficients at zero current velocity as function of wave direction and wave frequency (note: when using the QTF as input, the wave drift force coefficients are extracted from the diagonal)

  • wave directions

  • wave frequencies

  • water depth

  • current velocity

  • current direction

The obtained wave drift damping coefficients are:

They represent the change of wave drift forces in X and Y and moment in Z due to a velocity in X and Y.

2.13. Nonlinear modification of restoring and wave forces

In SIMO the model for wave-induced motion is based on some simplifying assumptions, in particular linearization of the hydrostatic, diffraction and radiation forces about the still water plane and the body’s hydrostatic equilibrium.

To partly overcome this deficiency, SIMO offers a way to account for the nonlinearity in the restoring force and the wave force in a quasistatic manner during the simulation. The basic idea is simple and straightforward: The body surface is divided into 4-node panel cells. At each time step, the location of each cell centroid is calculated in earth-fixed coordinates. The wave elevation at this cell centroid is also evaluated. If the cell centroid is above the mean waterline when the body is in its equilibrium position and below the wave surface at any time step during the simulation, the 3D local nonlinear correction in the restoring and wave forces is formed as the product of the panel area and the quasistatic pressure applied on the panel cell (mass density of the water \(\mathrm {\times }\) acceleration of gravity \(\mathrm {\times }\) vertical distance between the wave surface and the cell centroid). The local nonlinear correction is also included if the cell centroid is below the mean waterline when the body is in the equilibrium position, but above the wave surface at any time step during the simulation. The local nonlinear corrections are integrated to give the total nonlinear correction in the 6 degree-of-freedom generalized rigid body forces and moments in the body-fixed coordinate system.

Nonlinear correction of wave force and restoring force on body. Quasistatic pressure is computed over the wetted panels, and the resulting time-varying forces and moments are added as corrections to those obtained using frequency-dependent first-order excitation functions and wave drift functions.

The panel model need not cover the entire wetted hull, but only the part that will be partly wet or dry, see Nonlinear correction of wave force…​. The user specifies a minimum draft and a maximum draft for the panel model, which is read into SIMO from a text file of type GDF, which is used by WAMIT. Using few panels will reduce the computational time.

The wave elevation outside the hull can be that of the incoming, undisturbed wave or that corresponding to the diffracted wave based on first-order theory, if available.

However, the above-mentioned nonlinear modification has been partially accounted for in the WAMIT second order theory. This part of WAMIT second order force has to be removed in order to avoid double counting. The total second order hydrodynamic force from WAMIT can be written as

\[\boldsymbol{F}^{(2)}=\iint_{S_0}p^{(1)}\boldsymbol{n}^{(1)}\mathrm {d}S+\iint_{S_0}p^{(2)}\boldsymbol{n}^{(0)}\mathrm {d}S+\iint_{S_0}p^{(0)}\boldsymbol{n}^{(2)}\mathrm {d}S+\iint_{S-S_0}p^{(1)}\boldsymbol{n}^{(0)}\mathrm {d}S\]
\[\boldsymbol{M}^{(2)}=\iint_{S_0}p^{(1)}\boldsymbol{r}\times \boldsymbol{n}^{(1)}\mathrm {d}S+\iint_{S_0}p^{(2)}\boldsymbol{r}\times \boldsymbol{n}^{(0)}\mathrm {d}S+\iint_{S_0}p^{(0)}\boldsymbol{r}\times \boldsymbol{n}^{(2)}\mathrm {d}S+\iint_{S-S_0}p^{(1)}\boldsymbol{r}\times \boldsymbol{n}^{(0)}\mathrm {d}S\]

where \(\mathrm {p^{(0)}}\), \(\mathrm {p^{(1)}}\) and \(\mathrm {p^{(2)}}\) are the hydrostatic water pressure, first order and second order hydrodynamic water pressures while \(\mathrm {\boldsymbol{n}^{(0)}}\), \(\mathrm {\boldsymbol{n}^{(1)}}\) and \(\mathrm {\boldsymbol{n}^{(2)}}\) represent the mean, first order and second order normal vectors of the body surface expressed in the earth-fixed coordinate system. \(\mathrm {\boldsymbol{r}}\) is the position vector. \(\mathrm {S_0}\) denotes the mean wetted surface. \(\mathrm {S-S_0}\) is the difference between the instantaneous wetted surface and the mean wetted surface. The last terms in Equation (90) and Equation (91) can be expressed as

\[\iint_{S-S_0}p^{(1)}\boldsymbol{n}^{(0)}\mathrm {d}S=\frac{1}{2}\rho g\int_{WL}\boldsymbol{n}^{(0)}[\zeta^{(1)}_r]^2\{1-[n_z^{(0)}]^2\}^{-1/2}\mathrm {d}l\]
\[\iint_{S-S_0}p^{(1)}\boldsymbol{r}\times \boldsymbol{n}^{(0)}\mathrm {d}S=\frac{1}{2}\rho g\int_{WL}\boldsymbol{r}\times \boldsymbol{n}^{(0)}[\zeta^{(1)}_r]^2\{1-[n_z^{(0)}]^2\}^{-1/2}\mathrm {d}l\]

where \(\mathrm {\rho }\) and \(\mathrm {g}\) are the mass density of water and the acceleration of gravity, respectively. \(\mathrm {WL}\) stands for waterline. \(\mathrm {\zeta_r^{(1)}}\) is the first order relative wave elevation. \(\mathrm {n_z^{(0)}}\) represents the vertical component of the mean normal vector \(\mathrm {\boldsymbol{n}^{(0)}}\). The removal of Equation (92) and Equation (93) is taken care of in the program. If second-order terms are not included in the simulation these terms are not removed.

In the SIMO input manual, the nonlinear correction described above is referred to as "Nonlinear Buoyancy Correction"

2.14. Diffracted waves

Diffracted wave transfer functions for several diffracted wave points can be given for a body. These transfer functions are given in a similar manner as e.g. the wave to motion transfer functions. Diffracted wave transfer functions for wave elevation and global x-, y- and z-velocity may be given. If these TRF’s are included in a simulation, time series of diffracted waves may be generated, but the TRF’s will not affect the simulation results. If a TRF of velocity is specified, time series of both velocity and acceleration will be generated.

A body may be in the "wave shadow" of another, larger body. This effect can be included in SIMO by data group "Diffracted wave application" The body will then "feel" the diffracted waves, instead of the original "undisturbed" wave field.

There are two ways to specify that a body shall "feel" diffracted waves.

  1. For body types 1 or 2 (large volume bodies), the user may specify one or more "slender elements" and "fixed body elements" under data group DISTRIBUTED ELEMENT FORCE. Each of these elements may feel different diffracted waves by specifying a "body name" and "diffracted wave points". For slender elements, each end of the element refers to one "diffracted wave point". If the point numbers are identical, all parts of the element will "feel" the same wave particle motion series. Otherwise, linear interpolation is performed on the particle motion series. The special case of "fixed body elements" refers to one "diffracted wave point" only.

  2. For body type 3 (3-dof body), the effect of wave diffraction may be included in data group SMALL BODY HYDRODYNAMIC DATA by referring to a body and a "diffracted wave point".

Diffracted waves are not included in other force models. This means that in order to include diffraction effects in these models, e.g. "Second order wave drift coefficients" the coefficients have to be modified.

2.15. Current drag forces

2.15.1. Classic current coefficients

All body current forces are computed using the current velocity at the surface (\(\mathrm {z=0}\)). Traditionally, the viscous hull surge and sway force and yaw moment have been calculated based on current coefficients and the instantaneous magnitude of the translational relative velocity between the vessel and the fluid. The current drag forces are then expressed by:

\[q_{CU}^{(k)}(\alpha ,t)=C_1^{(k)}(\alpha )|u(t)|+C_2^{(k)}(\alpha )|u(t)|^2\]
\[|u|^2=(v_1-\dot {x}_1)^2+(v_2-\dot {x}_2)^2\]
\[\alpha =\arctan\frac{v_2-\dot {x}_2}{v_1-\dot {x}_1}\]

where:

velocity

current velocity

\(k\)

degree of freedom

\(C_1\)

linear current force coefficient

\(C_2\)

quadratic current force coefficient

\(u\)

relative velocity between low frequency body velocity and current velocity

\(\alpha\)

relative angle between direction of low frequency body velocity and current velocity

\(v_1,v_2\)

current velocity components at the surface

\(\dot{x}_1,\dot{x}_2\)

vessel velocity components in the vessel coordinate system

image323
Figure 2. Separation of cross-flow due to translation and rotation.

Usually, the linear coefficient, \(\mathrm {C_1}\) is not used. The model does not account for the effect from the yaw-induced cross-flow. This cross-flow has normally been included as a separate yaw damping term. The cross-flow based on a pure relative translational velocity and yaw velocity is illustrated in Separation of cross-flow…​.

Assuming that a quadratic yaw damping term has been included, the viscous yaw moment may be calculated as the sum of the contribution from the translational velocity and the yaw velocity

\[q^{(6)}(\alpha )=C_2^{(6)}(\alpha )|u|^2+b_{66}\dot {x}_6|\dot {x}_6|\]

where:

\(C_2^{(6)}(\alpha)\)

Current coefficient in yaw as a function of the incident heading, \(\alpha\)

\(b_{66}\)

Quadratic damping coefficient in yaw due to a pure yaw velocity

\(\dot{x}_6\)

Yaw velocity

Normally, the current coefficient in yaw, \(\mathrm {C_2^{(6)}(\alpha )}\), also includes the Munk moment which is an effect that can be derived from potential theory. When quadratic current coefficients are included in the simulation model, the Munk moment is excluded from the equations of motion. The Munk moment, denoted \(\mathrm {q_M(\alpha )}\), is given as

\[q_M(\alpha )=-\frac{1}{2}|u|^2[A_2-A_1]\sin(2\alpha )\]

where:

\(A_2\) Added mass in sway

\(A_1\)

Added mass in surge

It is, however, quite clear that the current velocity, the vessel translational velocity and the yaw velocity will influence the distribution of the cross-flow along the hull. Theoretically, the sum of these contributions will result in a linear cross-flow distribution along the hull.

Due to the quadratic nature of the viscous hull forces, the forces obtained from the vessel translation and the current and the forces obtained from a yaw-induced cross-flow can’t be separated and then added. Therefore a force model integrating up the sway force and yaw moment along the vessel due to the resulting instantaneous cross-flow has been established as strip model no 1 and strip model no. 2.

In order to get simple, analytical equations of the viscous sway force and yaw moment, the drag properties are assumed to be constant along the vessel.

For strip model no. 1:

\[C_d(x)=\frac{1}{2}\rho c_d(x)T(x)=\frac{C_2^{(2)}(90\deg)}{l}\]

For strip model no. 2:

\[C_d(x)=\frac{1}{2}\rho c_d(x)T(x)=\frac{C_2^{(2)}(\alpha )}{l\sin^2(\alpha )}\]

where:

\(c_d(x)\)

Non-dimensional drag coefficient of the vessel section

\(C_d(x)\)

Drag coefficient of the vessel section

\(T(x)\)

Vessel draught at the vessel section

\(l\)

Vessel length

\(\rho\)

Water density

The drag force of a strip, \(\mathrm {f_s(x)}\), located at a distance from the centre of motion is then given as

\[f_s(x)=C_d(x)(u_{r2}-x\dot {x}_6)^2\]

where

\[u_{r2}=|\boldsymbol{u}|\sin(\alpha )\]

The total sway force, \(\mathrm {q^{(2)}(\alpha )}\), is then given as

\[q^{(2)}(\alpha )=\int_{-l_2}^{l_1}C_d(x)(u_{r2}-x\dot {x}_6)^2\mathrm {d}x\]

while the yaw moment is given as

\[q^{(6)}(\alpha )=q_v^{(6)}(\alpha )+q_M(\alpha )=\int_{-l_2}^{l_1}C_d(x)(u_{r2}-x\dot {x}_6)^2x\mathrm {d}x+q_M(\alpha )\]

Generally, the two alternative models give increased damping moments and sway forces compared to the separated motion model.

The two-dimensional current coefficient as applied for strip model no.1 may deviate from the current coefficient for the actual heading of the relative velocity, leading to inaccurate quasistatic current sway force. Strip model no. 2, however, gives correct quasistatic current force, but may give numerical problems for small values of \(\mathrm {\sin(\alpha )}\).

When using the strip models in the program, the quadratic current coefficients in yaw are not used, and the Munk moment will be included in the equations of motion.

2.15.2. Position dependent quadratic current coefficients

This model is adapted for situations where large changes of positions and attitudes are expected during the simulation (a damaged condition for example). This model is only available for quadratic current drag coefficients (not available for linear current coefficients).

It is possible to define the current coefficients as a function of: - \(\mathrm {\alpha }\) : relative velocity direction in local coordinate system - \(\mathrm {z}\) : vertical position of the body in Earth-fixed coordinate system - \(\mathrm {\phi }\) : roll angle of the body - \(\mathrm {\theta }\) : pitch angle of the body

In that case, the coefficients are interpolated for a given body position and a relative current direction. The force is computed as in equation Equation (94).

The Munk moments are all assumed to be included in the quadratic current coefficients when using the advanced current coefficients formulation.

Calculations using strip theory are not available when position dependent coefficients are specified.

2.16. Distributed element force

2.16.1. Physical relevance

The distributed element force model is applicable to two different modelling features:

  • Long, slender elements

  • Concentrated, fixed elements (with zero extension)

Long, slender elements can be used to model jacket legs and bracings or spool pieces. A slender object such as a spool piece will normally have a 3D geometry consisting of a number of straight elements with different orientation, and can be modelled by a set of slender elements. End connectors can be modelled either by a short element or by a concentrated, fixed element.

The hydrodynamic properties of a spoolpiece can be assumed to be rotation symmetric. However, in order to handle partly covered or stiffened spoolpieces and other general long structures, an assumed rotation symmetry or assumed cylindrical shape would be too limiting for the model.

A fixed connection between a body and a small element, which can be assumed to give concentrated hydrodynamic forces, is useful for the modelling of various configurations, such as:

  • rigid connection between a framework-type seabed structure and four mudmat structures, one in each corner, for analysis of dynamics in the splash zone or during seabed penetration. The framework can be modelled as a 6-dof body and the mudmats as 3-dof bodies with position-dependent coefficients.

  • imaginary bodies attached to a floating vessel, in order to include non-linear deviations from force coefficients calculated by linear wave diffraction theory programs.

As an alternative model, a small 3-dof body can be rigidly connected to a large body by a set of springs having sufficiently high stiffness to avoid excitation of internal resonance. However, this will need extremely short time steps in the simulation.

2.16.2. Element definition

Both element types may have any orientation. The main body must be of BODY TYPE 1, i.e. 6-dof body.

Small-body theory is used to calculate forces both for the slender element and the concentrated fixed element. Stiff connection implies that all the "sub-body" (slender element or fixed, concentrated element) forces are calculated and transferred to the "main" body. Each slender element specified by the user is divided into NSTRIP strips with equal length.

2.16.3. Coordinate systems

Three coordinate systems are used,

  • XG, global coordinate system, XG(1), XG(2), XG(3)

  • XB, local body fixed coordinate system, XB(1), XB(2), XB(3)

  • XS, local strip coordinate system, XS(1), XS(2), XS(3)

All coordinate systems are orthogonal and right-handed. For XG and XB, see General Assumptions and Notation.

image337
Figure 3. Definition of local strip coordinate system.

Figure 3 illustrates a single strip, introducing the local strip coordinate system. The XS(1)-axis is directed along the longitudinal axis of the strip, and hence also the element, the XS(2)-axis is placed perpendicular to the XS(1)-axis in the local strip xy-plane and the XS(3)-axis is then defined from the definition of orthogonal right-handed coordinate systems.

2.16.4. Calculation of the mass matrix, M

Contribution from structural mass

The contribution from each of the slender elements to the mass matrix of the main body is calculated and added, in local body coordinates. The contribution from each element is calculated by adding the contributions from each strip, which has length \(\mathrm {L}\), mass per unit length, \(\mathrm {m}\), and coordinates \(\mathrm {(x,y,z)}\) in the local body coordinate system. The contribution from each strip within each element to the total mass matrix of the body will be

\[-\Delta \boldsymbol{M}(t) = \begin{bmatrix} 1&0&0&0&z&-y\\ 0&1&0&-z&0&x\\ 0&0&1&y&-x&0\\ 0&-z&y&y^2+z^2&-xy&-zx\\ z&0&-x&-xy&z^2+x^2&-yz\\ -y&x&0&-zx&-yz&x^2+y^2 \end{bmatrix} \cdot mL = \begin{bmatrix} \boldsymbol{I}&\boldsymbol{A}^T\\ \boldsymbol{A}&\boldsymbol{A}\boldsymbol{A}^T \end{bmatrix}\cdot mL\]

where \(\mathrm {\boldsymbol{I}}\) is a 3x3 unit matrix and

\[\boldsymbol{A}=\begin{bmatrix}0&-z&y\\\\z&0&-x\\\\-y&x&0\end{bmatrix}\]

The mandatory specified BODY MASS DATA for the main body can be zero (but not necessarily), and will be added to the contributions from the slender elements.

Contribution from added mass

The contribution from the added mass will be an extension to the above, since the added mass of each strip is different in the 3 directions.

Assume the following parameters:

\(\boldsymbol{V}\) Linear acceleration vector (main body)

\(\boldsymbol{W}\)

Rotation acceleration vector (main body)

\(\boldsymbol{R}\)

Radius vector to the strip

\(\boldsymbol{M}_a\)

Added mass matrix for the strip

\[\begin{array}{lllllll}\boldsymbol{V}=\begin{bmatrix}e\\f\\g\end{bmatrix}&&\boldsymbol{W}=\begin{bmatrix}p\\q\\r\end{bmatrix}&&\boldsymbol{R}=\begin{bmatrix}x\\y\\z\end{bmatrix}&&\boldsymbol{M}_a=\begin{bmatrix}a_x&0&0\\0&a_y&0\\0&0&a_z\end{bmatrix}\end{array}\]

The local acceleration of the strip will be

\[\boldsymbol{V}_L=\boldsymbol{V}+\boldsymbol{W}\times \boldsymbol{R}\]

The local force and moment can be expressed by

\[\begin{array}{l}\boldsymbol{F}_L=\boldsymbol{M}_a\cdot \boldsymbol{V}+\boldsymbol{M}_a\cdot (\boldsymbol{W}\times \boldsymbol{R})\\\\\boldsymbol{M}_L=\boldsymbol{R}\times \boldsymbol{F}_L=\boldsymbol{R}\times (\boldsymbol{M}_a\cdot \boldsymbol{V})+\boldsymbol{R}\times \big(\boldsymbol{M}_a\cdot (\boldsymbol{W}\times \boldsymbol{R})\big)\end{array}\]

Further expansion yields

\[\begin{array}{l}\boldsymbol{A}\cdot \boldsymbol{V}=\boldsymbol{R}\times (\boldsymbol{M}_a\cdot \boldsymbol{V})=\begin{bmatrix}-a_yzf+a_zyg\\a_xze-a_zxg\\-a_xye+a_yxf\end{bmatrix}=\begin{bmatrix}0&-a_yz&a_zy\\a_xz&0&-a_zx\\-a_xy&a_yx&0\end{bmatrix}\cdot \begin{bmatrix}e\\f\\g\end{bmatrix}\\\\\boldsymbol{B}\cdot \boldsymbol{W}=\boldsymbol{R}\times \big(\boldsymbol{M}_a\cdot (\boldsymbol{W}\times \boldsymbol{R})\big)=\begin{bmatrix}a_zy^2+a_yz^2&-a_zxy&-a_yxz\\-a_zxy&a_xz^2+a_zx^2&-a_xyz\\-a_yxz&-a_xyz&a_yx^2+a_xy^2\end{bmatrix}\cdot \begin{bmatrix}p\\q\\r\end{bmatrix}\end{array}\]

The resulting added mass contribution is

\[\Delta \boldsymbol{M}_a=\begin{bmatrix}\boldsymbol{M}_a&\boldsymbol{A}^T\\\boldsymbol{A}&\boldsymbol{B}\end{bmatrix}\]

The contribution from all strips will be added to the body mass matrix.

Depth-dependent hydrodynamic coefficient corrections may be specified for both slender and fixed elements. Added mass data for the main body may be zero (but not necessarily), and will be added to the contributions from the slender elements.

2.16.5. Calculation of external load, F

The model for external load on a slender element strip consists of three contributions:

\(\boldsymbol{F}_B\) buoyancy forces

\(\boldsymbol{F}_W\)

wave forces

\(\boldsymbol{F}_S\)

slamming forces

The resulting force on a slender element is accumulated force contribution from each strip. The body force is the sum of contributions from all elements.

Gravity forces

The gravity force acts in the global z-direction and is written

\[\boldsymbol{F}_{\mathrm {grav},G}=\begin{bmatrix}0\\0\\-mg\mathrm {d}S\end{bmatrix}\]
Buoyancy forces

The buoyancy force acts in the global z-direction, through the center of buoyancy. The load vector for a strip, expressed in the global coordinate system, can be written,

\[\boldsymbol{F}_{\mathrm {buoy},G}=\begin{bmatrix}0\\0\\\rho V_sg\mathrm {d}S\end{bmatrix}\]
Wave forces

The wave force on a strip acts through the center of buoyancy. The load vector for a strip, expressed in the local strip coordinate system can be written

\[\begin{equation} \boldsymbol{F}_{W,S} = \rho V_S\boldsymbol{a}_S + \boldsymbol{m}_h \circ \boldsymbol{a}_S + \boldsymbol{F}_{q,S} + \boldsymbol{C}_l \circ \boldsymbol{U}_{rel,S} \end{equation}\]

where the quadratic drag force is computed as,

\[\begin{equation} \boldsymbol{F}_{q,S} = \left[\begin{array}{c} C_{qx}U_{rel,sx}\left|U_{rel,sx}\right|\\ C_{qy}U_{rel,sy}\sqrt{U_{rel,sy}^{2}+U_{rel,sz}^{2}}\\ C_{qz}U_{rel,sz}\sqrt{U_{rel,sy}^{2}+U_{rel,sz}^{2}} \end{array}\right] \end{equation}\]

and,

\[\begin{equation} \boldsymbol{U}_{rel,S} = \boldsymbol{U}_{S}+\boldsymbol{v}_{S}-\boldsymbol{\dot{X}}_{S} \end{equation}\]

and we have that:

\(V_S\)

is submerged volume per length unit, calculated to \(z = 0\)

\(\boldsymbol{a}_S =\begin{bmatrix} a_{sx} & a_{sy} & a_{sz} \end{bmatrix}^T\)

is water particle acceleration components in local strip coordinate system

\(\boldsymbol{U}_{rel,S} = \begin{bmatrix} U_{rel,sx} & U_{rel,sy} & U_{rel,sz} \end{bmatrix}^T\)

is the relative water particle velocity components in local strip coordinate system (relative to strip velocity)

\(\boldsymbol{v}_S =\begin{bmatrix} v_{sx} & v_{sy} & v_{sz} \end{bmatrix}^T\)

is wave induced water particle velocity components in local strip coordinate system

\(\boldsymbol{\dot{X}}_S =\begin{bmatrix} \dot{X}_S & \dot{Y}_S & \dot{Z}_S \end{bmatrix}^T\)

is strip velocity components in local strip coordinate system

\(\boldsymbol{U}_S =\begin{bmatrix} U_{sx} & U_{sy} & U_{sz} \end{bmatrix}^T\)

is current flow velocity in local strip coordinate system

and \(\qquad \qquad \quad \circ\)

denotes elementwise multiplication

The first and second terms in equation Equation (114) contains the Froude-Krylov and diffraction force. The third term is the quadratic drag term of the Morison formula. The fourth term represents linear drag.

Slamming forces

The slamming force on a moving strip, can be related to the change in (added) mass with time.

Expressed in local strip coordinates, the slamming force can be expressed by

\[\boldsymbol{F}_{S,S} =-\frac{\partial \boldsymbol{m}_s}{\partial t}\boldsymbol{U}_{rel,S} =-\frac{\partial \boldsymbol{m}_s}{\partial h}\frac{\partial h}{\partial t}\boldsymbol{U}_{rel,S} =-\frac{\partial }{\partial h}\begin{bmatrix}m_{h,x(h)}&0&0\\0&m_{h,y(h)}&0\\0&0&m_{h,z(h)}\end{bmatrix}\frac{\partial h}{\partial t}\boldsymbol{U}_{rel,S}\]

where:

\(h\)

distance between instantaneous surface elevation and strip origin in global z-direction

2.16.6. Summary of calculation options

Two parameters are used to define different calculation options:

  • IVOL

    • \(\mathrm {=0}\): Gravity force and buoyancy force are omitted

    • \(\mathrm {=1}\): Gravity force and buoyancy force are included

  • IFOADD

    • \(\mathrm {=1}\): Forces integrated to wave surface: \(\mathrm {Z=\zeta}\)

    • \(\mathrm {=0}\): Forces integrated to still water level: \(\mathrm {Z=0}\)

A summary of the different force components and how they are integrated, is given in Table 4.1.

Table 4.1 Summary of force components
Force component IFOADD \(\mathrm {=1}\) IFOADD \(\mathrm {=0}\)

Inertia: \(\mathrm {F=(m+m_h)\ddot {x}}\)

\(\mathrm {a}\) to \(\mathrm {z=\zeta}\)

\(\mathrm {a}\) to \(\mathrm {z=\zeta}\)

Froude-Krylov: \(\mathrm {F=(\rho V+m_h)\ddot {\zeta}}\)

longitudinal

\(\mathrm {a,V}\) to \(\mathrm {z=0}\)

\(\mathrm {a,V}\) to \(\mathrm {z=0}\)

transverse

\(\mathrm {a,V}\) to \(\mathrm {z=\zeta}\)

\(\mathrm {a,V}\) to \(\mathrm {z=0}\)

Slamming: \(\displaystyle F=v_r\frac{\mathrm {d}m_h}{\mathrm {d}t}\)

\(\mathrm {a}\) to \(\mathrm {z=\zeta}\)

\(\mathrm {a}\) to \(\mathrm {z=\zeta}\)

Drag: latexmath:[\mathrm {F=C_lv_r+C_q

v_r

v_r}]

\(\mathrm {C_1,C_q}\) to \(\mathrm {z=\zeta}\)

\(\mathrm {C_1,C_q}\) to \(\mathrm {z=0}\)

Froude-Krylov adjustment: \(\mathrm {F=\rho gV}\)

\(\mathrm {V;z=0}\) to \(\mathrm {\zeta}\)

\(\mathrm {V;z=0}\) to \(\mathrm {\zeta}\)

Buoyancy \(\mathrm {F=\rho gV}\)

`Buoyancy' from \(\mathrm {z=0}\) to \(\mathrm {z=\zeta}\) is a correction to the Froude-Krylof force, for elements crossing the surface.

2.16.7. Wind force on slender element

The aerodynamic force on a slender element is calculated strip-wise. The element is divided into a number of slices, as indicated in Wind force on slender element…​. The fineness of the division is chosen by the user. It is assumed that the cross sectional geometry is constant and does not vary along the X-axis. Cartesian axes XYZ are defined as shown in Wind force on slender element…​.

Let wind load coefficients \(\mathrm {C_x^{wi}}\), \(\mathrm {C_y^{wi}}\) and \(\mathrm {C_z^{wi}}\) be associated with the three axial directions. Let \(\mathrm {v_x^{wi}}\), \(\mathrm {v_y^{wi}}\) and \(\mathrm {v_z^{wi}}\) be the components of relative wind velocity at the centre of one slice. Here relative wind velocity means the hypothetical wind speed at the location of the slice centre (i.e. the wind velocity that would be if there was no slender element) minus the velocity of the slice centre. The three force components on the slice are then modelled as

\[\begin{array}{l}\Delta F_x^{wi}=C_x^{wi}\:|v_x^{wi}|\:v_x^{wi}\Delta x\\\\\Delta F_y^{wi}=C_{xy}^{wi}\:v_{yz}^{wi}\:v_y^{wi}\Delta x\\\\\Delta F_z^{wi}=C_z^{wi}\:v_{yz}^{wi}\:v_z^{wi}\Delta x\end{array}\]

where \(\mathrm {\Delta x}\) is the length of the slice and \(\mathrm {v_{yz}^{wi}}\) is the component of the relative wind velocity in the YZ plane.

\[v_{yz}^{wi}=\sqrt{(v_y^{wi})^2+(v_z^{wi})^2}\]

The forces from all the strips are added to give the total force on the slender element. The corresponding components of force moment are obtained by multiplication by the respective moment arms from the defined point of reference, which may not be located as shown in Wind force on slender element…​.

Wind force on slender element. Decomposition of relative wind velocity vector shown.

Note that the coefficients \(\mathrm {C_x^{wi}}\), \(\mathrm {C_y^{wi}}\) and \(\mathrm {C_z^{wi}}\) are distributed coefficients with unit \([\mathrm {kN~\frac{(s/m)^2}{m}}]\) or equivalent.

The element is allowed to penetrate the water surface. In this case the aerodynamic force is computed on the dry part of the element. If a slice is partly above the water and partly below, the computation of relative wind velocity, force and moment is based on the mid-point of the part of the slice that is above the water surface. As for hydrodynamic force, the calculation will use the still water level or the instantaneous surface level, dependent on the parameter IFOADD, see Body Data Specification in the SIMO User Guide.

2.17. Specified forces

A number of excitation forces may be specified. Each force is defined by:

  • type (constant, harmonic, ramp (constant derivative))

  • coordinates

  • magnitude

  • attack method (degree of freedom or attack point and direction. The direction may be global or follow the body.)

2.18. Small-body hydrodynamic forces

The small volume body is described as a point with mass and volume. The motion is described by up to 3 degrees of freedom, all translations.

The hydrodynamic forces on the object are found according to Morison’s formula

\[q_h=\rho V'c_m\dot {u}_r+u_r\big(c_L+|u_r|c_Q\big)\]

where:

\(V'\)

submerged volume of the body

\(c_m\)

added mass coefficient

\(u_r, \dot{u}_r\)

relative velocity and acceleration between body and water particles

\(c_L, c_Q\)

linear and quadratic drag at velocity = 1

According to (Lehn and Øritsland, 1987) and (Sarpkaya and Isaacson, 1981), the dynamic equilibrium equation can be written

\[q=\frac{\partial (m\dot {x})}{\partial t}=\frac{\partial m}{\partial t}\dot {x}+m\frac{\partial \dot {x}}{\partial t}\]
\[\frac{\partial m}{\partial t}\dot {x}=\frac{\partial m}{\partial z}\frac{\partial z}{\partial t}\dot {x}=\frac{\partial m_h}{\partial z}\frac{\partial z}{\partial t}\dot {x}\]
\[m\frac{\partial \dot {x}}{\partial t}=m\ddot {x}\]

where:

\(m\)

mass and hydrodynamic mass of the body

\(m_h\)

hydrodynamic mass of the body

\(z\)

depth of the object relative to instantaneous sea surface

\(\rho\)

density of water

The dynamic equilibrium equation is written

\[(m_0+m_{h,i})\ddot {x}_i=q_i+q_{w,i}+q_{s,i}\]

where:

\('\)

denotes position-dependent parameters/coefficients

\(i\)

denotes degree of freedom (\(i=1,2,3\))

\(x,\dot{x},\ddot{x}\)

body position, velocity and acceleration

\(m_0\)

(dry) mass of the body

\(m_h\)

hydrodynamic mass of the body

\(q_i\)

hydrodynamic force in addition to \(m_{h,i} \ddot{x}_i\)

\(q_w\)

partially submerged body weight, \(q_w = m_0g - \rho V'g = m_0g - \rho V_0 c_vg\)

\(\rho\)

density of water

\(q_s\)

other static or time-dependent forces such as coupling forces

\(V_0\)

fully-submerged volume

\(g\)

acceleration of gravity

The hydrodynamic forces, \(\mathrm {q_i}\), are found according to the following equation

\[q_i=\big(\rho V'+m_{h,i}\big)\dot {u}_i+\big(c_L+|u_r|\:c_Q\big)u_{r,i}+\frac{\partial m_{h,i}}{\partial z'}\big(\dot {x}_3-\dot {x}_{3,W}\big)u_{r,i}\]

where:

\(c_L, c_Q\)

linear and quadratic drag for unit velocity

\(u_{r,i}\)

relative velocity between body and water particles in \(i\)’th direction

\(\dot{u}_i\)

wave particle acceleration in \(i\)’th direction

\(u_r\)

relative velocity, \(\displaystyle u_r = \sqrt{\sum_{i=1}^3 u_{r,i}^2}\)

\(z'\)

vertical displacement between object and sea surface

\(\dot{x}_{3,W}\)

velocity of particle on sea surface in vertical direction

2.19. Soil penetration forces

The soil penetration models apply for both skirted and piled structures. The soil penetration model is related to a small-volume body, which means that all rotations are excluded for the model.

There are three conditions which are considered in the modelling of soil resistance to seabed penetration. These are:

  1. The structure causes no pressure build-up such that the soil resistance against vertical movement is from vertical friction forces and resistance on the skirt/pile tip area.

  2. The structure has cells that are enclosed, with small ventilation channels. These channels allow water flow in and out of the cells. Suction may be applied after a predefined time. The net pressure force inside the cells adds to other accelerating forces.

  3. As 1, but the body is allowed to slide horizontally with a friction force in the horizontal direction, depending on the penetrating depth.

For skirted structures, soil resistance to penetration into the seabed is dependent on:

  • The degree of cyclic motion of the structure.

  • How much of the total skirt area is directly in contact with the sediments.

  • Whether the cells within the skirts are completely enclosed or not in soil around the structure’s perimeter.

For model 1 and 2 the penetrating body is arrested in the horizontal direction at landing. This is accomplished by setting the horizontal velocities, and the accelerating forces in horizontal direction, equal to zero. This implies a simplification, no sliding along the sea bed is included.

The body motion will be finally stopped when either:

  • the soil friction equals the net downward force from the structure, or

  • the body reaches its depth of full penetration

The soil penetration models should be combined with depth-dependent hydrodynamic coefficients to account for any proximity effects from the sea bed on added mass and damping properties.

The modelling of soil reaction forces includes:

  • Soil buoyancy force (if the soil can be regarded as a fluid)

  • Soil friction and bearing force

For structures with closed compartments, the modelling also includes:

  • Pressure forces in undisturbed soil, including suction.

  • Pressure buildup in a fractured compartment caused by hydrodynamic forces related to the underside of the closed compartments. These forces may give additional fractures.

2.19.1. Soil buoyancy force

The buoyancy due to higher density of the soil compared to water is included in the model as

\[F_b=(\rho _s-\rho _w)A_shg\]

where \(\mathrm {\rho _s}\) is mass density of soil, \(\mathrm {\rho _w}\) is density of water, \(\mathrm {A_s}\) is cross-section area of penetrating part of the structure, and \(\mathrm {h}\) is penetration.

2.19.2. Pressure buildup in a closed compartment

The volume flow is calculated as

\[q=-c\sqrt{p_c}\]

Where \(\mathrm {c}\) is the flow coefficient and \(\mathrm {p_c}\) is the net pressure inside the cavity, including any suction. \(\mathrm {V}\) is the total volume flow since the cavity became sealed and is calculated as

\[V=\int_{t_s}^tq\mathrm {d}t\]

where \(\mathrm {t_s}\) is the time since sealing occurred.

The compressibility of water (and soil) is considered. Once the skirt reaches the bottom of a new failure (or the sea bed, if no failure occurs at landing) the closed compartment is sealed, and the pressure build-up starts. The static force due to pressure inside the cavity can be approximated as

\[F=k(\frac{V}{A}-\delta )\]

where:

cavity

\(\delta\)

penetration depth, down from the sealing level

\(k\)

effective stiffness due to compressibility of water and soil inside the cavity

\(A\)

section area of cavity

The stiffness, \(\mathrm {k}\), can be approximated by

\[k=\frac{EA}{L}\]

where:

taken as the depth of full penetration)

\(E\)

elasticity module of water \((=2150~\mathrm{MPa})\)

\(L\)

height of cavity (constant value used during the entire penetration, taken as the depth of full penetration)

When the skirt tip reaches the level of the undisturbed soil, the structure is arrested until the pressure inside the skirt exceeds the soil capacity against failure. This capacity depends on penetration depth.

The water entry and exit of the closed compartment is a dynamic model which sets requirements to the time step. Subdivisions of time steps should possibly be used. The soil penetration model should only be used with the modified Euler integration method.

2.19.3. Pressure buildup in a fractured compartment

When the soil has fractured, the soil force on the skirt will be side friction against motion in both directions and (in case of a concrete structure) also skirt tip pressure, calculated as for the open compartment. When the skirt reaches the bottom of this new fracture, \(\mathrm {h_f}\), the skirt tip is in undisturbed soil and the model above applies.

During this phase, the pressure inside the cavity is calculated from the hydrodynamic forces related to the underside of the closed compartment, and new fractures may occur. These hydrodynamic forces are calculated as:

Drag:

\[\big(B_2(z)-B_{20}\big)\:|\dot {z}|\dot {z}+\big(B_1(z)-B_{10}\big)\dot {x}\]

Inertia:

\[\rho V\big(C_a(z)-C_{a0}\big)\ddot {z}\]

Where \(\mathrm {B_1(z)}\), \(\mathrm {B_2(z)}\) and \(\mathrm {C_a(z)}\) are depth-dependent hydrodynamic coefficients and \(\mathrm {B_{10}}\), \(\mathrm {B_{20}}\) and \(\mathrm {C_{a0}}\) are values valid for no sea bed interference.

Each single fracture will open up a ventilation channel along a part of the skirt periphery. The width and depth of the channels at fracture must be assessed by the user, and the input parameters (soil forces and hydrodynamic forces) given accordingly.

Depending on the above penetration forces and the hydrodynamic parameters valid for ventilation underneath the skirt, the damping of motion can increase sharply with depth.

2.19.4. Soil friction force

The soil reaction force is calculated as:

Downward motion:

\[F_d=F_{d0}(z)+F_1(z)\]

Upward motion:

\[F_u=F_{u0}(z)-F_1(z)\]

where

resistance, depth-dependent

), depth-dependent

skirt wall or pile friction

\(F_{d0}\)

Bearing force for DOWNward motion, typically representing skirt tip resistance, depth-dependent

\(F_{u0}\)

Bearing force for UPward motion, usually taken as a fraction of \(F_{d0}\), depth-dependent

\(F_1\)

Friction force for both DOWNward and UPward, typically representing skirt wall or pile friction

This force is directed opposite to the velocity of the penetrating structure, or if the velocity is zero, opposite to the net force (in this case not exceeding the net force).

3. Hydrodynamic interaction between bodies

This option enables the user to include hydrodynamic interaction effects between an arbitrary number of bodies. These effects may be significant in case the distance between the floating bodies is small.

The hydrodynamic interaction effects affect the first and second order wave excitation forces and the added mass- and damping-forces. The effects on the wave excitation forces are included as part of the body forces, see Body forces, while the effects on frequency-dependent added mass- and damping-forces are included in coupled retardation functions and coupled added mass at infinite frequency.

By introducing the coupled added mass at infinite frequency into the equation of motion, Equation (12), the inertia term is written

\[\begin{bmatrix}(\boldsymbol{m}+\boldsymbol{A}_\infty)_{i,i}&(\boldsymbol{A}_\infty)_{i,\,j}\\\\(\boldsymbol{A}_\infty)_{j,i}&(\boldsymbol{m}+\boldsymbol{A}_\infty)_{j,\,j}\end{bmatrix}\begin{bmatrix}\boldsymbol{\ddot {x}}_i\\\\\boldsymbol{\ddot {x}}_j\end{bmatrix}\]

where the indices \(\mathrm {i}\) and \(\mathrm {j}\) refer to body \(\mathrm {i}\) and body \(\mathrm {j}\). \(\mathrm {\boldsymbol{\ddot {x}}_i}\) and \(\mathrm {\boldsymbol{\ddot {x}}_j}\) are accelerations in 6 DOF referring to the respective body-fixed coordinate systems.

Symmetric mass properties are assumed, that is

\[(\boldsymbol{A}_\infty)_{j,i}=(\boldsymbol{A}_\infty)_{i,j}^T\]

The term associated with the retardation functions is written

\[\int_0^\tau\begin{bmatrix}\boldsymbol{h}+(t-\tau)_{i,i}&&\boldsymbol{h}+(t-\tau)_{i,\,j}\\\boldsymbol{h}+(t-\tau)_{j,i}&&\boldsymbol{h}+(t-\tau)_{j,\,j}\end{bmatrix}\begin{bmatrix}\boldsymbol{\dot {x}}_i\\\boldsymbol{\dot {x}}_j\end{bmatrix}\mathrm {d}\tau\]

The implemented formulation makes use of symmetric properties

\[\boldsymbol{h}(t-\tau)_{j,i}=\boldsymbol{h}(t-\tau)_{i,j}^T\]

4. Station-keeping forces

4.1. Anchor lines

The implementation of a catenary mooring line model in SIMO is based on the model used in the mooring analysis program MIMOSA. This includes both quasistatic analysis and a simplified dynamic analysis accounting for the effect of drag loading on the line. The line dynamics model in MIMOSA, however, operates in the frequency domain and is thereby linearized. In SIMO, the model is extended to time domain. The mooring lines are treated individually. They may therefore have completely different properties, like material, dimension, length etc. The mooring lines are composed of any number of segments, with the possibility of attaching buoys or clumpweights between them.

The mooring lines are assumed to form catenaries, modelled by the catenary equations.

The procedure for calculating the mooring line configuration is based on a "shooting method", or iteration on boundary conditions at one end in order to satisfy specified boundary conditions at the other. Using this procedure, a fairly accurate static equilibrium configuration for a multi segment line can be obtained with a minimum of computational effort.

The problem is solved in two dimensions. The line is assumed to remain in a vertical plane containing both the anchor, the touch down point and the fairlead. This means that the effect of transverse bottom friction is neglected so the line will rotate freely about the anchor.

The procedure is described in details in (Lie 1990).

For a pure quasistatic analysis the effects of transverse drag forces on the line are neglected. The total line tension and the angle of upper end from vertical are thereby determined by the location of the upper end relative to the anchor (both lying in the same vertical plane). A two dimensional line characteristics table is calculated.

At any later stage the line forces on a body may be found by the following procedure:

  1. The fairlead coordinate is transformed into the line plane.

  2. The total tension and direction are found by interpolation from the line characteristics table.

  3. The resulting force is transformed back to the body coordinate system.

Whenever line characteristics data are changed, for example due to winch running, the line characteristics table will be recalculated.

A Finite Element Model (FEM) is included as an option for tricky systems where the shooting method fails to converge.

Due to dynamic effects (i.e. effects due to velocity and acceleration of the line), the line tension may deviate from the tension predicted by quasistatic analysis. Transverse drag forces in particular will tend to increase the tensions relative to the quasistatic values.

A model for accounting for this effect has been developed by (Larsen and Sandvik, 1990). The model is based on four important assumptions.

  1. Only the tangential component of the top end motion is assumed to have any effect on the dynamic tension.

  2. The shape of the dynamic motion due to a tangential excitation is assumed to be equal to the change in static line geometry.

  3. Mass forces on the line are neglected.

  4. The elastic elongation of the line is determined quasistatically.

image400
Figure 4. Dynamic tension computation model.

The model is illustrated in Dynamic tension computation model. Moment equilibrium with respect to point \(\mathrm {P}\) yields

\[\int_0^S|\mathrm {d}\overrightarrow{f_c}\times \overrightarrow{r}|\mathrm {d}l=aT_{DC}\]

\(\mathrm {T_{DC}}\) is the tension caused by the drag loading on the line, and is given by

\[T_{DC}=k^E(x-u)-k^Gu\]

Manipulation with Equation (124) and Equation (125) leads to the following equation describing the line motion, \(\mathrm {u}\)

\[c^*\dot {u}|\dot {u}|+k^*u=k^Ex\]

where

\(c^*\) generalized line damping

\(k^*\)

generalized line stiffness

\(k^E\)

axial line stiffness, assumed to be constant

\(k^G\)

geometric catenary stiffness

\(u\)

displacement of the line

\(\dot{u}\)

velocity of the line

\(x\)

tangential motion excitation of the line end

The model may be understood by considering [Equivalent_geometric_model]. The coordinate \(\mathrm {x}\) expresses the imposed tangential motion of the upper line end, while the coordinate \(\mathrm {u}\) expresses the part of this motion that is compensated by geometric change of the line. The difference between \(\mathrm {x}\) and \(\mathrm {u}\) will thereby be the axial stretching of the line.

image405
Figure 5. Equivalent geometric model.

The coefficients for generalized line damping and stiffness \(\mathrm {c^*}\) and \(\mathrm {k^*}\) in Equation (141) are dependent on the position of the upper end point. As for the quasistatic case, a line characteristics table must be calculated. The characteristics table will comprise secant values for generalized damping and stiffness due to an excitation in the line tangent direction, as well as the tangent angle from vertical for each end point position.

Due to this position dependency of the coefficients, it is desirable to express Equation (126) on incremental form. The equilibrium of the line at time \(\mathrm {t_i}\) can be expressed

\[c_i^*\dot {u}_i|\dot {u}_i|+k_i^*u_i=k^Ex_i\]

After the time step \(\mathrm {\Delta t}\), the dynamic equilibrium will be

\[c_{i+1}^*\dot {u}_{i+1}|\dot {u}_{i+1}|+k_{i+1}^*u_{i+1}=k^Ex_{i+1}\]

By assuming constant coefficients within a time step and subtracting Equation (142) from

\[\Delta F_i^D+\Delta F_i^S=\Delta Q_i\]

Equation (144) expresses that the increment in drag and stiffness forces should balance the increment in loading due to top excitation. The different terms appear as follows

\[\Delta F_i^D=F_{i+1}^D-F_i^D=\overline{c^*}|\dot {u}_{i+1}|\dot {u}_{i+1}-|\dot {u}_i|\dot {u}_i\]
\[\Delta F_i^S=F_{i+1}^S-F_i^S=\overline{k^*}(u_{i+1}-u_i)=\overline{k^*}\Delta u_i\]
\[\Delta Q_i=Q_{i+1}-Q_i=k^E(x_{i+1}-x_i)=k^E\Delta x_i\]

where \(\mathrm {\displaystyle \overline{c^*}=\frac{1}{2}(c_{i+1}^*+c_i^*)}\) and \(\mathrm {\displaystyle \overline{k^*}=\frac{1}{2}(k_{i+1}^*+k_i^*)}\).

The equations are solved by an iterative procedure.

4.2. Thrusters

4.2.1. Thruster types

Thrusters are used for station-keeping, alone or in combination with mooring lines. Thrusters can be divided in two main categories, fixed-direction thrusters and rotatable thrusters, or azimuth thrusters. An azimuth thruster can usually rotate \(\mathrm {360^\circ}\) about the vessel’s vertical axis. In SIMO a fixed direction thruster is allowed to give a negative force, while an azimuth thruster always gives a force which is non-negative.

Further, a thruster can be open or ducted (nozzled). Although this distinction is reflected by the type designations for thrusters in SIMO, see Body Data Specification in the SIMO User Guide, it has no effect on the computation of thrust: A thruster’s behaviour is characterized by its coefficients of thrust and torque, which implicitly reflect whether it is open or ducted. In other words, the open/ducted property is not used internally inside SIMO (this may change in later program versions).

4.2.2. Coefficients of thrust and torque

A thruster’s thrust, \(\mathrm {T}\), and torque, \(\mathrm {Q}\), are given by

\[\begin{array}{l}T=n^2\rho D^4K_T(J)\\\\Q=n^2\rho D^5K_Q(J)\\\\\displaystyle J=\frac{v_a}{nD}\end{array}\]

Here \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) are dimensionless coefficients of thrust and torque, respectively. \(\mathrm {n}\) is the speed of rotation of the thruster shaft in revolutions per second, \(\mathrm {D}\) the screw diameter and \(\mathrm {\rho }\) the density of water. The thrust and torque coefficients are functions of the non-dimensional advance ratio, \(\mathrm {J}\). \(\mathrm {v_a}\) is the axial component of the relative velocity \(\mathrm {(\boldsymbol{v}_r)}\) of the water at the centre of the screw.

\[\begin{array}{l}\boldsymbol{v}_r=\boldsymbol{v}_{\mathrm {cur}}+\boldsymbol{v}_{\mathrm {orb}}-\boldsymbol{v}_T\\\\v_a=\boldsymbol{e}\cdot \boldsymbol{v}_r\end{array}\]

Here the vector \(\boldsymbol{v}_{\mathrm {cur}}\) is the current velocity, \(\boldsymbol{v}_{\mathrm {orb}}\) the orbital wave velocity and \(\mathrm {\boldsymbol{v}_T}\) the thruster velocity. \(\mathrm {\boldsymbol{e}}\) is the unit vector along the thruster axis, and the dot denotes scalar multiplication. Hence, \(\mathrm {v_a}\) is the projection of the vector velocity \(\mathrm {\boldsymbol{v}_r}\) on the thruster axis.

Figure 6 shows an example of thrust and torque coefficients for various values of the non-dimensional ratio \(\mathrm {P/D}\), where \(\mathrm {P}\) is the thruster’s pitch.

4.2.3. Reverse operation

A thruster which is intended primarily for rotation in one direction will have asymmetric blades that are designed favourably for that direction of rotation. Examples of this are an azimuth thruster and a ship’s main propeller. Reversing the direction of rotation of such a unit will produce less force at a given rotation speed.

Thrusters that are designed for two directions of rotation have symmetrical (e.g. elliptic) blade profiles. There may still be some asymmetry in thrust and torque owing to the existence of a transmission/gear box on one side of the thruster.

The asymmetry can be modelled by specifying separate coefficients of thrust and torque for reverse thruster operation. Alternatively and simpler, the coefficients for the forward direction can be used for reverse thrust after modification by a factor. This is offered as an alternative in the present model.

4.2.4. Data for thrust and torque coefficients

As there are four combinations of the signs of \(\mathrm {n}\) and \(\mathrm {v_a}\) in Equation (148) the thruster model is "4-quadrant". As pointed out above the coefficients for positive and negative \(\mathrm {n}\) will in general be different.

Hence two \(\mathrm {K_T}\) functions are required when the sign of the rotation speed can change, typically when the thruster is not of azimuthing type. The coefficients of thrust and torque can be specified in three ways:

  1. Separate \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) are given by the user for positive and negative rotation speeds \(\mathrm {n}\).

  2. \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) are given by the user for positive \(\mathrm {n}\) only, and user-chosen modification factors are applied for reverse rotation.

  3. Internally stored functions for \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) are used based on the \(\mathrm {P/D}\) ratio, which must be given by the user. The internal coefficients are given for positive \(\mathrm {n}\) only, and user-supplied modification factors are used for reverse operation.

In case 1 and 2 above, to specify \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) as functions of the advance ratio \(\mathrm {J}\), the user must enter a number of triplets \(\mathrm {\big\{J,~K_T(J),~K_Q(J)\big\}}\). \(\mathrm {J}\) can be positive or negative. A negative value of \(\mathrm {J}\) means that \(\mathrm {v_a}\) and \(\mathrm {n}\) have opposite signs. During simulation, SIMO uses linear interpolation between the given coefficient values. Outside the range of given \(\mathrm {J}\)’s, SIMO uses linear extrapolation.

The internal data (case 3) for \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) represent a nozzled thruster with four blades and an effective blade area ratio of \(\mathrm {0.7}\), see Figure 6.

As a simplification, \(\mathrm {K_T}\) and \(\mathrm {K_Q}\) need not be given separately for negative speeds of revolution \(\mathrm {n}\). Instead, reduction (modification) factors can be applied to the force when \(\mathrm {n<0}\). (For modelling freedom, modification factors are available for both the forward and the reverse direction of screw rotation.)

For tunnel thrusters the \(\mathrm {J}\)-dependence of the thrust coefficient is assumed not to be relevant. Hence, only \(\mathrm {K_T(0)}\) and \(\mathrm {K_Q(0)}\) are used. In general, different coefficient values are needed for the two directions of rotation.

image424
Figure 6. Built-in coefficients of thrust and torque for a number of pitch/diameter ratios. (Note that torque coefficients are multiplied by a factor of 10)

4.2.5. Thruster control

Speed control

Each thruster is controlled by a servo control system. The system consists of a variable frequency drive (VFD) and an electrical induction motor, see Figure 7 The task of the control system is to run the thruster at a given speed of rotation regardless of the variations in the thruster load. The input to the servo is the desired rotation speed \(\mathrm {n_0}\). The actual speed \(\mathrm {n}\) is measured and compared to \(\mathrm {n_0}\). If there is a difference, the power circuitry of the VFD will change the voltage and frequency of the current fed to the motor to make it apply a greater or lesser torque to the thruster.

When the thruster is a part of a DP system, \(\mathrm {n_0}\) is the rotation speed the DP controller wants the thruster to have, and which is sent as a command to the thruster servo. For conversion between thrust and rotation speed, the DP controller uses the "bollard pull" characteristics of the thruster.

image425
Figure 7. Feedback control system for thruster speed of rotation.

In Figure 7 \(\mathrm {n_0}\) is the desired rotation speed, \(\mathrm {n}\) is the actual speed (in revolutions per second). The variable frequency drive (VFD) consists of a control processor and power electronics. It compares \(\mathrm {n_0}\) and \(\mathrm {n}\) and feeds the motor with an electrical current of suitable voltage and frequency to make the motor apply the torque \(\mathrm {Q_M}\) necessary for keeping \(\mathrm {n}\) equal to \(\mathrm {n_0}\), or correcting \(\mathrm {n}\) if different from \(\mathrm {n_0}\). \(\mathrm {Q}\) is the load on the thruster.

For a given motor moment, \(\mathrm {Q_M}\), and load moment, \(\mathrm {Q}\), the angular acceleration of the thruster can be expressed by

\[2\pi I\dot {n}=Q_M-Q\]

Here \(\mathrm {I}\) is the effective moment of inertia of the thruster motor, thruster and transmission. \(\mathrm {Q}\) is expressed by Equation (148). The dot above \(\mathrm {n}\) denotes differentiation with respect to time. The factor \(\mathrm {2\pi }\) occurs since \(\mathrm {n}\) is in revolutions per second and not in radians per second. In steady state, \(\mathrm {Q_M=Q}\).

The servo controller calculates the motor torque as

\[Q_M=K_P(n_0-n)+\rho D^5K_Q(0)n_0^2\]

The controller has two terms. The last term is feedforward from the speed reference. The controller uses the thruster’s torque coefficient at \(\mathrm {J=0}\) to anticipate the load on the thruster. The first term is proportional control. The gain factor \(\mathrm {K_P}\) is used to compensate for errors in the anticipated load. Such errors may be caused by the thruster not operating under "bollard conditions". For instance, should the thruster ventilate, the load would be greatly diminished and cause motor runaway. The proportional control will prevent this from happening.

With \(\mathrm {K_P}\) set to zero, the controller will keep the motor torque constant. The rotation speed of the thruster will then vary according to the variation of the load \(\mathrm {Q}\). With a very large value of \(\mathrm {K_P}\), the rotation speed will be kept close to \(\mathrm {n_0}\), regardless of the load (as long as the motor can deliver the required torque).

Pitch variation, i.e. controllable pitch thruster, is not modelled in SIMO.

In addition to the parameters in the equations above, the performance of the thruster depends on the maximum torque, \(Q_M^{\mathrm {max}}\), the thruster can deliver. To make the thruster model simple to use, the parameters \(Q_M^{\mathrm {max}}\), \(\mathrm {I}\) and \(\mathrm {K_P}\) are not given directly, but derived from the following user-chosen parameters:

(J=0) )

\(F_{\mathrm{max}}\)

the maximum force the thruster can give in bollard pull condition \((J=0)\)

\(t_{\mathrm{rise}}\)

the minimum rise time

\(T_{CF}\)

the time constant of the servo controller

The maximum torque from the thruster motor drive is determined by SIMO as

\[Q_M^{\mathrm {max}}=D\frac{K_Q(0)}{K_T(0)}F_\mathrm {max}\]

The rise time, \(t_{\mathrm {rise}}\), of thrust is the time it takes to increase thrust from \(\mathrm {10~\%}\) to \(\mathrm {90~\%}\) of maximum when \(Q_M^{\mathrm {max}}\) is applied, as shown in Figure 9. The effective moment of inertia in Equation (150) is determined from the rise time as

\[I=0.67t_\mathrm {rise}\sqrt{\rho D^5K_Q(0)Q_M^\mathrm {max}}\]

The time constant of the controller, \(\mathrm {T_{CF}}\), is related to how quickly the controller reacts to a change in the load, \(\mathrm {Q_L}\). It is defined as

\[T_{CF}=\frac{2\pi I}{K_P}\]

The gain \(\mathrm {K_P}\) is inversely proportional to the time constant, i.e. a small value for \(\mathrm {T_{CF}}\) gives tight control - within the limit set by \(Q_M^\mathrm {max}\). The concept of time constant is illustrated in Figure 8, which shows how the ratio \(\mathrm {n/n_0}\) varies after a step change in the reference \(\mathrm {n_0}\).

image433
Figure 8. Illustration of time constant for a time constant value of 4 seconds

Figure 8 shows a step response. At \(t=T_C=4~\mathrm {s}\) the response value is \(\mathrm {63~\%}\) of the asymptotic value. The straight line has the initial slope of the response, \(\mathrm {=1/T_C}\).

image434
Figure 9. Quickest rising of rotation rate and thrust, beginning from zero.
Azimuth control

For rotatable thrusters the azimuthal direction is controlled by a servo controller. The model is simple and has two parameters: A time constant, \(\mathrm {T_{CD}}\), and a maximum value for azimuthal rotation speed, \(r_\mathrm {max}\).

\[\begin{array}{l}\displaystyle T_{CD}\frac{\mathrm {d}\psi}{\mathrm {d}t}=\psi_0-\psi\\\\\displaystyle |\frac{\mathrm {d}\psi}{\mathrm {d}t}|\leq r_\mathrm {max}\end{array}\]

Here, \(\mathrm {\psi}\) is the azimuthal angle and \(\mathrm {\psi_0}\) is the demanded angle from the DP controller. The model is of first order, which means that the azimuth change is "damping dominated", i.e. characterized by negligible inertia. This simplification is motivated solely by the wish to keep the number of parameters low.

4.2.6. Thrust reduction when thruster is close to surface

To model thrust reduction when thruster is close to the water surface, a reduction factor is used. In the extreme case air is drawn by the thruster, and the phenomenon is sometimes referred to as ventilation loss. The reduction factor is given as a function of the depth of the thruster’s centre relative to its diameter, as indicated by Figure 10 (note that in the figure a value of \(\mathrm {1.0}\) for the factor means no loss). No attempt is made to express the thrust reduction by a more intricate mathematical model: The user must supply the reduction factor as a table of points. The reduction factor is applied to both thrust and torque.

image437
Figure 10. Hypothetical example of surface proximity thrust reduction factor as a function of non-dimensional thruster submergence. Note that the reduction (loss) factor in the figure is defined such that a value of 1 represents no loss or reduction.

4.2.7. Direction-dependent loss factor for azimuth thrusters

When a thruster is installed on a vessel, its open-water characteristics may change due to various effects of interaction between the thruster and the hull. There may be friction losses when the thruster jet sweeps across the bottom of the hull, and there may be losses do to pressure if the jet hits a part of the hull, such as the pontoon of a semisubmersible platform. Further, the jet may be deflected due to interaction with the hull (Coanda effect).

These interaction effects may lead to a reduction in the net thruster force and a change in the direction of the force. For an azimuth thruster, these losses will in general depend on the orientation of the thruster, i.e. the losses will be direction-dependent. SIMO has no models for these losses, but offers a direction-dependent loss factor the user can specify as a table of numbers. The factor apply to the magnitude of the thrust, there is no correction for the directional change.

4.2.8. Resultant force and moment from the thruster system

The force from each thruster is modelled as a vector having the magnitude \(\mathrm {T}\) as given by Equation (148) and a direction equal to the geometrical direction of the thruster. For an azimuthing thruster the direction will vary, but the force vector will always exist in the vessel’s local xy-plane.

The moment vector from the thruster is the vector product of the thruster’s position vector in the vessel and the force vector. The total force and moment from the thruster system is obtained by simple vector summation.

4.3. Lift and drag force (rudder)

It is possible to model the effect of a rudder behind a propeller by associating a lifting surface unit (LSU) with a thruster. The force from the LSU will depend on the water velocity induced by the thruster, the current velocity and the velocity of the LSU.

The thruster-induced velocity (see (Lehn, 1992)) is calculated as

\[\begin{array}{l}\displaystyle U_0=\sqrt{\frac{F}{\rho A}}\\\\\displaystyle U_\mathrm {mean}=\frac{U_0}{(0.89+0.149h)\sqrt{2}}\end{array}\]

\(F\)

the thruster force

\(\rho\)

water density

\(A\)

area calculated from thruster diameter

\(U_0\)

velocity at duct outlet at centreline of propeller race

\(h\)

distance from thruster to LSU

\(U_\mathrm{mean}\)

mean velocity at location of LSU (race where \(U = 0.5 \: U_\mathrm{centreline}\))

The force coefficients of the LSU are given as quadratic force coefficients of drag, lift and moment as functions of rudder angle.

A LSU can be also used to model a general fin. In this case the user has to specify a dummy thruster giving zero force to associate the LSU with this thruster.

4.4. Force-elongation with/without hysteresis

image440
Figure 11. Force-displacement relationship.

Any force-elongation relationship may be specified as a station-keeping force. The curves for increasing and decreasing elongation may be different, corresponding to hysteresis. Damping may be specified as a force proportional to any exponent of the velocity of the attack point.

The purpose of this model is to describe axial forces and shear forces in a general way. This model may for example be used for a fender system or for a hawser or anchor line.

The axial force, \(\mathrm {F_a}\), is written

\[F_a=F_{as}(d)+C_a(d)|\dot {d}|^p\mathrm {sign}(\dot {d})\]

where:

\(F_{as}(d)\) is the axial stiffness force

\(C_a(d)\)

is the axial damping coefficient

This is a general model that may incorporate a wide range of restoring force and damping models.

In Figure 12 a hysteresis model is shown. This model is obtained when \(\mathrm {p=0}\).

In Figure 13 a linear damping model is shown. This model is obtained when \(\mathrm {p=1}\). A force-displacement relationship is obtained by \(\mathrm {\boldsymbol{C}\equiv\boldsymbol{0}}\).

Two versions of the contact force may be specified:

Fixed contact points

image443
Figure 12. Hysteresis model, p = 0.
image444
Figure 13. Linear damping model, p = 1.

In this case, the distance \(\mathrm {d}\) is measured along a line connecting a point on the body with a fixed point in the global coordinate system.

The direction of the force is updated according to the instantaneous positions.

4.5. Docking device

Docking devices may assist during the final precision manoeuvre before landing or connection of a structure to a fixed point on a body or on the seabed. These devices are arranged in pairs (guide funnel and guide posts or docking cone/cylinder and guide pin).

The orientation is vertical for all typical load-handling cases. However, in order to utilise the model also for i.e. horizontal coupling between bodies, the orientation can be defined by the direction cosines of the axis of the guiding cone/cylinder.

Model characteristics and limitations

The model specified in the following will have the following characteristics:

  1. The model can be applicable both to coupling between bodies and to positioning of a body relative to a fixed point.

  2. Rotation symmetric stiffness around the longitudinal axis is assumed. The coupling force between the interacting bodies (or between a body and a fixed point) acts normal to the internal surface of the cylinder (which can be conical).

  3. When the tip of the circular pin enters into the cylinder opening, the transverse relative motion between the two will be limited by the difference in radius at the respective longitudinal distance, \(\mathrm {\Delta Z}\), see Figure 14. Outside this limit, structural stiffness will give a restoring force directed towards the centre.

  4. It will be detected if the entry into the cylinder opening is successful or not. Until a successful entry is made, the model will give zero forces.

  5. The forces in the axial direction of the cylinder surface is zero if friction forces are neglected.

image445
Figure 14. Longitudinal distance.

The longitudinal distance, \(\mathrm {\Delta Z}\), is shown in Figure 14.

Calculation method

The calculation of guiding force is made as follows:

Calculate the position in the global coordinate system of the entering unit, \(\mathrm {\boldsymbol{X}_1}\), and the end of the guiding cylinder, \(\boldsymbol{X}_{\mathrm {cyl}}\).

Find the vector from the cylinder end to the end of the pile

\[\boldsymbol{V}=\begin{bmatrix}X_1(1)-X_{\mathrm {cyl}}(1)\\\\X_1(2)-X_{\mathrm {cyl}}(2)\\\\X_1(3)-X_{\mathrm {cyl}}(3)\end{bmatrix}\]

The unity vector in axial direction of the cylinder is \(\mathrm {\boldsymbol{U}}\).

The pile entry into the cylinder is when the scalar product between the two vectors \(\mathrm {\boldsymbol{V}}\) and \(\mathrm {\boldsymbol{U}}\) changes from negative to positive sign, and the distance from the cylinder centre is smaller than the defined maximum radius of entry. After a successful entry, the distance between the cylinder axis and \(\mathrm {\boldsymbol{X}_1}\) is found from:

\[\boldsymbol{D}=\boldsymbol{V}-\boldsymbol{U}\cdot (\boldsymbol{U}\cdot \boldsymbol{V})\]

The axial motion into the cylinder is found from

\[\Delta L=\boldsymbol{U}\cdot \boldsymbol{V}\]

The contact force is then calculated by interpolation between the specified force tables, as a function of \(\mathrm {\Delta Z}\) and \(\mathrm {D}\). The force components are found by scaling with \(\mathrm {D(i)/|D|}\).

When the guide pin enters into a conical section of the guiding cylinder, the radial force component is found as described above. The axial component of the force is found as

\[F_A=-F_R\cdot \frac{\partial R_0}{\partial Z}\]

where \(\mathrm {\displaystyle \frac{\partial R_0}{\partial Z}\:(<0)}\) is the rate of radius reduction with distance into the cone. At each \(\mathrm {\Delta Z}\), \(\mathrm {R_0}\) is the largest radius in the force characteristics giving zero force.

4.6. Fender

A fender is here defined as a contact element, either:

  • between a body and a globally fixed point ("positioning element"), or

  • between two bodies ("coupling element")

The implemented model will give:

  • a compressive force parallel to the normal vector of the sliding plane

  • a friction force along the sliding plane

The fender can either be point symmetric, giving friction in all directions along the sliding plane, or it can be defined as a roller. In the latter case only motion parallel to the rotation axis results in friction. Possible applications are sketched in Figure 15. Any fender type can be selected for the shown alternatives.

image452
Figure 15. Fender alternatives, examples.

4.6.1. Coordinates

A berthing ("positioning") fender can either be attached to a globally fixed point or to the moving body. The other contact point of the fender is then on the defined fender (friction) plane.

Similarly, a ("coupling") fender between two bodies is attached to one of the bodies. The fender plane will then be on the other body.

The attachment point is given in local body coordinates if the fender is attached to a body, else in global coordinates.

A point and a normal vector, both given in the relevant coordinate system define the fender plane. The orientation and the dimension of the fender plane in the local x- and y-direction also have to be specified. Table 1 summarises the alternatives.

The force from a fender will be zero if the projection of the fender point on the fender plane is outside the rectangle of the fender plane and also if the fender point enters the fender plane from "below".

Table 1. Fendering options
Fendering option Positioning Coupling

Fender attached to:

Fixed point

Body

Body 1

Attachment point

\(\mathrm {X_{FG}}\)

\(\mathrm {X_{FB}}\)

\(\mathrm {X_{FB1}}\)

Fender plane point

\(\mathrm {X_{SB}}\)

\(\mathrm {X_{SG}}\)

\(\mathrm {X_{SB2}}\)

Normal vector, fender plane

\(\mathrm {X_{NB}}\)

\(\mathrm {X_{NG}}\)

\(\mathrm {X_{NB2}}\)

4.6.2. Force Calculation

Normal force

The compressive normal force is found by interpolation, from a specified relation between distance and force and from the specified internal damping. The distance between the two points \(\mathrm {X_F}\) and \(\mathrm {X_S}\) projected into the normal vector of the sliding plane is used in this interpolation. In the following it is assumed that the coordinates are transformed to the coordinate system the fender is attached to, confer Table 4.2.

The distance is given as

\[\displaystyle \boldsymbol{r}=\begin{bmatrix}X_F(1)-X_S(1)\\X_F(2)-X_S(2)\\X_F(3)-X_S(3)\end{bmatrix}\]

The normal vector is given as

\[\displaystyle \boldsymbol{n}=\begin{bmatrix}X_N(1)\\X_N(2)\\X_N(3)\end{bmatrix}\cdot \frac{1}{\sqrt{X_N(1)^2+X_N(2)^2+X_N(3)^2}}\]

The projected distance is given as

\[R=\boldsymbol{r}\cdot \boldsymbol{n}\]

The distance vector is given as

\[\boldsymbol{R}=R\cdot \boldsymbol{n}\]

The updated contact point on the sliding plane, \(\mathrm {\boldsymbol{S}}\), neglecting shear deformation of the fender, is given as

\[\boldsymbol{S}=\boldsymbol{r}-\boldsymbol{R}\]

The compressive fender force is found from

\[\boldsymbol{F}_n=-(f(R)+c|\dot {R}|^e\cdot \frac{\dot {R}}{|\dot {R}|})\cdot \boldsymbol{n}\]

where:

\(f(R)\) fender characteristics

\(c\)

damping constant

\(\dot{R}\)

deformation velocity

\(e\)

specified exponent

Friction force

Point fender

As mentioned previously, the fender can either have a direction-independent friction coefficient, or it can be defined as a roller, which gives zero friction in the direction perpendicular to the rotation axis. The rotation axis must be parallel to the fender plane.

Shear stiffness of the fender is introduced, in order to make the model as realistic as possible, and to avoid numerical instability. The shear stiffness is assumed direction-independent, and thus the shear deformation will be in the direction of the friction force.

The friction force model is sketched in Figure 16.

image460
Figure 16. Friction force model.

Point symmetric fender

For a point symmetric fender, the friction force is given by the following, if the fender moves along the sliding plane.

\[\boldsymbol{F}_f=-\mu|\boldsymbol{F}_n|\frac{\boldsymbol{s}}{|\boldsymbol{s}|}\]

where \(\mathrm {\boldsymbol{s}}\) is the sliding motion along the plane (since the last time step).

The shear deformation of the fender is expressed by the shear stiffness, \(\mathrm {k_s}\), and the in-plane fender force is given by

\[\boldsymbol{F}=k_s\boldsymbol{s}_s\]

The fender will be stuck to the plane until the in-plane force exceeds the static friction.

\[\begin{array}{c}\displaystyle \boldsymbol{F}\leq -\mu_s|\boldsymbol{F}_n|\frac{\boldsymbol{s}_s}{|\boldsymbol{s}_s|}\\\\\mathrm {or}\\\\\displaystyle |\boldsymbol{s}_s|\leq \frac{\mu_s}{k_s}|\boldsymbol{F}_n|\end{array}\]

When the fender slides, the shear deformation will be

\[\boldsymbol{s}_s=\frac{\mu}{k_s}\boldsymbol{F}_f\]

The `undeformed' fender position is found by subtracting the shear deformation vector from the current position of the contact point. Movement of the contact point leads to sliding if the movement of the point is larger than the upper limit for the shear deformation, \(s_{s,\mathrm {max}}\). The shear deformation vector is updated every time step.

Roller fender

For a roller fender, the force is calculated from:

the unit vector defining the roller axis, \(\mathrm {\boldsymbol{v}}\), the sliding motion parallel to the roller axis,

\[\boldsymbol{s}_f=\boldsymbol{v}(\boldsymbol{s}\cdot \boldsymbol{v})\]

the friction force,

\[\boldsymbol{F}_f=\mu|\boldsymbol{F}_n|\frac{\boldsymbol{s}_f}{|\boldsymbol{s}_f|}\]

and the total fender force,

\[\boldsymbol{F}_F=\boldsymbol{F}_n+\boldsymbol{F}_f\]

The force at the other end of the fender will have the same magnitude and point of attack, and opposite direction, and will be transformed to its respective coordinate system. The point of attack is

\[\boldsymbol{X}=\boldsymbol{X\!F}'+\boldsymbol{R}'\]

where both vectors are transformed to the coordinate system of the "owner" of the attachment point.

If the fender does not slide on the plane, a higher friction coefficient, \(\mathrm {\mu_s}\), is used for the static friction, defining the sliding start.

Spherical fender

The spherical fender (as shown to the right in Fender alternatives, examples) is characterised by contact starting at a positive distance from the fender point.

The calculation of compressive force and friction force will be as for the point-shaped fender, with one exception: if the contact point moves, the fender either rolls on the friction plane or it slides. This problem is treated as follows:

  • Calculate the location of the present contact points (on the plane and on the fender) in the previous time step (global coordinates), assuming that they are both fixed to its respective body (or at a global position).

  • Calculate the motion vector from the previous time step to the present time step, for both points (zero if globally fixed).

  • Find the relative motion vector as the difference between the two motion vectors, and the component of this vector parallel to the fender plane.

4.6.3. Summary of the model and its restrictions

  • Both "positioning" and "coupling" fenders.

  • Only one sliding plane, without limitations, given by a point and its normal vector.

  • The sliding plane can be on either side of the fender.

  • Compression damping included.

  • Shear stiffness and deformation of the fender included.

  • Static friction not necessarily equal to sliding friction.

  • Roller fender, with zero friction along a defined direction, or fixed fender, with equal friction coefficient in all directions.

  • Spherical fender, where the compressive force starts at a specified distance from the fender center.

4.7. Bumpers

Application

The bumper element model is used to model contact force between a body and a globally fixed cursor or bumper, or contact forces between bodies. The model is particularly useful in the analysis of offshore installation operations, where deflectors / bumper bars are used to guide a module to its correct position and to protect existing equipment from impact damages.

The bumper element model can be applied for global positioning of a body, or as a contact coupling between bodies.

image471

Model

The bumper element model defines a pair of lines (i.e. bumper bars) with specified location and length. One line is always fixed to a body while the other line can either be fixed to another body or have a globally fixed location. The contact force between the lines (i.e. bars) is defined through a specified distance / force relation, which can be non-linear and include damping. Sliding friction between the bumper bars is not included.

Assume that the end points of the two lines have coordinates \(\mathrm {A_1}\) , \(\mathrm {A_2}\) and \(\mathrm {B_1}\), \(\mathrm {B_2}\), respectively, defined in the same coordinate system. The respective vectors along the two lines:

\[\begin{array}{l}\boldsymbol{a}=\{A_2\}-\{A_1\}\\\\\boldsymbol{b}=\{B_2\}-\{B_1\}\end{array}\]

A vector between the two lines (from \(\mathrm {B_1}\) to \(\mathrm {A_1}\)),

\[\boldsymbol{c}=\{A_1\}-\{B_1\}\]
image474

The distance between the two lines and the distance vector is found from

\[\begin{array}{l}d=\boldsymbol{c}\cdot \boldsymbol{e}_n\\\\\displaystyle \boldsymbol{e}_n=\frac{\boldsymbol{a}\times \boldsymbol{b}}{|\boldsymbol{a}\times \boldsymbol{b}|}\\\\\boldsymbol{d}=d\boldsymbol{e}_n\end{array}\]

If the distance between the lines, \(\mathrm {d}\), is smaller than the distance at contact, \(\mathrm {d_0}\), and if the contact point is located between the end points of both lines, the contact force between the elements, found by interpolation in the specified force characteristics table, will be parallel to \(\mathrm {\boldsymbol{d}}\).

The coordinates of the contact point, \(\mathrm {P}\), are found as

\[\{P\}=\{A_1\}-\boldsymbol{d}+p\boldsymbol{a}=\{B_1\}+r\boldsymbol{b}\]

which can be rewritten to

\[p\boldsymbol{a}-\boldsymbol{d}=r\boldsymbol{b}-\boldsymbol{c}\]

There is contact between the lines if \(\mathrm {0\leq p\leq 1}\) and \(\mathrm {0\leq r\leq 1}\).

In addition, the ends of all bumper elements are treated as if there are half-spheres at the ends. This means that there will be contact between the two bumpers shown in Bumper contact.

bumpercontact
Figure 17. Bumper contact.

5. Body components

The data group BODY COMPONENTS can optionally be specified together within the data defining each body, and includes all information about equipment at the attachment points of coupling elements that can be changed during an analysis. One feature is winch running.

5.1. Application

The BODY COMPONENT group contains a description of attachment points on a body:

  • Name and location (in local body coordinates).

  • Point characteristics (whether the end of a coupling wire is fixed at the attachment point or if the body slides along the wire span, cf. Guidewires).

  • Winch data (maximum acceleration and speed, wire lengths paid out and remaining on the drum, sequences of winch running).

Winch running can be applied to the following types of coupling elements:

Table 2. Coupling types for winch running

SIMPLE WIRE COUPLING

Wire cross section stiffness (EA) is given (Winch running also in previous versions)

MULTIPLE WIRE COUPLING

Several wire segments connected to a branch point (Winch running also in previous versions)

FIXED ELONGATION COUPLING

Spring with specified length - force relation (Winch running modifies reference length)

USER DEFINED COUPLING

For a positioning line (CATENARY SYSTEM DATA) the parameter XWINCH is run length of winch. This winch is a fictitious winch, and the input is unchanged. When the mooring line tension is calculated by interpolation on the line characteristics, (XOFF - XWINCH) is used as horizontal distance to anchor where XOFF is the physical horizontal distance between fairlead and anchor.

5.2. Guidewires

Brief description

An arbitrary number of guidewires can be defined and used for leading a module towards a foundation, or towards guideposts.

The guidewire model is applicable both on 3D bodies and 6D bodies. On 3D bodies each of the lines passes through one point/ring (modelling the guide funnel), and for a 6D body two points/rings can be defined for each guidewire, representing the end openings on the respective guide funnel. The ring(s) are specified by local position, diameter and orientation. A principal sketch is found in Guidewire alternatives. The ends of the guide funnel are given in local body coordinates, and these two points define the centre of the guide funnel.

The force from the guidewire on the body at each end of the guide funnel is shown in the same figure.

image480
Figure 18. Guidewire alternatives.

The model is also suitable for describing e.g. the pull-down of a buoyant module, using a wire through a sheave on the sea bed, see Pull-down.

image481
Figure 19. Pull-down.

5.3. Example

An example with use of coupling points and a guidewire is shown in System with 4 coupling points.

image482
Figure 20. System with 4 coupling points.

Both VESSEL_A and VESSEL_B are of body type 4 which means that they will not move during a simulation. The coupling point Point_on_VESSEL_A has a winch named WINCH1. This starts to run at \(t=10~\mathrm {s}\) and stops at \(t=160~\mathrm {s}\) with a velocity of \(0.8~\mathrm {m/s}\). The initial position of body LOAD is set to \(-58~\mathrm {m}\) and the static equilibrium position at \(\mathrm {t=0}\) is found to be \(-63.5~\mathrm {m}\). For a simulation with length \(200~\mathrm {s}\) the vertical position is plotted in Vertical motion of body LOAD.

image483
Figure 21. Vertical motion of body LOAD.

Because body LOAD slides along the coupling line without friction, the body does not get any motion in other degrees of freedom.

The SIMO "sys-file" for the system in System with 4 coupling points is:

'************************************************************
SYSTEM DESCRIPTION SIMO
'************************************************************
' txsys, 3 lines
Test of guide points
2006-03-01 9:03
May 2002
'LENUNI TIMUNI MASUNI GRAV  RHOW  RHOA
 m      s      Mg     9.81  1.025 0.00125
'DEPTH  DIRSLO   SLOPE
 1000.
'************************************************************
BODY DATA SPECIFICATION
'*************************************** V E S S E L _ A ****
'CHBDY
 VESSEL_A
Text line # 1
Text line # 2
Text line # 3
'IBDTYP
   4
'=============================
BODY LOCATION DATA
'=============================
'  Xglob    Yglob    Zglob    Phi    Theta     Psi
   0.0       0.0       0.     0.0     0.0      0.0
'====================================== BODY COMPONENTS ====
BODY COMPONENTS
'==============
COUPLING POINT
'--------------
 Attachment point for coupling line
 at body VESSEL_A
'Coupling point of type FIXED
'  chcopo           cpltyp   xcpl    ycpl    zcpl
  Point_on_VESSEL_A  FIXED    0.      0.      0.
COUPLING WINCH
'  chcowi    cwico
  WINCH1      NRUN
'  waccm  wvelm   druma  druml   nrun
   .1     1.      100.    200.    1
'  tstart   tstop   runvel
     10.     160.    .8
'***********************************************************
BODY DATA SPECIFICATION
'************************************** V E S S E L _ B ****
'CHBDY
 VESSEL_B
Text line # 1
Text line # 2
Text line # 3
'IBDTYP
   4
'=============================
BODY LOCATION DATA
'=============================
'  Xglob     Yglob    Zglob    Phi    Theta    Psi
   100.      0.0       0.      0.0     0.0     0.0
'************************************************************
BODY DATA SPECIFICATION
'************************************************************
'CHBDY
 LOAD
Text line # 1
Text line # 2
Text line # 3
'IBDTYP
   1
'=============================
BODY LOCATION DATA
'=============================
'  Xglob    Yglob    Zglob    Phi    Theta    Psi
   50.      0.0      -58.     0.0    0.0      0.0
'=============================
BODY MASS DATA
'=============================
'txmass, 2 lines
Example
2006-03-01 9:03
'xcog        ycog       zcog
 0.          0.0        0.
'-----------------------------
 MASS COEFFICIENTS
'-----------------------------
'   rm     rixx     riyx    riyy     rizx    rizy    rizz
  1000.    10000.   0.000   40000.   0.000   0.000   40000.
'-----------------------------
 GRAVITY FORCE INCLUDED
'-----------------------------
LINEAR DAMPING
txt1
txt2
'  dl1     dl2     dl3     dl4     dl5     dl6
   50.      0.      0.      0.      0.      0.
    0.     50.      0.      0.      0.      0.
    0.      0.     50.      0.      0.      0.
    0.      0.      0.    200.      0.      0.
    0.      0.      0.      0.    200.      0.
    0.      0.      0.      0.      0.    200.
' ====================================== BODY COMPONENTS ====
BODY COMPONENTS
' ==============
'2 coupling points of type GUIDE
COUPLING POINT
' --------------
First guide point
at body LOAD
' chcopo         cpltyp   xcpl   ycpl   zcpl
  Guide1_at_LOAD  guide    -14.   0.     10.
' dia   dv1   dv2   dv3
  .11     2.    0.    0.
' ==============
COUPLING POINT
' --------------
Second guide point
at body LOAD
' chcopo         cpltyp   xcpl   ycpl   zcpl
  Guide2_at_LOAD  guide    14.   0.     10.
' dia   dv1   dv2   dv3
  .11    2.    0.    0.
'**********************************************************
   COUPLING DATA
'**********************************************************
CPL1
' ===================
Coupling between
VESSEL_A and VESSEL_B
' ===================
SIMPLE WIRE COUPLING
' ===================

'  chbdy1    xbdy1    ybdy1    zbdy1
'  VESSEL_A     0.       0.       0.
 Point_on_VESSEL_A
' ------------------------------------ Load
'  chbdy2    xbdy2    ybdy2    zbdy2
  VESSEL_B     0.       0.       0.
' ------------------------------------
'  nguide
    2
 Guide1_at_LOAD
 Guide2_at_LOAD
'   EA    RLEN      FLEXC       DAMPSW
  8.8E4   147.       0.          1000.
'  ifmoco   ftime   btens
     0        0.    100.
'**********************************************************
END

5.4. Heave compensator

A heave compensator (HC) is used to minimize the motion of a suspended load, by paying out or taking in crane wire and thus compensating for the vertical crane motion.

The HC has been implemented as an attribute to the BODY COMPONENT / COUPLING POINT of the type FIXED.

To a coupling point equipped with HC, it is possible to connect coupling elements of the following types:

  • "Fixed elongation coupling"

  • "Simple wire coupling"

  • "Multiple wire coupling"

  • "Liftline coupling"

The numerical model of the HC is to be based on a method previously developed by Hans Berntsen, SINTEF. The intention is to facilitate modelling of an active HC, where the vertical motion of the suspended load is minimized, preferably so that the HC piston moves around its mean-stroke position.

It is not possible to combine the active HC with a tensioner, and is not possible to "de-couple" a coupling point to which a HC is attached.

The modelled active heave compensator is a cylinder type device, where oil is supplied in the pressure side of the piston, controlled based on measured vertical position of the crane top. A principal sketch is shown in Sketch of heave compensator system. In the present implementation a detailed model for internal function of the compensator is not included. The following features are disregarded:

  • Mass of the moving internal components (piston, rod, sheave, etc.)

  • Internal oil flow resistance

  • Friction between piston and cylinder and in sheaves

  • Total volume and elasticity of the oil

  • Large motion of the crane top due to crane manoeuvres or ballasting is not compensated.

image484
Figure 22. Sketch of heave compensator system
Table 3. Key parameters

\(\mathrm {Z_{cr}}\)

Instantaneous crane top displacement from initial static position \(\mathrm {[m]}\) .

\(\mathrm {Z_{cr0}}\)

Initial / static crane top position \(\mathrm {[m]}\)

\(\mathrm {Z_{p}}\)

Piston position \(\mathrm {[m]}\) . It is assumed that the piston moves symmetric ally about its half-stroke position, where \(\mathrm {Z_p=0}\)

\(\mathrm {N_{wc}}\)

Number of wire parts at compensator

\(\mathrm {N_{w}}\)

Number of wire parts between crane top and hook as modelled in SIMO

\(\mathrm {A_{c}}\)

Compensator piston area \(\mathrm {[m^2]}\)

\(\mathrm {S}\)

Cylinder stroke, allowable motion range of the piston \(\mathrm {[m]}\)

\(\mathrm {q}\)

Oil flow into / from cylinder \(\mathrm {[m^3/s]}\)

\(\mathrm {FACTOR}\)

Fraction of the crane motion that shall be compensated

\(A_{\mathrm {MAX}}\)

Clipping amplitude of piston motion

5.4.1. Controller

The controller uses two methods for heave compensation:

  • Feedback loop PD control

  • Forward-feed control with a Kalman-type filter

Feedback loop control of compensator piston motion

Difference between desired piston position and present position:

\[\Delta Z_{p_k}=Z_{cr}\frac{N_w}{N_{wc}}-Z_p\]

Differentiated and LP filtered:

\[\Delta \dot Z_{p_k}=\Delta \dot Z_{p_{k-1}}+\eta _f[\frac{(\Delta Z_{p_k}-\Delta Z_{p_{k-1}})}{dt}-\Delta \dot Z_{p_{k-1}}]\]

where:

\(\eta_f\)

\(=1-\exp{(-\displaystyle \frac{\mathrm{d} t}{T_f} )}\)

\(T_f\)

time constant in LP filter

\(dt\)

time step \(\mathrm{[s]}\)

Reference from feed-back control:

\[U_k=K_p(\Delta Z_{p_k}-T_d\Delta \dot Z_{p_k})\]

where:

\(K_p\) gain feedback loop

\(T_d\)

feedback derivative time

Feed-forward loop control, using estimated crane velocity and acceleration

Vertical crane velocity (if \(\mathrm {Z_{cr}}\) is not reduced by \(\mathrm {FACTOR}\) or clipped by \(\mathrm {A_{MAX}}\) ):

\[\dot Z_{cr_{k}}=\frac{Z_{cr_{k}}-Z_{cr_{k-1}}}{dt}\]

Vertical crane acceleration:

\[\frac{\ddot Z_{cr_{k}}=(Z_{cr_{k}}+Z_{cr_{k-2}}-2Z_{cr_{k-1}})}{dt^2}\]

Reference from feed-forward control:

\[U_{f_{k}}=U_{f_{k-1}}+\eta _{f}[K_f(\dot Z_{cr_{k}}+T_{fd}\ddot Z_{cr_{k}})-U_{f_{k-1}}]\]

where:

\(K_f\) gain in forward loop

\(T_{fd}\)

feed-forward derivative time

Control of oil flow

Total control reference velocity:

\[U_{t_{k}}=U_k+U_{f_{k}}\]

Capacity utilization found from low-pass filtered velocity:

\[\phi _{k}=\phi _{k-1}+\eta _{q}(U_{t_{k}}-\phi _{k-1})\]

Limitation: \(\mathrm {-1\le \phi _{k}\le 1}\)

\(\eta_q\) \(=1-\exp (-dt/tq)\)

\(T_q\)

time constant for hydraulic valve

Oil flow into the cylinder:

\[q_{k}=K_q(\phi _{k}+\phi _{k-1})\]

\(K_q\)

hydraulic valve characteristics

Piston position:

\[Z_{p_{k}}=Z_{p_{k-1}}+dt\frac{q_{k}}{A_c}\]

Limitation: \(\mathrm {|Z_{p_{k}}|\le S/2}\)

The length adjustment of the crane wire is the output from the compensator model, to SIMO:

\[\Delta L_{k}=Z_{p_{k}}\frac{N_w}{N_{wc}}\]

5.4.2. Selection of compensator parameters

Cylinder cross section area

If we neglect friction the mean compensator force must counteract the weight of the suspended load, \(\mathrm {W_{l}}\) , the hook, \(\mathrm {W_{h}}\) , and the wire below the crane top, \(\mathrm {W_{w}}\) :

\[F_c=\frac{N_{wc}}{N_w}(W_l+W_h+W_w)=\frac{N_{wc}}{N_w}\sum\,W=pA_c\]

A suitable piston area can be estimated by assuming a mean pressure \(p=150\,\mathrm {bar}=15000\,\mathrm {kN/m^2}\) .

Maximum oil flow and valve characteristics

\[q_\mathrm {max}=2K_q=A_c[\displaystyle \frac{(Z_{p_k}-Z_{p_{k-1}})}{dt}]_\mathrm {max}=\displaystyle A_c\frac{N_w}{N_{wc}}[\dot Z_{cr}]_\mathrm {max}\]
\[K_q=\displaystyle \frac{A_c}{2}\frac{N_w}{N_{wc}}[\dot Z_{cr}]_\mathrm {max}\]

\([\dot Z_{cr}]_{\mathrm {max}}\) = maximum compensated crane velocity

Feed-forward gain

\[K_f=\frac{A_c}{K_q}=\displaystyle 2\frac{N_{wc}}{N_w}\frac{1}{[\dot Z_{cr}]_\mathrm {max}}\]

Time constants for valve and forward feed loop

Recommended values: \(\mathrm {T_{fd}=T_q=0.8}\) seconds

Feed-back and LP filter parameters

These include:

\(K_p\) feedback gain

\(T_d\)

differentiation time constant

\(T_{fl}\)

low-pass filter time constant

It is recommended to adjust these parameters through series of test runs, as the performance will depend on the suspended load and the lifting system. Suggested start values: \(\mathrm {K_p=3,T_d=0,T_{fl}=0.2}\) .

A simulation time step longer than \(\mathrm {T_{fl}/5}\) will reduce the compensator performance, and is not recommended.

An example of input is given below:

BODY COMPONENT

' ============ Coupling point 2 at body TOP

COUPLING POINT

' ------------

Coupling point TOP\_CP2

at body TOP

' chcopo type x y z

TOP\_CP2 FIXED 0. 0. 0.

'-------------------------------------------------

MOTION COMPENSATOR

'cmoco

HC\_S

'itype

1

'limod factor ampmax ehla

1 1.0 10. 0

'nwire nwirec stroke acyl

2 2 4. .02

' hckp hctd hckf hctfd hckq hctq tf1

3.0 0.0 .75 0.8 0.0069 .8 0.2

5.5. Tensioner

For the couplings Simple wire coupling and Fixed elongation coupling, a TENSIONER may be given. The tensioner is a passive pneumatic hydraulic cylinder, where the supplied pressure holds the mean tension in the pipe, preferrably so that the piston moves around its mean-stroke position. The holding force at the mid-stroke position can be adjusted by adding or removing gas/oil at the pressure side.

The input may ge given as shown below:

' ===================================

BODY COMPONENTS

' ===================================

COUPLING POINT

'------------------------------------

Crane winch

Txt\_2

'chcopo cpltyp xcpl ycpl zcpl

COPO1 FIXED 0. 0. 0.

'------------------------------------

COUPLING WINCH

'------------------------------------

'chcowi cwico

Winch\_1 NRUN

' waccm wvelm druma druml nrun

0.5 1. 850. 1000. 0

'

TENSIONER

' ftens df\_dt stifft stroke IHLA

3200. 5. 80. 8. 1

' hla\_name

HLA\_TENSI

image489
Figure 23. The effect of a tensioner on a characteristics.

6. Coupling forces

For the couplings "Simple wire coupling", "Multiple wire coupling" and "Lift line coupling", the parameters FLEXC and DAMPSW may be given.

FLEXC (inverse of stiffness) is the flexibility in the wire attachment point. This flexibility should account both for flexibility in the construction and elasticity in the wire between the winch and the used wire attachment point.

DAMPSW is material damping in the line. The value can normally be set to \(\mathrm {1-2\%}\) of \(\mathrm {EA}\) , where \(\mathrm {E}\) is the modulus of elasticity and \(\mathrm {A}\) is the cross-sectional area.

6.1. Simple wire coupling

The simple wire coupling is modelled as a linear spring according to :

\[\Delta {l}={\frac{T}{k}}\]

where :

\(\Delta {l}\) elongation

\({T}\)

wire tension

\({k}\)

effective axial stiffness

The effective axial stiffness is given by :

\[\frac{1}{k}=\frac{1}{EA}+\frac{1}{k_0}\]

where :

\({E}\)

modulus of elasticity

\({A}\)

cross-section area

\({l}\)

unstretched wire length (may be variable with respect to time)

\(\frac{1}{{k_0}}\)

connection flexibility (crane flexibility)

\(\Delta {l}\)

change in elongation line

Knowing the position of each line end, the elongation and thereby the tension may be determined.

Material damping is included as:

\[F=\frac{C_w\Delta l}{l\Delta t}\]

6.2. Multiple wire coupling (sling module)

The possibility of coupling several wire segments in a common branch point is incuded. In this way wire systems including slings may be modelled.

All wire segments will have one end fastened in a body and the other in the common branch point. Axial stiffness properties of each segment, \(\mathrm {\boldsymbol{k_{i}}}\) , are found in the same way as for the single wire coupling. To determine the location of the branch point, an iteration procedure must be used.

    1. A location of the branch point, \(\mathrm {\boldsymbol{X_{b}}}\) , is assumed.

    1. Knowing the location of the body fixed segment ends, and the stiffness property of each segment, the tension in each segment and the resulting force on the branch point, \(\mathrm {\boldsymbol{F_{b}}}\) , may be calculated.

    1. If the branch point is in equlibrium \(\mathrm {(|\boldsymbol{F_{b}}|=0)}\) , goto (g).

    1. Establish the global stiffness matrix for the branch point, \(\mathrm {\boldsymbol{K}}\) .

    1. The corrective motion of the branch point is found by :

\[\Delta \boldsymbol{X_b}=\boldsymbol{K^{-1}F_b}\]
    1. Correct \(\mathrm {\boldsymbol{X_{b}}}\) and goto (b).

    1. Calculation finished.

6.3. Force-elongation characteristic

The model is identical to that described as station-keeping force. As a coupling force, any force-elongation relationship may be specified. The curves for increasing and decreasing elongation may be different in order to model hysteresis effects. Damping may be specified as a force proportional to any exponent of the relative velocity of the end points. This model is also used for passive motion compensators.

6.4. Lift Line Coupling

The lift line model gives a fast and quite accurate model for prediction of the total static and dynamic line end forces and load offset. This coupling type requires that the line is mainly vertical. The net weight of the body at the lower end of the line will be applied as a force on the line. It is therefore not possible to have two coupling lines down to the same body. The simplified lift line is modelled by the nominal Young’s modulus (\(\mathrm {E}\) ), diameter, transverse and longitudinal drag coefficient and the mass per meter length. The model gives zero compressive forces. The basic model is a straight line. A modification to the straight-line assumption is made to take into account the transfer of horizontal motion and forces from the crane to the lifted body due to WF motion, LF motion and transient motions (moving vessel).

Current forces and the damping forces of the line acting on the body and the load are calculated along a straight line. The line weight and buoyancy and the dynamic line mass forces are also calculated along this line.

The load offset is calculated using the cross-flow principle for the drag force components, \(\mathrm {q_{lw,d,i}}\) normal to the line:

\[q_{lw,d,x}=\frac{1}{2}\rho C_dA|\boldsymbol{v_r}|\nu_{rx}\]

and similarly in y- and z-direction. The density of water is denoted \(\mathrm {\rho }\), \(\mathrm {C_d}\) is the drag coefficient and \(\mathrm {\boldsymbol{v_r}=[\nu_{rx},\nu_{ry}\nu_{rz}]^T}\) is the relative velocity between the line and the water particles. The cross-flow principle is based on the normal flow velocity. The velocity is decomposed into a normal velocity and a tangential velocity along the lift line, hence vertical forces will be present for horizontal velocities on non-vertical lift lines even for relative small deflection angles from the vertical axis. Initially the lift line stretching is taken into account by a linear force elongation relationship assuming uniform strain along the lift line:

\[q_{lw,el,\mathrm {initial}}=\frac{EA}{L}\Delta l\]

where \(\mathrm {L}\) is the total line length, and \(\mathrm {\Delta l}\) is the elongation of the line.

A straight-line model overestimates the transfer of the horizontal motion from the surface vessel to the load. A modification to the straight-line assumption is made in order to take into account the transfer of the horizontal motion and forces from the crane to the lifted load due to the WF motion. A vibrating string analogy is used in order to estimate the length of the lift line that is influenced by the crane vessel WF motion, and to express that this length increases with increasing tension. By using the vibrating string analogy, the wavelength \(\mathrm {\lambda_1}\) of the first order motion mode is given by:

\[\lambda_1=T_1\sqrt{\frac{|q_{lw,\mathrm {top}}|}{m'}}\]

where \(q_{lw,\mathrm {top}}\) is the top tension, and \(\mathrm {m'}\) , is the mass per unit length of the lift line. The first order mode period \(\mathrm {T_1}\) is substituted by the zero-crossing period of the wave spectrum \(\mathrm {T_z}\) . Part of the line or the whole submerged line length is exposed to motions according to the vibrating string analogy. \(\mathrm {\lambda_1/4}\) is used as the basis for the calculation of the forces on the part of the submerged line length that is influenced by the WF crane vessel motion. The vessel-induced WF velocity of the line from the vessel decreases linearly down along the lift line. Hence, the drag force due to this motion is quadratically decreased due to the quadratic nature of the drag forces. These velocities and the other velocity components acting on the line are superposed. Alternatively, due to the vibrating string analogy it can be argued that the motions should have a periodical shape (a quarter of a cosine wavelength), but the linearly decreasing velocity is assumed here. The assumed shape of the vessel induced WF motion neglects any other velocity components acting on the line, a highly questionable but useful method.

This simplified model has been extended to include the weight and buoyancy of the line and also calculate the dynamic top and bottom tension due to line mass forces. Lower end of the line is considered as a free end. The top tension \(q_{lw,w,\mathrm {top}}\) due to the weight and buoyancy for a line fully submerged in water and with uniform weight distribution along the line, as shown in Line elongation and tension due to weight and buoyancy for a fully submerged line with uniform weight distribution., is:

\[q_{lw,w,\mathrm {top}}=(m'-\rho A)gL=w_w'L\]

where \(\mathrm {g}\) is the acceleration of gravity and \(\mathrm {w_w'}\) is the unit weight in water. The line elongation \(\mathrm {\Delta l_w}\) due to the weight and buoyancy is:

\[\Delta l_w=\int_{0}^{L}\!{\frac{w_w'z}{EA}}\textrm{d}{z}=\frac{1}{2}\frac{w_w'L^2}{EA}\]
Line elongation and tension due to weight and buoyancy for a fully submerged line with uniform weight distribution.

image512

Line elongation and tension due to weight and buoyancy for a partly submerged line with uniform weight distribution.

image513

In Line elongation and tension due to weight and buoyancy for a partly submerged line with uniform weight distribution., the line is partly submerged, and the top tension and the line elongation due to the weight and buoyancy is then:

\[\begin{array}{lll}q_{lw,w,\mathrm {top}}&=m'g(L-L_w)+w_w'L_w=m_a'(L-L_w)+w_w'L_w&(4.178)\\\\\Delta l_w&=\displaystyle \frac{1}{2}\frac{w_w'L_w^2}{EA}+\frac{w_w'L_w(L-L_w)}{EA}+\frac{1}{2}\frac{w_a'(L-L_w)^2}{EA}&(4.179)\end{array}\]

where \(\mathrm {w_a'}\) is the unit weight in air of the line and \(\mathrm {L_w}\) is the submerged line length.

Line elongation and tension due to linearly distributed end accelerations of the line mass.

image519

Considering linearly distributed acceleration along the line (local z-axis from the line end), as shown in Figure 24, the tension due to the accelerations of the line from \(\mathrm {\ddot z_C}\) (crane acceleration) and \(\mathrm {\ddot z_L}\) (load acceleration) is:

\[\begin{array}{lll}q_{lw,a}(z)&=\displaystyle m'(\ddot z_Lz+\frac{1}{2}(\ddot z_C-\ddot z_L)\frac{z^2}{L})&(4.180)\\\\q_{lw,a}(z=L)&=\displaystyle q_{lw,a,\mathrm {top}}=\frac{1}{2}m'(\ddot z_C+\ddot z_L)L&(4.181)\\\\q_{lw,a}(z=0)&=0&(4.182)\end{array}\]

The elongation \(\mathrm {\Delta l_a}\) due to the accelerations of the line is then:

\[\Delta l_a=\frac{m'}{EA}\int_{0}^{L}\!{(\ddot z_Lz+\frac{1}{2}(\ddot z_C-\ddot z_L)\frac{z^2}{L})}\textrm{d}{z}=\frac{m'L^2}{EA}(\frac{1}{6}\ddot z_C+\frac{1}{3}\ddot z_L)\]

6.5. Ratchet coupling

The rachet coupling can take tension or compressive force. The force will start to act after a specified time. The force consists of 3 components: a constant force, a stiffness force, and a damping force. The constant force component (which can be set to zero) acts from \(\mathrm {t=0}\) . For a ratchet with compressive force, the damping will only act when the distance is decreasing and for a "tension ratchet" the damping will only act when the distance is increasing.

An example with a tension ratchet is shown in Figure 24. From \(\mathrm {t=0}\) to point A, the force in the ratchet will be equal to the specified constant force. The distance in point A, will be used as a reference distance for the stiffness force until point C is reached. The distances in point C and D are new distance references. In the green part of the curve, the force is equal to the constant force component.

It is possible to specify an event of ratchet failure. It will then fail if a specified force is exceeded. When the element fails, the tension is set to zero.

image527
Figure 24. Example of the force in a tension ratchet.

If zero is specified both for the stiffness and the (constant) minimum force, the ratchet element can be used to model hydraulic/pneumatic cylinders with direction-dependent damping.

An example of a ratchet coupling is shown in Figure 25.

image528
Figure 25. Simulation example

Body A has a forced sinusoidal motion with amplitude \(1\mathrm {m}\) and period \(10\mathrm {s}\) .

Body B can slide without friction along the x-axis. Body B has a mass of \(1000\,\mathrm {tons}\) and a linear damping of \(10000\,\mathrm {kNs/m}\) .

The resulting motion and force is shown in Example of the motions and ratchet force in the simulation example.

Example of the motions and ratchet force in the simulation example

image529

6.6. Stiff couplings

Stiff couplings such as articulated joints and hinges may couple a body to the ground. Such elements may be modelled by extremely high stiffness.

7. Dynamic positioning

The dynamic positioning module is a control module with the following input and output:

Input:

  • Position measurement

  • Wind measurement

  • Anchor line force measurement

  • Thrust measurement

Output:

  • Desired resultant forces and moment from the thrusters

A PID controller and a Kalman filter-based controller exist. In addition it is possible to use an external dynamic positioning controller that can be coded by SITEF Ocean as a shared library (DLL), or use the internal SINTEF Ocean DP-system (MTDP) which is used for model tests in the SINTEF Ocean Basin.

7.1. PID controller

The controller converts position and velocity errors into a demand for thrust forces to correct the errors. A decoupling approach allows PID control parameters to be specified separately for surge, sway and yaw.

PID characteristics

The control law of the PID (Proportional + Integral + Derivative) controller is

\[\begin{array}{ll}F_{T0}(t)&=K_D\dot \varepsilon (t)+K_p\varepsilon (t)+K_I\displaystyle \int_{0}^{t}\!{\varepsilon (\tau)}\textrm{d}{\tau}\\\\\varepsilon (t)&=x_0-x(t)\\\\\dot \varepsilon (t)&=\dot x_0(t)-\dot x(t)\end{array}\]

where

\(x(t)\) filtered position

\(x_0(t)\)

desired (target) position

\(\varepsilon (t)\)

position error

\(F_{T0}\)

desired control force from thrusters

\(K_D\)

velocity feedback gain

\(K_P\)

position feedback gain

\(K_I\)

integral feedback gain

Figure 26 shows an example of the PID controller’s frequency response from the position error \(\mathrm {\varepsilon }\) to control force \(\mathrm {F_{T0}}\) . It is seen that for low frequencies (below \(\mathrm {\omega _{I}}\) ) the controller acts like an integrator. In the medium frequency range it has proportional action, and for high frequencies (above \(\mathrm {\omega _{D}}\) ) the controller behaves like a differentiator. With increasing frequency, the phase angle of the controller goes from \(-90\mathrm {deg}\) to \(+90\mathrm {deg}\) . The relations between the frequencies \(\mathrm {\omega _{I}}\) and \(\mathrm {\omega _{V}}\) and the feedback gains are

\[\begin{array}{ll}K_p&=K_P'(1+\displaystyle \frac{\omega _I}{\omega _D})\\\\K_I&=K_P'\omega _I\\\\K_D&=\displaystyle \frac{K_P'}{\omega _D}\end{array}\]

where \(\mathrm {K_P'}\) is the level of the horizontal asymptote in the figure. Usually \(\mathrm {\omega _I\ll\omega _D}\) , so that \(\mathrm {K_P'\approx K_P}\)

As \(\mathrm {\omega }\) approaches zero the controller gain approaches infinity due to the integral term. As a result, there will be no static (mean) position error - which is the purpose of integral action. Without the integral term (\(\mathrm {K_{I}=0}\) or equivalently, \(\mathrm {\omega _{I}=0}\) ) the static error would be \(x_{\mathrm {stat}}=F_{\mathrm {ext}}/K_{P}\) , where \(F_{\mathrm {ext}}\) is the static component of the external disturbing force.

The purpose of the derivative action (i.e. the velocity feedback) is to provide phase lead. As is seen from Figure 26, the derivative term causes the phase to approach \(+90\,\mathrm {deg}\) . In general, a certain amount of phase lead will be required in order to obtain stability.

Filtering

It appears from Figure 26 that the PID’s gain also approaches infinity for high frequencies. This will make the system very sensitive to high frequency noise in the position and velocity signals as this noise would be amplified and put through to the thrusters. For this reason, high frequency components are removed from the signals by filtering before they are fed to the controller.

There are two filters. One serves the double purpose of 1) low-pass filtering the position measurement to reduce the amount of measurement noise, and 2) differentiating the filtered position. The result is the input \(\mathrm {x}\) and and its derivative. Frequency response of low-pass filter (solid lines) and combined low-pass filter/differentiator (dashed) shows the frequency response of the low-pass-filter/differentiator, which has the cut-off frequency \(\mathrm {\omega _{C}}\) as the only parameter. The cut-off frequency is the frequency at which the attenuation starts. Above this frequency the attenuation is \(20\,\mathrm {dB}\) (10 times) per frequency decade.

image535
Figure 26. Frequency response of a PID controller (Dashed lines are asymptotes)
Frequency response of low-pass filter (solid lines) and combined low-pass filter/differentiator (dashed)

image536

image537
Figure 27. Frequency response of wave filter.

In Frequency response of low-pass filter (solid lines) and combined low-pass filter/differentiator (dashed), the frequency axis is scaled as \(\mathrm {\omega /\omega _{c}}\) where \(\mathrm {\omega _{c}}\) is the cut-off frequency.

The other filter is a wave filter. Its task is to remove the wave frequency components of the position measurement signal to prevent them from entering the thrusters. This is because thrusters usually have no chance of counteracting first-order wave forces. [tm-fig-freq-wave-filter] shows the frequency response of a wave filter for different values of the strength parameter \(\mathrm {\eta _{WF}}\).

The wave filter is of band stop type. It is characterised by a frequency parameter and a strength parameter. The frequency response of the filter is given in [tm-fig-freq-wave-filter] for a number of values of the strength parameter \(\mathrm {\eta _{WF}}\) , to be chosen in the range from \(\mathrm {0}\) to \(\mathrm {1}\) . The maximum attenuation of the wave filter is \(16.5\mathrm {dB}\) or about \(\mathrm {7}\) times. It is intended that optimum filter performance should be achieved simply by equating the filter frequency with the peak frequency of the wave spectrum. As a consequence of the typical asymmetry of wave spectra, the frequency parameter is chosen to the left of the filter’s centre frequency.

It is perfectly possible to filter off the wave frequency motion using the low-pass filter only. It is just a question of choosing the cut-off frequency low enough. The wave filter, however, is better suited for this task, since it gives more attenuation for less phase lag than the low-pass filter. The wave filter does not reduce any high frequency measurement noise. If such noise exists and is a problem, the cut-off frequency of the low-pass filter must be chosen low enough to give a sufficient suppression of it. In any circumstance the cut-off frequency should not be chosen any lower than necessary, in order to introduce as little negative phase as possible into control loop. On the other hand, the cut-off frequency should be chosen well below one half of the sampling rate of the measurements.

Mechanical analogy

In a DP system, the position and velocity feedback gains can be interpreted as stiffness and damping coefficients, respectively, and this mechanical analogy can be useful in determining these gains. However, due to the influence of the integral term and the filters, the analogy should not be pulled too far. Both the integral term and the filters may affect the system’s frequency response substantially, and the controlled system’s behaviour should therefore not be judged on the basis of \(\mathrm {K_P}\) and \(\mathrm {K_D}\) alone. A procedure for selecting the control parameters is included below.

Cross-coupling effects

The above description applies to each of the three controllers for surge, sway and yaw. When the vessel’s centre of mass is chosen as the point to be positioned and the mooring system (if any) is fairly symmetric, the effects of cross-coupling between the motions will in general be small. In this case it may be good enough to control the motions independently of each other. If, on the other hand, a point near the vessel’s bow is to be positioned, large coupling effects must be expected. In this case effective control may be difficult to achieve with independent controllers.

There are three kinds of cross-coupling effects: through mass forces, damping forces and mooring forces. SIMO uses the vessel’s mass matrix (dry mass + hydrodynamic mass) to "couple together" the independently calculated force outputs from the three PIDs. This happens automatically without participation from the user. In this manner compensation for the effect of mass coupling, which is well defined, is achieved. It is expected that also the damping crosscoupling is partially compensated for in this manner, so that approximate decoupled control of any point on a non-moored ship can be obtained. For a moored ship however, the characteristics of the mooring system must be taken into account when the PID parameters are selected.

Procedure for tuning the PID controller

Tuning a control system is no exact science. The following procedure should serve as a guideline.

  1. For each of the three motions surge, sway and yaw, determine position and velocity feedback gains \(\mathrm {K_{P}}\) and \(\mathrm {K_{D}}\) that give the wanted stiffness and damping. Set the integral gains \(\mathrm {K_{I}}\) to zero. Turn off the wave frequency filter (\(\mathrm {\eta _{WF}=0}\) ) and set the cutoff frequency of the lowpass filter to \(\mathrm {1/5}\) of the sampling rate.

  2. From an initial position of non-equilibrium and with no dynamic environmental forces acting, simulate the transient motion towards equilibrium (To avoid complications caused by possible cross-coupling effects, it is safest to do this for one motion at a time, simulating pure transients in surge, sway and yaw). The transients should be non-oscillating. If they are not, the damping coefficients must be reduced.

  3. Make a simulation with dynamic environmental excitation. Inspect the thruster outputs to see if the amount of wave frequency components is admissible. If it is, proceed to 5. If not, select the appropriate value for the strength parameter of the wave filter according to [tm-fig-freq-wave-filter].

  4. With the selected wave filter, simulate transients again. The wave filter introduces negative phase in the control loop, which reduces stability. If any transient is unstable or shows oscillatory behaviour that does not die out rapidly (i.e. during one or two half cycles), stability must be improved. If it is unacceptable to reduce the strength of the wave filter, the remedy is to: 1) make the controller phase more positive at the frequency of oscillation, and/or 2) reduce the controller gain at the same frequency. Remedy (1) is carried out by reducing \(\mathrm {K_{P}}\) while keeping \(\mathrm {K_{D}}\) constant. Try this choice first, for example by halving the value of \(\mathrm {K_{P}}\) . Remedy (2) is effectuated by reducing \(\mathrm {K_{P}}\) and \(\mathrm {K_{D}}\) by the same amount (e.g. halving them). Try this remedy if (1) fails or is insufficient.

  5. If integral action is desired, determine the upper frequency limit \(\mathrm {\omega _{I}}\) for this action, i.e. the frequency above which the integral action becomes insignificant (cf. [tm-fig-freq-wave-filter]) Then set \(\mathrm {K_{I}=w_{I}\cdot K_{P}}\) and simulate transients. The integral action also deteriorates stability. If stability is unsatisfactory, \(\mathrm {K_{I}}\) must be reduced.

  6. If it is known that the amount of measurement noise (the present version of SIMO does not simulate such noise) necessitates harder low-pass filtering of the position measurements, the cut-off frequency of the filter must be reduced from the value set according to paragraph 1 above. This may then lead to an adjustment of the controller gains as outlined in paragraph 4, until a satisfactory transient motion is obtained for all three motion variables.

The above procedure reflects the standard knowledge from control theory that for a feedback system to be stable, the open-loop gain must be less than \(1(0\,\mathrm {dB})\) at the crossover frequency (i.e. the frequency at which the open-loop phase is \(-180\,\mathrm {deg}\) .)

7.2. Kalman filter

The Kalman filter module was originally developed at SINTEF, division of Automatic Control during the period 1975-77, and has later been maintained and slightly modified by MARINTEK.

The system is based on an optimal state estimator which is running in parallel with the physical process (the SIMO vessel). With feedback from estimated positions and velocities, the estimator acts as an adaptive controller.

The mathematical model of the vessel is divided in three main sections:

  • Model of wave-induced motions, high frequency (HF) model.

  • Vessel motions in calm water, low frequency (LF) model.

  • Model for slowly varying forces from waves and current.

Slowly varying forces due to current and second order wave forces can be modelled by two methods:

  • Estimation of global current velocity components

  • Estimation of local bias forces

In both cases, other unknown forces will also be estimated. The bias esimation represents the integrator effect of the controller.

Input data to the Kalman filter routines are:

\(x_\mathrm{ref}^G\)

global x-reference

\(y_\mathrm{ref}^G\)

global y reference

\(\psi_\mathrm{ref}\)

heading reference

\(x_\mathrm{loc}^L\)

x-coordinate of point on vessel to be positioned

\(y_\mathrm{loc}^L\)

y-coordinate of point on vessel to be positioned

\(d_x\)

longitudinal quadratic drag coefficient

\(d_y\)

transverse quadratic drag coefficient

\(d_\psi\)

rotational quadratic drag coefficient

\(v_\mathrm{typ}\)

typical velocity

\(m_x\)

body mass and added mass in surge

\(m_y\)

body mass and added mass in sway

\(m_\psi\)

body inertia and added inertia in yaw

\(R(3, \, 4)\)

Kalman filter gain matrix

\(G(3, \, 2)\)

linear controller matrix

Each time step, the following measurements are input to the Kalman filter:

\(F_{Mx}^L\) measured longitudinal forces

\(F_{My}^L\)

measured transverse forces

\(F_{M\psi}\)

measured moment

\(x_M^G\)

measured global x-position of body origin

\(y_M^G\)

measured global y-position of body origin

\(\psi_M\)

measured heading

The output from the routines are:

\(F_{Tx}^L\) longitudinal thrust demand

\(F_{My}^L\)

transverse thrust demand

\(F_{T\psi}\)

thrust moment demand

The Kalman filter has the following state variables:

\(x_{LF}^G\) global LF x-position

\(y_{LF}^G\)

global LF y-position

\(\psi_{LF}\)

LF heading

\(\dot x_{LF}^G\)

global LF x-velocity

\(\dot y_{LF}^G\)

global LF y-velocity

\(\dot \psi_{LF}\)

LF heading rate

\(v_{cx}^G\)

North current velocity if current is estimated

\(v_{cy}^G\)

West current velocity if current is estimated

\(F_{Bx}^L\)

local bias x-force if bias is estimated

\(F_{By}^L\)

local bias y-force if bias is estimated

\(F_{B\psi}\)

bias moment

\(x_{HF}^L\)

local HF x-position (surge)

\(y_{HF}^L\)

local HF y-position (sway)

\(\psi_{HF}\)

HF heading (yaw)

\(\dot x_{HF}^L\)

local HF x-velocity

\(\dot y_{HF}^L\)

local HF y-velocity

\(\dot \psi_{HF}\)

HF heading rate

\(\omega_x\)

surge oscillation frequency

\(\omega_y\)

sway oscillation frequency

\(\omega_\psi\)

yaw oscillation frequency

In addition, the following important help variables are used:

\(\delta x^L\) position error surge

\(\delta y^L\)

position error sway

\(\delta \psi\)

heading error

\(\varepsilon_x^L\)

innovation in surge

\(\varepsilon_y^L\)

innovation in sway

\(\varepsilon_\psi\)

innovation in yaw

The differences between the position measurements from the vessel and estimated measurements in the Kalman filter are called innovation processes. State variables in the estimator are updated by multiplications of innovations with Kalman filter gains.

The first step is calculation of the innovations.

\[\begin{array}{ll} \varepsilon _x^G&=x_M^G-x_{LF}^G\\\\ \varepsilon _y^G&=y_M^G-y_{LF}^G\\\\ \varepsilon _x^G, \varepsilon _y^G &\rightarrow \varepsilon _x^L,\varepsilon _y^L\\\\ \varepsilon _x^L&=\varepsilon _x^L-x_{HF}^L\\\\ \varepsilon _y^L&=\varepsilon _y^L-y_{HF}^L\\\\ \varepsilon _\psi&=\psi_M-\psi_{LF}-\psi_{HF} \end{array}\]

Updating LF state variables:

\[\begin{array}{ll}\delta x_{LF}^L&=R_{11}\varepsilon _x^L\\\\ \delta \dot x_{LF}^L&=R_{12}\varepsilon _x^L\\\\ \delta y_{LF}^L&=R_{21}\varepsilon _y^L\\\\ \delta \dot y_{LF}^L&=R_{22}\varepsilon _y^L\\\\ \delta x_{LF}^L,\delta y_{LF}^L&\rightarrow\delta x_{LF}^G,\delta y_{LF}^G\\\\ \delta \dot x_{LF}^L,\delta \dot y_{LF}^L&\rightarrow\delta \dot x_{LF}^G,\delta \dot y_{LF}^G\\\\ x_{LF}^G&=x_{LF}^G+\delta x_{LF}^G\\\\ \dot x_{LF}^G&=\dot x_{LF}^G+\delta \dot x_{LF}^G\\\\ y_{LF}^G&=y_{LF}^G+\delta y_{LF}^G\\\\ \dot y_{LF}^G&=\dot y_{LF}^G+\delta \dot y_{LF}^G\\\\ \psi_{LF}&=\psi_{LF}+R_{31}\varepsilon _\psi\\\\ \dot \psi_{LF}&=\dot \psi_{LF}+R_{32}\varepsilon _\psi\\\\ \end{array}\]

Updating HF state variables:

\[\begin{array}{ll}x_{HF}^L&=x_{HF}^L+R_{14}\varepsilon _x^L\\\\ y_{HF}^L&=y_{HF}^L+R_{24}\varepsilon _y^L\\\\ \psi_{HF}&=\psi_{HF}+R_{34}\varepsilon _\psi\\\\ \end{array}\]

Updating bias estimates. If current is estimated:

\[\begin{array}{ll}\delta v_{cx}^L&=R_{13}\varepsilon _x^L\\\\ \delta v_{cy}^L&=R_{23}\varepsilon _y^L\\\\ \delta v_{cx}^L,\delta v_{cy}^L&\rightarrow\delta v_{cx}^G,\delta v_{cy}^G\\\\ v_{cx}^G&=v_{cx}^G+\displaystyle \frac{\delta v_{cx}^G}{\mathrm {max}(|v_{cx}^G|,\,v_\mathrm {typ})}\\\\ v_{cy}^G&=v_{cy}^G+\displaystyle \frac{\delta v_{cy}^G}{\mathrm {max}(|v_{cy}^G|,\,v_\mathrm {typ})}\end{array}\]

If bias force is estimated:

\[\begin{array}{ll}F_{Bx}^L&=F_{Bx}^L+R_{13}\varepsilon _x^L\\\\ F_{By}^L&=F_{By}^L+R_{23}\varepsilon _y^L\\\\ \end{array}\]

In both cases:

\[F_{B\psi}=F_{B\psi}+R_{33}\varepsilon _\psi\]

Then the frequencies of high frequency motion are updated according to (Balchen et al., 1980)

Thrust demand is computed as feedback from the LF state variables and feed-forward from measured forces (wind) and estimated forces (wave drift and current).

Feedback from LF state variables:

\[\begin{array}{ll}x_\mathrm {loc}^L,y_\mathrm {loc}^L&\rightarrow x_\mathrm {loc}^G,y_\mathrm {loc}^G\\\\ \delta x^G&=x_{LF}^G+x_\mathrm {loc}^G-x_\mathrm {ref}^G\\\\ \delta y^G&=y_{LF}^G+y_\mathrm {loc}^G-y_\mathrm {ref}^G\\\\ \delta \psi&=\psi_{LF}-\psi_\mathrm {ref}\\\\ \delta x^G,\delta y^G&\rightarrow\delta x^L,\delta y^L\\\\ \dot x_{LF}^G,\\dot y_{LF}^G&\rightarrow\delta \dot x^L,\delta \dot y^L\\\\ F_{Tx}^L&=G_{11}\delta x^L+G_{12}\delta \dot x^L\\\\ F_{Ty}^L&=G_{21}\delta y^L+G_{22}\delta \dot y^L\\\\ F_{T\psi}&=G_{31}\delta \psi+G_{32}\psi_{LF}\\\\ \end{array}\]

If current is estimated:

\[\begin{array}{ll}v_{cx}^G,v_{cy}^G&\rightarrow v_{cx}^L,v_{cy}^L\\\\ F_{Dx}^L&=d_x\|v_{cx}^L|v_{cx}^L\\\\ F_{Dy}^L&=d_y\|v_{cy}^L|v_{cy}^L\end{array}\]

If bias force is estimated:

\[\begin{array}{ll}F_{Dx}^L=F_{Bx}^L\\\\ F_{Dy}^L=F_{By}^L\end{array}\]

The calculated thrust demands are compensated for estimated and measured bias:

\[\begin{array}{ll}F_{Tx}^L&=F_{Tx}^L-F_{Dx}^L-F_{Mx}^L\\\\ F_{Ty}^L&=F_{Ty}^L-F_{Dy}^L-F_{My}^L\\\\ F_{T\psi}&=F_{T\psi}-F_{B\psi}-F_{M\psi}\\\\ \end{array}\]

The LF model consists of three discrete differential equations defined by the vessel’s mass and damping coefficients. If current is estimated:

\[\begin{array}{ll}F_{Dx}^L&=-d_x|v_{cx}^L|v_{cx}^L\\\\ F_{Dy}^L&=-d_y|v_{cy}^L|v_{cy}^L\end{array}\]

If bias force is estimated:

\[F_{D\psi}=-d_\psi|\psi_{LF}|\psi_{LF}\]

In both cases:

\[\begin{array}{ll} F_{Dx}&=-2d_xv_\mathrm {typ}\dot x_{LF}^L+F_{Bx}^L\\\\ F_{Dy}&=-2d_yv_\mathrm {typ}\dot y_{LF}^L+F_{By}^L \end{array}\]

The accelerations are expressed by:

\[\begin{array}{ll}\ddot x_{LF}^L&=\displaystyle \frac{1}{m_x}(F_{Dx}^L+F_{Tx}^L+F_{Mx}^L)\\\\ \ddot y_{LF}^L&=\displaystyle \frac{1}{m_y}(F_{Dy}^L+F_{Ty}^L+F_{My}^L)\\\\ \ddot \psi_{LF}&=\displaystyle \frac{1}{m_\psi}(F_{D\psi}+F_{T\psi}+F_{M\psi})\end{array}\]

The HF model for each degree of freedom is an adaptive harmonic oscillator. The outputs from the oscillators are estimates of wave-induced motion in the vessel coordinate system.

\[\begin{array}{ll}\ddot x_{HF}^L&=-\omega _x^2x_{HF}^L\\\\ \ddot y_{HF}^L&=-\omega _y^2y_{HF}^L\\\\ \ddot \psi_{HF}&=-\omega _\psi^2\psi_{HF}^L\end{array}\]

The accelerations are integrated to get a priori estimates of global LF-states and local HF-states.

Initial tuning of the Kalman filter-based controller may be performed as described below.

The controller gains:

\[\begin{array}{ll}\displaystyle G(1,\,1)=-\frac{2\pi }{T_x}m_x\qquad\qquad&G(1,\,2)=-2\xi _x\sqrt{G(1,\,1)m_x}\\\\ \displaystyle G(2,\,1)=-\frac{2\pi }{T_y}m_y\qquad\qquad&G(2,\,2)=-2\xi _y\sqrt{G(2,\,1)m_y}\\\\ \displaystyle G(3,\,1)=-\frac{2\pi }{T_\psi}m_\psi\qquad\qquad&G(3,\,2)=-2\xi _\psi\sqrt{G(3,\,1)m_\psi}\end{array}\]

where \(\mathrm {T_x}\), \(\mathrm {T_y}\) and \(\mathrm {T_\psi}\) are the chosen/wanted natural periods for surge, sway and yaw, respectively. The damping factors are given as \(\mathrm {\xi _x}\) , \(\mathrm {\xi _y}\) and \(\mathrm {\xi _\psi}\), and they are typically set to a value of 0.7 which indicates a damping level of 70% of the critical damping.

Depending on the natural periods of the motions (low frequency) and on the wave frequency motions, cut-off frequencies \(\mathrm {\omega _c}\) for filtering the wave frequency motions are normally chosen. Typically, the cut-off periods for surge, sway and yaw of \(\mathrm {T_c=2\pi /\omega _c}\) are given values between \(25\mathrm {s}\) and \(50\mathrm {s}\) .

Then Kalman filter coefficients (for initial tuning) may be chosen as:

\[\begin{array}{ll}R(1,\,1)=2\xi _x\omega _{c,x}\qquad&R(1,\,2)=\omega _{c,x}^2\\\\ R(2,\,1)=2\xi _y\omega _{c,y}\qquad&R(2,\,2)=\omega _{c,y}^2\\\\ R(3,\,1)=2\xi _\psi\omega _{c,\psi}\qquad&R(3,\,2)=\omega _{c,\psi}^2\end{array}\]

7.3. Velocity PI controller

When the simulator is run in manual mode, the user may use an input device, such as a joystick, to give a reference in velocity. The controller is similar to the controllers with position reference, but the velocity controller uses errors in velocity instead of errors in position.

The controller converts velocity errors into a demand for thruster forces to correct the errors. A decoupling approach allows the PI control parameters to be specified separately for surge, sway and yaw, since these are dominant when the reference point is not given in the centre of gravity.

PI characteristics

The control law of the PI (Proportional + Integral) controller is

\[\begin{array}{ll}F_{T0}&=K_P\varepsilon (t)+K_I\displaystyle \int_{0}^{t}\!{\varepsilon (\tau)}\textrm{d}{\tau}\\\\ \varepsilon (t)&=\dot x_0(t)-\dot x(t)\end{array}\]

where

\(\dot x(t)\) filtered velocity

\(\dot x_0(t)\)

desired (target) velocity

\(\varepsilon(t)\)

velocity error

\(F_{T0}\)

wanted control force from thrusters

\(K_P\)

velocity feedback gain

\(K_I\)

integral feedback gain

The PI controller may be tuned similarly to the tuning of the PID controller in Section PID controller. The damping term \(\mathrm {K_D}\) is omitted because feedback in acceleration is very sensitive to measurement noice. The wave filter is the same as the one described in Section PID controller.

8. Guidance

Guidance systems for marine vessels are used to generate a reference trajectory for the vessel to follow. The trajectories may be generated by means of waypoints.

image604
Figure 28. Waypoints

8.1. Waypoint Tracking

A waypoint system consists of a waypoint generator with a human interface.

Waypoints are usually specified using Cartesian coordinates:

\[Wpt.pos=\{(x,y,z)_{0},(x,y,z)_{1},\dot s(x,y,z)_{N}\}\]

Other waypoint properties one may specify are speed, heading, turning radius, etc.

Several algorithms exist to generate a path that the vessel can follow.

8.2. Straight lines and circular arcs

When generating paths using waypoints, one often uses straight lines and circle arcs to connect the waypoints, as seen in [tm-fig-waypoints]. The turning radius may also be given as a function of the vessel speed and a maximum rate of turn.

8.3. Line of Sight

The Line-Of-Sight concept is a geometric principle to set the desired heading towards the next waypoint. In its simplest form, the desired yaw angle is calculated by the formula:

\[\psi_d=\tan^{-1}(\frac{y_\mathrm {los}-y}{x_\mathrm {los}-x})\]

Where \((x_{\mathrm {los}},y_{\mathrm {los}})\) is the position of the waypoint and \(\mathrm {(x,y)}\) the vessels (or vessel reference ) current position. The disadvantage is that the cross-track error is not reduced to zero before the waypoint is reached. Since the cross track error is often the most important error to minimize, an addition to this simple algorithm is needed. One augmentation is to use a yaw reference equal to the bearing from the current position to the intersection of the desired track and the second line as seen in the Figure 29.

image606
Figure 29. Circle of acceptance

The switching to a new waypoint is done when the ship has entered a circle of acceptance.

Other ways to generate reference paths from waypoints may be to use spline or polynomial interpolation methods. These may be utilized in such a way that all the reference signals are continuous, which they are not when using turning circles.

8.4. Weather vaning strategy (CHREF=GBCIRCLE)

The goal of this dynamic positioning strategy is to maintain a given distance between the vessel and a point P given in the global coordinate system. The position of P ("circle center position") is given by the user as well as the distance R (circle radius) to be maintained between the vessel and P. The reference circle is defined by its center P and its radius R.

This strategy is particularly suitable for cases where the DP vessel has to maintain a constant distance from a given point while the direction of the weather may change.

At each time step, the algorithm updates the heading reference so that the vessel will be given the order to point its bow towards the circle center P. The coordinates of the reference point are updated so that it lays on the circle’s perimeter, while minimizing the distance to the vessel (or the point on the vessel).

Figure 30 shows some possible commanded positions:

image607
Figure 30. Command positions

9. Articulated structures

The articulated structure functionality has been implemented as an extension of the body type 4 of SIMO, and can for instance be used for the modelling of cranes. It consists of members rigidly connected to each other in a master slave fashion. Each slave can have a predescribed relative motion relative to the slave.

9.1. Mathematical basis

In the following, (G) refers to the global coordinate system XG, (B) refers to the local body-fixed coordinate system XB of the member in question, and (M) refers to the XB system of the member’s master body. The state of the articulated structure members is given by relative coordinates between a member body and its master. Assume that member \(\mathrm {i}\) has the following configuration:

  • Attachment point \(\mathrm {\boldsymbol{r_{AP,i}^{(B)}}=[x_{AP,i},y_{AP,i},z_{AP,i}]^{(B)}}\) for connection to the master body, referenced in the coordinates of the local coordinate system.

  • Attachment point \(\mathrm {\boldsymbol{r_{AP,i}^{(M)}}=[x_{AP,i},y_{AP,i},z_{AP,i}]^{(M)}}\) for connection to the master body, referenced in the coordinates of the master body.

  • Centre of mass \(\mathrm {\boldsymbol{r_{g,i}^{(B)}}=[x_{g,i},y_{g,i},z_{g,i}]^{(B)}}\)

  • Member mass \(\mathrm {m_{i}}\) and moment of inertia \(\mathrm {I_{c,i}}\) about the centre of gravity.

The translation of member \(\mathrm {i}\) can now be facilitated by a shift in the attachment point seen from the master reference system:

\[\Delta \boldsymbol{r_{AP,i}^{(M)}}=[\Delta x_{AP_i},\delta y_{AP_i},\delta z_{AP_i}]\]

and a rotation of member \(\mathrm {i}\) about the attachment point is given by modifying the rotation matrix \(\mathrm {\boldsymbol{\Lambda_{(L)}^{(M)}}}\) (see Equations of Motion). Any vector \(\mathrm {\boldsymbol{\chi}}\) represented in the local coordinate system by \(\mathrm {\boldsymbol{x^{(B)}}}\) is transferred to the closest master’s coordinate system or the global coordinate system by:

\[\begin{array}{ll}\boldsymbol{x^{(M)}}&=\boldsymbol{\Lambda_{(L)}^{(M)}x^{(B)}}\\\\ \boldsymbol{x^{(G)}}&=\boldsymbol{\Lambda_{(L)}^{(G)}x^{(B)}}\\\\\end{array}\]

Furthermore, the mass of the member adds directly to the total mass of the master body. For the moment of inertia we need to apply the parallel axes theorem for transferring the moment of inertia between reference systems. It is convenient to define the system transformation matrix \(\mathrm {\boldsymbol{H}}\) as given e.g. by (Fossen, 2002), which implicitly facilitates this transformation:

\[\boldsymbol{H_{(O)}^{(P)}}=\begin{bmatrix}\boldsymbol{I_{3\times 3}}&\boldsymbol{R}^T(\boldsymbol{r_{OP}})\\\boldsymbol{0_{3\times 3}}&\boldsymbol{I_{3\times 3}}\end{bmatrix}\]

Here the cross product operator is given by:

\[\boldsymbol{R}(\boldsymbol{x})=\begin{bmatrix}0&-x_3&x_2\\x_3&0&-x_1\\-x_2&x_1&0\end{bmatrix}\]

Having the inertia properties \(\mathrm {\boldsymbol{M_{i,O}^{(L)}}}\) of an articulated structure member defined relative to the local origin \(\mathrm {O}\), we first have to change the reference system by using the rotation matrix

\[\boldsymbol{M_{i,O}^{(M)}}=\boldsymbol{\Lambda_{(L)}^{(M)}M_{i,O}^{(L)}}(\boldsymbol{\Lambda_{(L)}^{(M)}})^T\]

Now the inertia properties of the system may be referred relative to a different point \(\mathrm {P}\) simply by use of the system transformation matrix:

\[\boldsymbol{M_{i,P}^{(M)}}=\boldsymbol{H_{(O)}^{(P)}M_{i,O}^{(M)}}(\boldsymbol{H_{(O)}^{(P)}})^T\]

Forces may be transformed to the new reference point by use of the same matrix:

\[\boldsymbol{F_{i,P}^{(M)}}=(\boldsymbol{H_{(O)}^{(P)}})^T(\boldsymbol{\Lambda_{(L)}^{(M)}})^T\boldsymbol{F_{i,O}^{(L)}}\]

The vector going from point \(\mathrm {P}\) (the master’s origo) to point \(\mathrm {O}\) (the local origo) is found in the master’s coordinate system as

\[\boldsymbol{r_{PO}^{(M)}}=\boldsymbol{\Lambda_{(L)}^{(M)}}(-\boldsymbol{r_{AP}^{(B)}})+\boldsymbol{r_{AP}^{(M)}}\]

If needed, the reaction forces due to motion of the crane members can now be calculated by

\[\boldsymbol{F_r}=\frac{\mathrm {d}\boldsymbol{p}}{\mathrm {d}t}=\frac{\mathrm {d}(m\boldsymbol{\nu})}{\mathrm {d}t}=m\frac{\partial \boldsymbol{\nu}}{\partial t}+\boldsymbol{\nu}\frac{\partial m}{\partial t}\]

For small velocities, the second term will be negligible compared to the first one. The forces may then be estimated based on the masses and the accelerations alone.

The velocity and acceleration of each articulated structure member can be computed by applying the expression for the time derivative of a vector in one reference system \(\mathrm {xyz}\) , expressed in an other reference system \(\mathrm {XYZ}\), (Hellan, 1969):

\[\frac{\mathrm {d}\boldsymbol{x}}{\mathrm {d}t}=\boldsymbol{\omega }\times \boldsymbol{x}+\frac{\mathrm {d}\boldsymbol{x_\mathrm {rel}}}{\mathrm {d}t}\]
Table 4. Alternatives for the definition of motion sequences
SEQTYPE Start time Stop time \(\mathrm {\boldsymbol{\Delta \chi}}\) Velocity Maximum Accel.

1

x

x

x

x

2

x

x

x

x

3

x

x

x

x

The total velocity of a body may then be found as:

\[\boldsymbol{\nu}=\boldsymbol{\nu_0}+\boldsymbol{\nu_\mathrm {rel}}+\boldsymbol{\omega }\times \boldsymbol{r}\]

where \(\mathrm {\boldsymbol{\nu_0}}\) is the translational velocity of the \(\mathrm {xyz}\) system with respect to the \(\mathrm {XYZ}\) system, \(\boldsymbol{\nu_\mathrm {rel}}\) is the translational velocity of the body with respect to the \(\mathrm {xyz}\) system and \(\mathrm {\boldsymbol{\omega }\times \boldsymbol{r}}\) is the instant circular speed of the body placed at position \(\mathrm {r}\) relative to the origin of \(\mathrm {xyz}\) . Applied to our system with member bodies we may find the global velocity \(\mathrm {\boldsymbol{v_i^{(G)}}}\) and acceleration \(\mathrm {\boldsymbol{a_i^{(G)}}}\) of member \(\mathrm {i}\) as:

\[\begin{array}{ll}\boldsymbol{v_i^{(G)}}&=\boldsymbol{v_M^{(G)}}+\boldsymbol{\Lambda_{(M)}^{(G)}v_i^{(G)}}+\boldsymbol{\omega _M^{(G)}}\times (\boldsymbol{\Lambda_{(M)}^{(G)}r_{PO}^{(M)}})\\\\&=\boldsymbol{v_M^{(G)}}+\boldsymbol{\Lambda_{(M)}^{(G)}\dot r_{PO}^{(M)}}+\boldsymbol{\omega _M^{(G)}}\times (\boldsymbol{\Lambda_{(M)}^{(G)}r_{PO}^{(M)}})\\\\\end{array}\]

and

\[\boldsymbol{\omega _i^{(G)}}=\boldsymbol{\omega _M^{(G)}}+\boldsymbol{\Lambda_{(M)}^{(G)}\omega _L^{(M)}}\]

where the local system is rotating with instantaneous angular speed \(\mathrm {\boldsymbol{\omega _M^{(G)}}}\) relative to the master system, and the master system is rotating with instantaneous angular speed \(\mathrm {\boldsymbol{\omega _L^{(G)}}}\) relative to the global system. Using the same expression for the time derivative of vectors, the acceleration is found to be

\[\begin{array}{ll}\boldsymbol{a_i^{(G)}}&=\boldsymbol{a_M^{(G)}}+\boldsymbol{\Lambda_{(M)}^{(G)}a_i^{(M)}}+\boldsymbol{\alpha _M^{(G)}}\times (\boldsymbol{\Lambda_{(M)}^{(G)}r_{PO}^{(M)}})+2\boldsymbol{\omega _M^{(G)}}\times (\boldsymbol{\Lambda_{(M)}^{(G)}v_i^{(M)}})\\\\&+\boldsymbol{\omega _M^{(G)}}\times (\boldsymbol{\omega _M^{(G)}}\times (\boldsymbol{\Lambda_{(M)}^{(G)}r_{PO}^{(M)}}))\end{array}\]

and

\[\boldsymbol{\alpha _i^{(G)}}=\boldsymbol{\alpha _M^{(G)}}+\boldsymbol{\Lambda_{(M)}^{(G)}\alpha _L^{(M)}}\]

where \(\mathrm {\boldsymbol{\alpha }}\) is the angular acceleration. These formulas can be used recursively for computing positions, velocities and accelerations of all members when the relative motion between neighbouring members is known.

9.2. Computation of motion sequences

The motion of the articulated structure members can be defined in one of the following ways (see also Table 4):

  1. By giving start time, distance to move, velocity at steady state and ramp-up acceleration

  2. By giving start time, stop time, velocity at steady state and ramp-up acceleration

  3. By giving start time, stop time, distance to move and ramp-up acceleration

Other combinations could also be thought of, but are not considered as being relevant.

Sequences with ramps of constant acceleration

We assume intervals of constant acceleration \(\mathrm {\pm a}\) , velocity \(\mathrm {\nu}\) and retardation \(\mathrm {\mp a}\) , where \(\mathrm {a}\) is a positive scalar. The upper signs then refer to the case where the velocity \(\mathrm {\nu}\) is positive and the lower signs to the case where it is negative. The following equation gives the relation between position, velocity and acceleration throughout the sequence:

\[\Delta x=\Delta t_a^2a+\Delta t_\nu\nu\]

where \(\mathrm {\Delta t_a=|\nu|/a}\) and \(\mathrm {\Delta t_\nu=(t_b-t_a)-2\Delta t_a}\) are the length of time intervals for constant acceleration and constant velocity, respectively.

In case of sequence type 1 we need to find the stop time, which is given by the expression:

\[t_b=t_a+\frac{\Delta x}{\nu}+\frac{|\nu|}{a}\]

as \(\mathrm {a}\) is always positive. Furthermore sequence type 3 gives the need of calculating the steady-state velocity, which may be found as

\[\nu=\mathrm {sign}(\Delta x)a(\Delta t-\sqrt{(\Delta t)^2-\frac{4}{a}|\Delta x|})\]

Sequences with continuous acceleration

With continuous and stepwise straight velocity functions, we get discontinuities in the acceleration function. This can be avoided by e.g. defining the startup ramps (of length \(\mathrm {\Delta }\)t) as velocity cosine functions,

\[\nu(t)=\frac{1}{2}\nu_s[1-\cos(\frac{\pi }{\Delta t}t)]\]

where \(\mathrm {\nu_s}\) is the steady-state velocity. This gives an acceleration equal to

\[a(t)=\frac{1}{2}\nu_s\frac{\pi }{\Delta t}\sin(\frac{\pi }{\Delta t}t)\]

Then the acceleration has a maximum value equal to \(\displaystyle a_\mathrm {max}=|(1/2)\cdot \nu_s(\pi /\Delta t)|\) which is a factor \(\mathrm {\pi /2}\) larger than what we get with constant acceleration ramps. In order to ensure smooth start and stop of sequences we may just replace the straight velocity ramps with the smooth ones, keeping the other parameters the same. In computations using equations Equation (245) to Equation (247) we must replace \(\mathrm {a}\) by \(2a/\pi\) .

Multistep sequences

Each sequence is defined with a start time, stop time and velocity. In case the user specifies several motion sequences for the same body, the sequences will just be added. This is the same as defined for the winch sequence format. See Figure 32 for an example.

image647
Figure 31. Example of a simple sequence
image648
Figure 32. Multistep sequence