1. Station-keeping forces

1.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 1. 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 the equation for the dynamic equilibrium equation and hydrodynamic forces 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 2. Equivalent geometric model.

The coefficients for generalized line damping and stiffness \(\mathrm {c^*}\) and \(\mathrm {k^*}\) in Equation (3) 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 force_models/body.adoc#tm_eq_4_99 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 (4) from

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

Equation (6) 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.

1.2. Thrusters

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

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

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

1.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 (10) 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 3.

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

1.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 4 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 4. Feedback control system for thruster speed of rotation.

In Figure 4 \(\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 (10). 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 6. The effective moment of inertia in Equation (12) 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 5, which shows how the ratio \(\mathrm {n/n_0}\) varies after a step change in the reference \(\mathrm {n_0}\).

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

Figure 5 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 6. 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.

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

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

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

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

1.4. Force-elongation with/without hysteresis

image440
Figure 8. 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 9 a hysteresis model is shown. This model is obtained when \(\mathrm {p=0}\).

In Figure 10 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 9. Hysteresis model, p = 0.
image444
Figure 10. 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.

1.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 11. 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 11. Longitudinal distance.

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

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.

1.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 12. Any fender type can be selected for the shown alternatives.

image452
Figure 12. Fender alternatives, examples.

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

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

image460
Figure 13. 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.

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

1.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 14. Bumper contact.