1. Body forces
1.1. Mass forces
The mass forces are described by the product \(\mathrm {\boldsymbol{m\boldsymbol{\ddot {x}}}}\) where \(\mathrm {\boldsymbol{m}}\) is the (constant) mass matrix.
1.2. Timedependent mass
The data group TimeDependent 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:

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

Filling liquid into one or a set of tanks with specified geometry. Free surface effects can then be accounted for. In this case, one TDM may be modelled:  in terms of
PORTIONS
. Up to 5PORTIONS
are allowed. EachPORTION
has a simple geometry (CONE
orBOX
). The TDM’s location and orientation relative to the body (vessel) needs to be specified.  in terms of a set of tanks. Each tank are represented by a mesh of triangular panels described in a STL (STereoLithography) file. The maximum number of ballast tanks is unconstrained. 
A mass that moves relative to the body (future extension)
In all cases, the effect of the force terms \(\frac{\mathrm{d}m}{\mathrm{d}t}\cdot v\) and \(\frac{\mathrm{d}I}{\mathrm{d}t}\cdot \omega\) are neglected (mass \(\mathrm {m}\) and inertia \(\mathrm {I}\)). This means that forces that are due to relative velocity between the TDM and the body are assumed to be negligible. In each time instant in a simulation, the mass matrix of the body will be updated and the gravity force due to the TDM will be calculated. In cases with tank filling, the centre of gravity of the liquid will also be updated. No sloshing forces are accounted for in the model and the surface of the liquid is always horizontal.
1.2.1. Operational relevance
A model of timevarying mass is useful for simulation of operations involving ballasting or deballasting, 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 timedependent mass. If a TDM without geometry is used, the TDM mass may be negative.
1.2.2. Point mass model
Model characteristics and limitations
The model enables specification of one or more point masses on a body. Each point mass may vary with time. The location of the point is fixed in the local body coordinate system. The mass point can be a mass point or a tank with specified geometry. In both cases the mass is taken as an addition to the system mass specification.
the mass will always contribute to gravity force. 
Model description
The instantaneous point mass, \(\mathrm {m(t)}\), is found by using the prescribed mass rate. During static analysis the specified value at time zero \(\mathrm {m(t=0)}\), is used.
A point mass, \(\mathrm {m(t)}\), located at coordinates \(\mathrm {(x,y,z)}\) in local (body fixed) system, will give the following change in the body mass matrix.
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.
1.2.3. Tank with simple geometry
When a TDM with geometry is specified and the volume of one portion
(portion B) shall be subtracted from portion A (IADD= 1
), then no
part of portion B may be outside portion A. There is no check in SIMO
that this event does not occur. If the geometry of a TDM is defined by
adding portions, the portions must not be overlapping (the intersection
must be zero).
For a TDM with geometry, trim of the TDM in the range \(\mathrm {(89\deg<}\) trim \(\mathrm {<91\deg)}\) is not allowed. If this becomes a problem, the orientation of the tank xaxis relative to the vessel should be modified.
1.2.4. Set of tanks with general geometry (ballast tank system)
When a TDM is defined as a ballast tank system, each tank geometry is represented by a triangular mesh described in a STL file. The number of tanks is unconstrained.
Computing the mass properties of a tank
A tank containing an amount of fluid contributes to the gravity force
and to the mass matrix of the body. In order to compute this
contribution, we assume that we know where the free surface of the fluid
is relative to the tank. SIMO
first proceeds to sort each panel of the
mesh as follows:

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

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

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

if all the three vertices are under the free surface line, the panel is included in the "contributing panels".
Thanks to this cutting procedure, the mesh can be coarse as long as it describes the geometry correctly.
For each "contributing" panel, the following integrals are computed analytically:

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

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

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

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

\(\mathrm {\boldsymbol{X}}\) is the position vector in the Earthfixed coordinate system

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

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

