Force models

Assume that the equations of motion can be written

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

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

In the frequency domain, the equation is written

or

Using the following relations,

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

Using the inverse Fourier transform gives

Physically, values of \(\mathrm {\boldsymbol{h}(t-\tau)}\) for, \(\mathrm {t<0}\), i.e. before the "experiment" started, are zero. Also, causality implies that \(\mathrm {\boldsymbol{h}(t-\tau)\equiv0}\) for \(\mathrm {\tau>t}\).

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

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

or similarly

Note that this transform pair is identical to the Fourier transform defined in General Assumptions and Notation.

Using the fact that \(\mathrm {\boldsymbol{c}(\omega )=\boldsymbol{c}(-\omega )}\) and \(\mathrm {\boldsymbol{a}(\omega )=\boldsymbol{a}(-\omega )}\) gives

From causality, \(\mathrm {\boldsymbol{h}(\tau)=0}\) for \(\mathrm {\tau<0}\); the process can not have any memory effect of the future. This means that the two parts in the integral, Equation (15) must be opposite for \(\mathrm {\tau<0}\) and identical for \(\mathrm {\tau>0}\), or mathematically:

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

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

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

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

The relations in Equation (17) and 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.

As an alternative to solving the whole differential equation Equation (1) in time domain by use of the retardation function, the motions can be separated in a high frequency part and a low frequency part.

The high frequency motions are solved in the frequency domain, which require the motions to be linear responses to waves. This means that \(\mathrm {\boldsymbol{D}_2}\), the quadratic damping, is set to be zero and \(\mathrm {\boldsymbol{K}}\) constant.

The exciting force is separated in a high frequency part, \(\mathrm {\boldsymbol{q}^{(1)}}\) and a low frequency part, \(\mathrm {\boldsymbol{q}^{(2)}}\).

The position vector can then be separated,

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

In the frequency domain, using Equation (34), we have

where \(\mathrm {\boldsymbol{H}_1}\) is the first order transfer function between excitation force and wave elevation and \(\mathrm {\boldsymbol{X}_{HF}}\) is the first order transfer function between motion and wave elevation.

The low frequency motions are solved in time domain for which the dynamic equilibrium equation is written

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

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:

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.

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.

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

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

Wave forces may be split into three contributions:

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.

where:

ringing moment

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

where:

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.

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

where:

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

where:

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

where the time average is

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

Assuming that

the slowly-varying force can be approximated by

The double sum is circumvented using the relations

and

Thus, the slowly varying force may be written

where

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

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

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

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

\(\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

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

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

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

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.

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

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 (61) and Equation (62) can be expressed as

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 (63) and Equation (64) 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"

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.

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.

Two different methods are possible:

Simplified method

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

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

and

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

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.

Newman’s method

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

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 relative velocities in direction surge, sway and yaw, respectively. The wave drift damping force coefficients in the matrix above can be expressed by

where:

)

direction (i) due to velocity in direction (j) in bi-chromatic waves with frequencies ( _m ), ( _n ) and directions ( _m ), ( _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.

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

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

where:

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

where:

The dynamic equilibrium equation is written

where:

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

where:

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:

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

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

For structures with closed compartments, the modelling also includes:

Soil buoyancy force

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

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.

Pressure buildup in a closed compartment

The volume flow is calculated as

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

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

where:

cavity

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

where:

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.

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:

Inertia:

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.

Soil friction force

The soil reaction force is calculated as:

Downward motion:

Upward motion:

where

resistance, depth-dependent

), depth-dependent

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

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:

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.

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

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

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

where

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.

The coefficients for generalized line damping and stiffness \(\mathrm {c^*}\) and \(\mathrm {k^*}\) in Equation (116) 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 (101) on incremental form. The equilibrium of the line at time \(\mathrm {t_i}\) can be expressed

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

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

Equation (119) 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

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.

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

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.

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

where:

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

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.

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:

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

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:

The axial motion into the cylinder is found from

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

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.

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

The implemented model will give:

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.

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.

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:

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

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

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

which can be rewritten to

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.

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

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

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.

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.

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.

An example with use of coupling points and a guidewire is shown in 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.

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:

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:

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:

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:

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

where :

The effective axial stiffness is given by :

where :

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

Material damping is included as:

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.

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.

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:

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:

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:

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:

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:

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:

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

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:

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

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.

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.

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.

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

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

where

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

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.

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.

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

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:

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

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:

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

The output from the routines are:

The Kalman filter has the following state variables:

In addition, the following important help variables are used:

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.

Updating LF state variables:

Updating HF state variables:

Updating bias estimates. If current is estimated:

If bias force is estimated:

In both cases:

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:

If current is estimated:

If bias force is estimated:

The calculated thrust demands are compensated for estimated and measured bias:

The LF model consists of three discrete differential equations defined by the vessel’s mass and damping coefficients. If current is estimated:

If bias force is estimated:

In both cases:

The accelerations are expressed by:

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.

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:

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:

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

where

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.

A waypoint system consists of a waypoint generator with a human interface.

Waypoints are usually specified using Cartesian coordinates:

Other waypoint properties one may specify are speed, heading, turning radius, etc.

Several algorithms exist to generate a path that the vessel can follow.

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.

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:

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.

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.

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:

The translation of member \(\mathrm {i}\) can now be facilitated by a shift in the attachment point seen from the master reference system:

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:

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:

Here the cross product operator is given by:

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

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:

Forces may be transformed to the new reference point by use of the same matrix:

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

If needed, the reaction forces due to motion of the crane members can now be calculated by

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):

The total velocity of a body may then be found as:

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:

and

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

and

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.

The motion of the articulated structure members can be defined in one of the following ways (see also Table 4):

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:

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:

as \(\mathrm {a}\) is always positive. Furthermore sequence type 3 gives the need of calculating the steady-state velocity, which may be found as

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,

where \(\mathrm {\nu_s}\) is the steady-state velocity. This gives an acceleration equal to

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 (220) to Equation (222) 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.