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

Model characteristics and limitations

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

Model description

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

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

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

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

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

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

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

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

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

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

Computing the mass properties of a tank

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

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

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

where:

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

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

Computing the position of the fluid free surface

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

Permeability factor

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

Damaged tank

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

Limitations

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

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

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

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

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

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

where:

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

The force is calculated according to the following formula

where:

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

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

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

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

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

The average drift force in irregular waves is

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

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

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

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

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

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

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

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

Wave drift forces on a floating body depends on the relative velocity between the body and the ambient current. By linearizing the wave drift force about a relative velocity \(\mathbf{v}_{r,0}\) we can write:

The last term is here referred to as the wave drift damping force, and \(\mathbf{B}\) is a time varying matrix computed from wave drift damping coefficients and a wave realization. SIMA/SIMO can either linearize about zero, \(\mathbf{v}_{r,0} = \mathbf{0}\), or about the current velocity, \(\mathbf{v}_{r,0} = \mathbf{v}_c\). Wave drift damping coefficients computed by WADAM is based on linearization about \(\mathbf{0}\). Wave drift damping coefficients originating from Aranhas formula will on the other hand be more accurate when linearizing about \(\mathbf{v}_c\). When a current is present and the linearization is done about \(\mathbf{0}\), the term "wave drift damping" is misleading since the current will also introduce additional excitation implicitly. When the linearization is done about \(\mathbf{v}_c\), this additional excitation must be included explicitly as an additional excitation term which modifies the zero-current wave drift force (otherwise SIMA will give a warning). This is relevant for Aranhas model, as further described in Wave-current interaction correction of QTFs using extended Aranha method.

When the linearization is done about \(\mathbf{0}\), the wave drift damping force in the three horizontal degrees of freedom can be written as,

where \(v_{ic}\) and \(v_{ib}\) is the current velocity and body velocity, respectively, both decomposed in the low-frequency body related coordinate system. This coordinate system have a vertical z-axis, with x and y-axes following the low-frequency yaw motion of the body.

When the linearization is done about \(\mathbf{v}_c\), the wave drift damping force in the three horizontal degrees of freedom can be written as,

The coefficients in the wave drift damping term is time varying and can be written as:

which is completely analogous to equation Equation (11) for the slowly varying second-order wave force model. The calculation of these terms are carried out in a manner as described in Simplified wave drift force model, based on Newmans method. As a consequence, only the diagonal terms \(B_{mm,ij}\) are needed, and it is these which are referred to as the (frequency dependent) wave drift damping coefficients.

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

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

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.

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

The current components inline and transverse to the wave direction becomes:

where \(\theta_w = \theta - \beta\).

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

Here, the frequency of encounter can be written as,

and the direction of encounter can be written as:

It should here be noted that the use of \(\omega_{e}\) and \(\beta_{e}\) on the left hand side of Equation (42) is a consequence of the fact that the wave is measured in a stationary point, as discussed in Wave-current correction of second order wave drift force and wave drift damping using an extended Aranhas formula. The discrete frequencies and directions used in the zero-current wave-body interaction analysis (e.g. by WADAM) to compute \(Q_{0}\) is further assumed to correspond to frequencies and direction of encounter when a current is introduced. The corresponding wave frequencies and directions, as measured in a point following the current, denoted \(\omega\) and \(\beta\) and used in the calculation of the scaling factor \(B\), therefore need to he computed by solving Equation (43) and Equation (44) for the unknown \(\omega\) and \(\beta\).

\(C_p\) and \(C_g\) appearing in Equation (43) and Equation (44) is the phase and group velocity of the waves (finite or infinite water depth), respectively.

The scaling factor is based on the geometric average of the Aranha scaling factors for the two interacting wave components:

The Aranha scaling factor can be written as,

where,

The dispersion relation in the finite water depth case reads,

which gives a phase velocity and group velocity:

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

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

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

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

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

After some calculation, we obtain:

and

where:

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

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

The obtained wave drift damping coefficients are:

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

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

where:

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

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

where:

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

where:

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

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

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

For strip model no. 1:

For strip model no. 2:

where:

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

where

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

while the yaw moment is given as

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Three coordinate systems are used,

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

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

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

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

Two parameters are used to define different calculation options:

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

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

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

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

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

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

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

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

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.

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.

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.

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