\(\mathrm {V_{m2}}\) is the contribution of the panel to the total second volume moment
By summing up the contribution of each "contributing panel", the total fluid volume is computed, as well as the position of the center of volume. The two latter quantities are used to compute the total gravity force and moment, using the density of the fluid contained in the tank.
The second volume moments are used to build the mass matrix of the tank, which is transposed to the bodyfixed 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.
1.3. Hydrodynamic reaction forces
Based on first order potential wave theory, the hydrodynamic reaction forces are described by addedmass coefficients, \(\mathrm {\boldsymbol{A}}\), and damping coefficients, \(\mathrm {\boldsymbol{C}}\). Both addedmass and damping are frequencydependent. They are included in the retardation functions, as described in Solution by convolution integral or the firstorder motion transfer functions as described in Separation of motions.
1.4. Damping forces
Linear damping is expressed by \(\mathrm {\boldsymbol{D}_1\boldsymbol{\dot {x}}}\), where \(\mathrm {\boldsymbol{D}_1}\) is the linear damping matrix.
Quadratic damping is expressed by \(\mathrm {\boldsymbol{D}_2\boldsymbol{f}(\boldsymbol{\dot {x}})}\) where \(\mathrm {\boldsymbol{D}_2}\) is the quadratic damping matrix and \(\mathrm {\boldsymbol{f}}\) is a vector function where each element is given by \(\mathrm {f_i=\dot {x}_i\dot {x}_i}\).
1.5. Hydrostatic stiffness
1.5.1. Linear stiffness matrix model
The force is expressed by \(\mathrm {\boldsymbol{K(xx_0)}}\) where \(\mathrm {\boldsymbol{K}}\) is the hydrostatic stiffness matrix and \(\mathrm {\boldsymbol{x_0}}\) the stiffness reference position.
1.5.2. Nonlinear stiffness model
When using this formulation, the hydrostatic force model will return the complete buoyancy force (and not the restoring forces as it is the case when using linear stiffness matrix). The buoyancy force is computed based on the part of the hull geometry volume which is under the calm water line (waves are ignored).
The geometry of the hull is described by a STL (STereoLithography) file.
The mesh of the geometry is then based on triangular panels. At each
time step, SIMO
will update the position of the body and proceed to a
cut of the mesh as follows:

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

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

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

if all the three vertices are under the calm water line, the panel is included in the "contributing panels".
Thanks to this cutting procedure, it is not necessary to refine the mesh about the calm water line. The mesh can be coarse as long as it describes the geometry correctly.
For each "contributing" panel, the following integrals are computed analytically:

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

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

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

\(\mathrm {\boldsymbol{X}}\) is the position vector in the Earthfixed coordinate system

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

\(\mathrm {V_{m1}}\) is the contribution of the panel to the total volume center
By summing up the contribution of each "contributing panel", the total submerged volume is computed, as well as the position of the center of buoyancy. The two latter quantities are used to compute the total buoyancy forces and moments which are transformed to the bodyfixed coordinate system.
1.6. Aerodynamic forces
The wind force is calculated based on the instantaneous wind and body velocities.
1.6.1. Classic wind coefficients
The force is calculated according to the following formula
where:
\(j\) 
degree of freedom 
\(C\) 
wind force coefficient for the instantaneous relative direction 
\(v\) 
relative wind speed between body and wind 
\(\alpha\) 
relative velocity direction in local coordinate system 
The low frequency body motion is used to calculate the relative wind speed. Optionally, the velocity of the body may be omitted when calculating relative velocity.
1.6.2. Position dependent wind coefficients
This model is adapted for situations where large changes of positions and attitudes are expected during the simulation (a damaged condition for example).
It is possible to define the wind coefficients as a function of:

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

\(\mathrm {z}\) : vertical position of the body in Earthfixed coordinate system

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

\(\mathrm {\theta }\) : pitch angle of the body
In that case, the coefficients are interpolated for a given body position and a relative wind direction. The force is computed as in equation Equation (3).
1.7. Wave excitation forces
Wave forces may be split into three contributions:

firstorder forces that oscillate with wave frequencies

secondorder mean, rapidly varying, and slowly varying wave drift forces.

higherorder ringing forces.
It is assumed that the wave forces can be formally expressed by the following equation. The two first terms are a Volterra series truncated after the quadratic term, see (Hasselmann, 1966) and (Dalzell, 1976), while the last term is the timedependent ringing load.
where:
ringing moment
\(q^{(1)}(t)\) 
timedependent firstorder wave force 
\(q^{(2)}(t)\) 
timedependent secondorder wave force 
\(q^{(R)}(t)\) 
timedependent thirdorder ringing force or third and fourthorder ringing moment 
\(h^{(1)}\) 
linear impulse response function 
\(h^{(2)}\) 
secondorder impulse response function 
It is assumed that the kernels \(\mathrm {h^{(1)}}\) and \(\mathrm {h^{(2)}}\) are smooth, absolutely integrable and possess Fourier transforms as follows
where:
\(H^{(1)}\)  firstorder transfer function 

\(H^{(2)}\) 
secondorder transfer function 
The physical interpretation of firstorder transfer functions is well established in linear response analysis.
The interpretation of the secondorder transfer function is not obvious. A simple example taken from (Dalzell, 1976) shows that the secondorder 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 secondorder forces is relatively small compared to firstorder 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 secondorder 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 firstorder forces.
The ringing response of TLPs with high natural heave and roll/pitch periods \((35~\mathrm {s})\) or large GBSs in deep water with high fundamental modal period \((3.56~\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.
1.8. Firstorder wave excitation forces
First order wave forces are described in the frequency domain as a transfer function between wave elevation and force
where:
\(\boldsymbol{H}^{(1)}\) 
complex first order transfer function 
\(\tilde{\zeta} (\omega)\) 
complex harmonic wave component 
1.9. Secondorder wave forces
The secondorder wave force for unidirectional irregular waves can be expressed in the form
where:
\(\tilde{\zeta}_m\) 
complex Fourier component of sea surface elevation with frequency \(\omega_m\) and direction \(\beta\) 
\(H_{mn}^{(2+)} \equiv H^{(2+)} (\omega_m, \omega_n, \beta, \beta)\) 
secondorder transfer function for the sumfrequency force 
\(H_{mn}^{(2)} \equiv H^{(2)} (\omega_m, \omega_n, \beta, \beta)\) 
secondorder transfer function for the differencefrequency force 
SIMO can include full second order forces using the efficient method presented in (Stansberg, 1994), which circumvents (without loss of accuracy) the double sum on every time step. Full second order forces is not yet available for multidirectional (shortcrested) seas. The simplified wave drift force model (presented in the next section) can also be used for multidirectional seas.
1.10. Simplified wave drift force model
The Newman method can be used to calculate secondorder 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 secondorder transfer function.
Neglecting the sumfrequency force in equation Equation (11), the slowlyvarying secondorder force, also referred to as the wave drift force, is expressed by
where the time average is
The main contribution to the slowlyvarying force is assumed to be associated with frequency combinations for which, \(\mathrm {\omega _m\approx \omega _n}\) or
Assuming that
the slowlyvarying 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 lowpass 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 pregeneration 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 differencefrequency force transfer functions in both parts of Equation (10), the force is
As in Equation (13) 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 (15), can then be seen from
It is possible to precalculate wave drift forces for a number of headings before time domain simulation starts, and interpolate between these time series in time domain.
It is possible to correct the wave drift force coefficients for the presence of wavecurrent interaction effects using Aranhas method. This is described in more detail in section Wavecurrent correction of second order wave drift force and wave drift damping using an extended Aranhas formula.
1.10.1. Second order wave forces in shortcrested seas
The following method is implemented in order to include the effect of shortcrested seas on wave drift forces.
The average drift force in irregular waves is
in which the average, or unidirectional 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 shortcrested 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 )}\).
1.10.2. Second order wave forces calculated with FFT and COS method in shortcrested seas
Long crested waves
The calculation is based on the complex wave Fourier components \(\mathrm {X(\omega )}\) and the real drift force coefficients \(\mathrm {C(\omega )}\) for the actual direction.
The drift force is calculated from the complex coefficients
FFT and COSINE simulation method will give the same result.
FFT method shortcrested seas
For each frequency, SIMO will calculate \(\mathrm {Z(\omega )}\) from equivalent coefficients \(\mathrm {X(\omega )}\) and \(\mathrm {C(\omega )}\)
1.11. Wave drift damping
Two different methods are possible:

Newman’s method. This method uses frequency dependent wave drift damping coefficients based on the change of QTF diagonal (wave drift force coefficients) with change in relative velocity.

A simplified model based on the derivatives of nondimensional average wave drift forces with respect to relative velocity. This method uses a frequencyindependent wave drift damping coefficient and should be used with care.
1.11.1. Newman’s method
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 zerocurrent wave drift force (otherwise SIMA will give a warning). This is relevant for Aranhas model, as further described in Wavecurrent 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 lowfrequency body related coordinate system. This coordinate system have a vertical zaxis, with x and yaxes following the lowfrequency 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 secondorder 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.
1.11.2. Simplified method
The wave drift damping model is based on the derivatives of nondimensional 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 sinksource 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.
1.12. Wavecurrent correction of second order wave drift force and wave drift damping using an extended Aranhas formula
The following subsections documents the calculation of

wave drift forces for a specified current condition (velocity and direction), and,

wave drift damping coefficients for a specified current condition (velocity and direction)
using Aranhas method (Aranha, J. A. P. (1994, 1996)). Aranhas method was extended to finite waterdepth by Malenica et. al. (1995). We are also using an heuristic extension of Aranhas method for the offdiagonal elements of the QTF, as explained in the following.
In this section we use \(\omega_e\) and \(\beta_e\) to denote the frequency of encounter and the direction of encounter, respectively. SIMO is intended for simulation of stationary bodies, and the encounter effect is therefore due to the presence of a current.
Waves are always measured in a stationary point in SIMO, and as a consequence the wave spectrum models defined in Waves takes frequency and direction of encounter as input, although here the symbols \(\omega\) and \(\beta\) is used without subscript \(e\). 
Moreover, all transfer functions which are used by SIMO have a frequency parameter which corresponds directly to the frequency of oscillation of the corresponding response, and therefore need to be interpreted as a frequency of encounter when a current is present. Also here, the symbol \(\omega\) is generally used without subscript \(e\), except for in this section.
1.12.1. Wavecurrent interaction correction of QTFs using extended Aranha method
Consider a stationary floating body subjected to a current speed \(U\) with a goingto direction of \(\theta\) relative to the body xaxis and waves propagating in direction \(\beta\) relative to the body xaxis. 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 wavecurrent 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 Wavecurrent correction of second order wave drift force and wave drift damping using an extended Aranhas formula. The discrete frequencies and directions used in the zerocurrent wavebody 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:
1.12.2. Wave drift damping coefficients calculation
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:

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

wave directions

wave frequencies

water depth

current velocity

current direction
The obtained wave drift damping coefficients are:

B11, B21, B61 from Equation (56)

B12, B22, B62 from Equation (57)
They represent the change of wave drift forces in X and Y and moment in Z due to a velocity in X and Y.
1.13. Nonlinear modification of restoring and wave forces
In SIMO the model for waveinduced 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 4node panel cells. At each time step, the location of each cell centroid is calculated in earthfixed 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 degreeoffreedom generalized rigid body forces and moments in the bodyfixed coordinate system.
Nonlinear correction of wave force and restoring force on body. Quasistatic pressure is computed over the wetted panels, and the resulting timevarying forces and moments are added as corrections to those obtained using frequencydependent firstorder 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 firstorder theory, if available.
However, the abovementioned 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 earthfixed coordinate system. \(\mathrm {\boldsymbol{r}}\) is the position vector. \(\mathrm {S_0}\) denotes the mean wetted surface. \(\mathrm {SS_0}\) is the difference between the instantaneous wetted surface and the mean wetted surface. The last terms in Equation (63) and Equation (64) 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 (65) and Equation (66) is taken care of in the program. If secondorder 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"
1.14. Diffracted waves
Diffracted wave transfer functions for several diffracted wave points can be given for a body. These transfer functions are given in a similar manner as e.g. the wave to motion transfer functions. Diffracted wave transfer functions for wave elevation and global x, y and zvelocity 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.

For body types 1 or 2 (large volume bodies), the user may specify one or more "slender elements" and "fixed body elements" under data group
DISTRIBUTED ELEMENT FORCE
. Each of these elements may feel different diffracted waves by specifying a "body name" and "diffracted wave points". For slender elements, each end of the element refers to one "diffracted wave point". If the point numbers are identical, all parts of the element will "feel" the same wave particle motion series. Otherwise, linear interpolation is performed on the particle motion series. The special case of "fixed body elements" refers to one "diffracted wave point" only. 
For body type 3 (3dof body), the effect of wave diffraction may be included in data group
SMALL BODY HYDRODYNAMIC DATA
by referring to a body and a "diffracted wave point".
Diffracted waves are not included in other force models. This means that in order to include diffraction effects in these models, e.g. "Second order wave drift coefficients" the coefficients have to be modified.
1.15. Current drag forces
1.15.1. Classic current coefficients
All body current forces are computed using the current velocity at the surface (\(\mathrm {z=0}\)). Traditionally, the viscous hull surge and sway force and yaw moment have been calculated based on current coefficients and the instantaneous magnitude of the translational relative velocity between the vessel and the fluid. The current drag forces are then expressed by:
where:
\(k\) 
degree of freedom 
\(C_1\) 
linear current force coefficient 
\(C_2\) 
quadratic current force coefficient 
\(u\) 
relative velocity between low frequency body velocity and current velocity 
\(\alpha\) 
relative angle between direction of low frequency body velocity and current velocity 
\(v_1,v_2\) 
current velocity components at the surface 
\(\dot{x}_1,\dot{x}_2\) 
vessel velocity components in the vessel coordinate system 
Usually, the linear coefficient, \(\mathrm {C_1}\) is not used. The model does not account for the effect from the yawinduced crossflow. This crossflow has normally been included as a separate yaw damping term. The crossflow based on a pure relative translational velocity and yaw velocity is illustrated in Separation of crossflow….
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:
\(C_2^{(6)}(\alpha)\) 
Current coefficient in yaw as a function of the incident heading, \(\alpha\) 
\(b_{66}\) 
Quadratic damping coefficient in yaw due to a pure yaw velocity 
\(\dot{x}_6\) 
Yaw velocity 
Normally, the current coefficient in yaw, \(\mathrm {C_2^{(6)}(\alpha )}\), also includes the Munk moment which is an effect that can be derived from potential theory. When quadratic current coefficients are included in the simulation model, the Munk moment is excluded from the equations of motion. The Munk moment, denoted \(\mathrm {q_M(\alpha )}\), is given as
where:
\(A_2\)  Added mass in sway 

\(A_1\) 
Added mass in surge 
It is, however, quite clear that the current velocity, the vessel translational velocity and the yaw velocity will influence the distribution of the crossflow along the hull. Theoretically, the sum of these contributions will result in a linear crossflow 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 yawinduced crossflow 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 crossflow 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:
\(c_d(x)\) 
Nondimensional drag coefficient of the vessel section 
\(C_d(x)\) 
Drag coefficient of the vessel section 
\(T(x)\) 
Vessel draught at the vessel section 
\(l\) 
Vessel length 
\(\rho\) 
Water density 
The drag force of a strip, \(\mathrm {f_s(x)}\), located at a distance from the centre of motion is then given as
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 twodimensional 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.
1.15.2. Position dependent quadratic current coefficients
This model is adapted for situations where large changes of positions and attitudes are expected during the simulation (a damaged condition for example). This model is only available for quadratic current drag coefficients (not available for linear current coefficients).
It is possible to define the current coefficients as a function of:  \(\mathrm {\alpha }\) : relative velocity direction in local coordinate system  \(\mathrm {z}\) : vertical position of the body in Earthfixed 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.
1.16. Distributed element force
1.16.1. Physical relevance
The distributed element force model is applicable to two different modelling features:

Long, slender elements

Concentrated, fixed elements (with zero extension)
Long, slender elements can be used to model jacket legs and bracings or spool pieces. A slender object such as a spool piece will normally have a 3D geometry consisting of a number of straight elements with different orientation, and can be modelled by a set of slender elements. End connectors can be modelled either by a short element or by a concentrated, fixed element.
The hydrodynamic properties of a spoolpiece can be assumed to be rotation symmetric. However, in order to handle partly covered or stiffened spoolpieces and other general long structures, an assumed rotation symmetry or assumed cylindrical shape would be too limiting for the model.
A fixed connection between a body and a small element, which can be assumed to give concentrated hydrodynamic forces, is useful for the modelling of various configurations, such as:

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

imaginary bodies attached to a floating vessel, in order to include nonlinear deviations from force coefficients calculated by linear wave diffraction theory programs.
As an alternative model, a small 3dof 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.
1.16.2. Element definition
Both element types may have any orientation. The main body must be of
BODY TYPE 1
, i.e. 6dof body.
Smallbody theory is used to calculate forces both for the slender
element and the concentrated fixed element. Stiff connection implies
that all the "subbody" (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.
1.16.3. Coordinate systems
Three coordinate systems are used,

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

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

XS
, local strip coordinate system,XS(1)
,XS(2)
,XS(3)
All coordinate systems are orthogonal and righthanded. 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 xyplane and the XS(3)
axis is then defined from the
definition of orthogonal righthanded coordinate systems.
1.16.4. Calculation of the mass matrix, M
Contribution from structural mass
The contribution from each of the slender elements to the mass matrix of the main body is calculated and added, in local body coordinates. The contribution from each element is calculated by adding the contributions from each strip, which has length \(\mathrm {L}\), mass per unit length, \(\mathrm {m}\), and coordinates \(\mathrm {(x,y,z)}\) in the local body coordinate system. The contribution from each strip within each element to the total mass matrix of the body will be
where \(\mathrm {\boldsymbol{I}}\) is a 3x3 unit matrix and
The mandatory specified BODY MASS DATA
for the main body can be zero
(but not necessarily), and will be added to the contributions from the
slender elements.
Contribution from added mass
The contribution from the added mass will be an extension to the above, since the added mass of each strip is different in the 3 directions.
Assume the following parameters:
\(\boldsymbol{V}\)  Linear acceleration vector (main body) 

\(\boldsymbol{W}\) 
Rotation acceleration vector (main body) 
\(\boldsymbol{R}\) 
Radius vector to the strip 
\(\boldsymbol{M}_a\) 
Added mass matrix for the strip 
The local acceleration of the strip will be
The local force and moment can be expressed by
Further expansion yields
The resulting added mass contribution is
The contribution from all strips will be added to the body mass matrix.
Depthdependent hydrodynamic coefficient corrections may be specified for both slender and fixed elements. Added mass data for the main body may be zero (but not necessarily), and will be added to the contributions from the slender elements.
1.16.5. Calculation of external load, F
The model for external load on a slender element strip consists of three contributions:
\(\boldsymbol{F}_B\)  buoyancy forces 

\(\boldsymbol{F}_W\) 
wave forces 
\(\boldsymbol{F}_S\) 
slamming forces 
The resulting force on a slender element is accumulated force contribution from each strip. The body force is the sum of contributions from all elements.
Gravity forces
The gravity force acts in the global zdirection and is written
Buoyancy forces
The buoyancy force acts in the global zdirection, through the center of buoyancy. The load vector for a strip, expressed in the global coordinate system, can be written,
Wave forces
The wave force on a strip acts through the center of buoyancy. The load vector for a strip, expressed in the local strip coordinate system can be written
where the quadratic drag force is computed as,
and,
and we have that:
\(V_S\) 
is submerged volume per length unit, calculated to \(z = 0\) 
\(\boldsymbol{a}_S =\begin{bmatrix} a_{sx} & a_{sy} & a_{sz} \end{bmatrix}^T\) 
is water particle acceleration components in local strip coordinate system 
\(\boldsymbol{U}_{rel,S} = \begin{bmatrix} U_{rel,sx} & U_{rel,sy} & U_{rel,sz} \end{bmatrix}^T\) 
is the relative water particle velocity components in local strip coordinate system (relative to strip velocity) 
\(\boldsymbol{v}_S =\begin{bmatrix} v_{sx} & v_{sy} & v_{sz} \end{bmatrix}^T\) 
is wave induced water particle velocity components in local strip coordinate system 
\(\boldsymbol{\dot{X}}_S =\begin{bmatrix} \dot{X}_S & \dot{Y}_S & \dot{Z}_S \end{bmatrix}^T\) 
is strip velocity components in local strip coordinate system 
\(\boldsymbol{U}_S =\begin{bmatrix} U_{sx} & U_{sy} & U_{sz} \end{bmatrix}^T\) 
is current flow velocity in local strip coordinate system 
and \(\qquad \qquad \quad \circ\) 
denotes elementwise multiplication 
The first and second terms in equation Equation (87) contains the FroudeKrylov and diffraction force. The third term is the quadratic drag term of the Morison formula. The fourth term represents linear drag.
Slamming forces
The slamming force on a moving strip, can be related to the change in (added) mass with time.
Expressed in local strip coordinates, the slamming force can be expressed by
where:
\(h\) 
distance between instantaneous surface elevation and strip origin in global zdirection 
1.16.6. Summary of calculation options
Two parameters are used to define different calculation options:

IVOL

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

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


IFOADD

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

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

A summary of the different force components and how they are integrated, is given in Summary of force components.
Force component  IFOADD \(\mathrm {=1}\) 
IFOADD
\(\mathrm {=0}\) 

Inertia: \(\mathrm {F=(m+m_h)\ddot {x}}\) 
\(\mathrm {a}\) to \(\mathrm {z=\zeta}\) 
\(\mathrm {a}\) to \(\mathrm {z=\zeta}\) 
FroudeKrylov: \(\mathrm {F=(\rho V+m_h)\ddot {\zeta}}\) 

longitudinal 
\(\mathrm {a,V}\) to \(\mathrm {z=0}\) 
\(\mathrm {a,V}\) to \(\mathrm {z=0}\) 
transverse 
\(\mathrm {a,V}\) to \(\mathrm {z=\zeta}\) 
\(\mathrm {a,V}\) to \(\mathrm {z=0}\) 
Slamming: \(\displaystyle F=v_r\frac{\mathrm {d}m_h}{\mathrm {d}t}\) 
\(\mathrm {a}\) to \(\mathrm {z=\zeta}\) 
\(\mathrm {a}\) to \(\mathrm {z=\zeta}\) 
Drag: \(\mathrm {F=C_lv_r+C_q \vert v_r\vert v_r}\) 
\(\mathrm {C_1,C_q}\) to \(\mathrm {z=\zeta}\) 
\(\mathrm {C_1,C_q}\) to \(\mathrm {z=0}\) 
FroudeKrylov adjustment: \(\mathrm {F=\rho gV}\) 
\(\mathrm {V;z=0}\) to \(\mathrm {\zeta}\) 
\(\mathrm {V;z=0}\) to \(\mathrm {\zeta}\) 
Buoyancy \(\mathrm {F=\rho gV}\) 
\(\mathrm {V}\) to \(\mathrm {z=0}\) 
\(\mathrm {V}\) to \(\mathrm {z=0}\) 
Buoyancy
from \(\mathrm {z=0}\) to
\(\mathrm {z=\zeta}\) is a correction to the FroudeKrylof
force, for elements crossing the surface.
1.16.7. Wind force on slender element
The aerodynamic force on a slender element is calculated stripwise. 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 Xaxis. 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 coefficients \(\mathrm {C_x^{wi}}\), \(\mathrm {C_y^{wi}}\) and \(\mathrm {C_z^{wi}}\) are distributed coefficients with unit \([\mathrm {kN~\frac{(s/m)^2}{m}}]\) or equivalent. 
The element is allowed to penetrate the water surface. In this case the
aerodynamic force is computed on the dry part of the element. If a slice
is partly above the water and partly below, the computation of relative
wind velocity, force and moment is based on the midpoint 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.
1.17. Specified forces
A number of excitation forces may be specified. Each force is defined by:

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

coordinates

magnitude

attack method (degree of freedom or attack point and direction. The direction may be global or follow the body.)
1.18. Smallbody hydrodynamic forces
The small volume body is described as a point with mass and volume. The motion is described by up to 3 degrees of freedom, all translations.
The hydrodynamic forces on the object are found according to Morison’s formula
where:
\(V'\) 
submerged volume of the body 
\(c_m\) 
added mass coefficient 
\(u_r, \dot{u}_r\) 
relative velocity and acceleration between body and water particles 
\(c_L, c_Q\) 
linear and quadratic drag at velocity = 1 
According to (Lehn and Øritsland, 1987) and (Sarpkaya and Isaacson, 1981), the dynamic equilibrium equation can be written
where:
\(m\) 
mass and hydrodynamic mass of the body 
\(m_h\) 
hydrodynamic mass of the body 
\(z\) 
depth of the object relative to instantaneous sea surface 
\(\rho\) 
density of water 
The dynamic equilibrium equation is written
where:
\('\) 
denotes positiondependent parameters/coefficients 
\(i\) 
denotes degree of freedom (\(i=1,2,3\)) 
\(x,\dot{x},\ddot{x}\) 
body position, velocity and acceleration 
\(m_0\) 
(dry) mass of the body 
\(m_h\) 
hydrodynamic mass of the body 
\(q_i\) 
hydrodynamic force in addition to \(m_{h,i} \ddot{x}_i\) 
\(q_w\) 
partially submerged body weight, \(q_w = m_0g  \rho V'g = m_0g  \rho V_0 c_vg\) 
\(\rho\) 
density of water 
\(q_s\) 
other static or timedependent forces such as coupling forces 
\(V_0\) 
fullysubmerged volume 
\(g\) 
acceleration of gravity 
The hydrodynamic forces, \(\mathrm {q_i}\), are found according to the following equation
where:
\(c_L, c_Q\) 
linear and quadratic drag for unit velocity 
\(u_{r,i}\) 
relative velocity between body and water particles in \(i\)’th direction 
\(\dot{u}_i\) 
wave particle acceleration in \(i\)’th direction 
\(u_r\) 
relative velocity, \(\displaystyle u_r = \sqrt{\sum_{i=1}^3 u_{r,i}^2}\) 
\(z'\) 
vertical displacement between object and sea surface 
\(\dot{x}_{3,W}\) 
velocity of particle on sea surface in vertical direction 
1.19. Soil penetration forces
The soil penetration models apply for both skirted and piled structures. The soil penetration model is related to a smallvolume 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:

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

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

As 1, but the body is allowed to slide horizontally with a friction force in the horizontal direction, depending on the penetrating depth.
For skirted structures, soil resistance to penetration into the seabed is dependent on:

The degree of cyclic motion of the structure.

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

Whether the cells within the skirts are completely enclosed or not in soil around the structure’s perimeter.
For model 1 and 2 the penetrating body is arrested in the horizontal direction at landing. This is accomplished by setting the horizontal velocities, and the accelerating forces in horizontal direction, equal to zero. This implies a simplification, no sliding along the sea bed is included.
The body motion will be finally stopped when either:

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

the body reaches its depth of full penetration
The soil penetration models should be combined with depthdependent hydrodynamic coefficients to account for any proximity effects from the sea bed on added mass and damping properties.
The modelling of soil reaction forces includes:

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

Soil friction and bearing force
For structures with closed compartments, the modelling also includes:

Pressure forces in undisturbed soil, including suction.

Pressure buildup in a fractured compartment caused by hydrodynamic forces related to the underside of the closed compartments. These forces may give additional fractures.
1.19.1. 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 crosssection area of penetrating part of the structure, and \(\mathrm {h}\) is penetration.
1.19.2. 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 buildup starts. The static force due to pressure inside the cavity can be approximated as
where:
cavity
\(\delta\) 
penetration depth, down from the sealing level 
\(k\) 
effective stiffness due to compressibility of water and soil inside the cavity 
\(A\) 
section area of cavity 
The stiffness, \(\mathrm {k}\), can be approximated by
where:
taken as the depth of full penetration)
\(E\) 
elasticity module of water \((=2150~\mathrm{MPa})\) 
\(L\) 
height of cavity (constant value used during the entire penetration, taken as the depth of full penetration) 
When the skirt tip reaches the level of the undisturbed soil, the structure is arrested until the pressure inside the skirt exceeds the soil capacity against failure. This capacity depends on penetration depth.
The water entry and exit of the closed compartment is a dynamic model which sets requirements to the time step. Subdivisions of time steps should possibly be used. The soil penetration model should only be used with the modified Euler integration method.
1.19.3. Pressure buildup in a fractured compartment
When the soil has fractured, the soil force on the skirt will be side friction against motion in both directions and (in case of a concrete structure) also skirt tip pressure, calculated as for the open compartment. When the skirt reaches the bottom of this new fracture, \(\mathrm {h_f}\), the skirt tip is in undisturbed soil and the model above applies.
During this phase, the pressure inside the cavity is calculated from the hydrodynamic forces related to the underside of the closed compartment, and new fractures may occur. These hydrodynamic forces are calculated as:
Drag:
Inertia:
Where \(\mathrm {B_1(z)}\), \(\mathrm {B_2(z)}\) and \(\mathrm {C_a(z)}\) are depthdependent 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.
1.19.4. Soil friction force
The soil reaction force is calculated as:
Downward motion:
Upward motion:
where
resistance, depthdependent
), depthdependent
skirt wall or pile friction
\(F_{d0}\) 
Bearing force for DOWNward motion, typically representing skirt tip resistance, depthdependent 
\(F_{u0}\) 
Bearing force for UPward motion, usually taken as a fraction of \(F_{d0}\), depthdependent 
\(F_1\) 
Friction force for both DOWNward and UPward, typically representing skirt wall or pile friction 
This force is directed opposite to the velocity of the penetrating structure, or if the velocity is zero, opposite to the net force (in this case not exceeding the net force).