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. Time-dependent mass The data group Time-Dependent Mass (TDM) allows the user to model the effects of a mass that can vary both in magnitude and in position. A TDM can be used to model: 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 5 PORTIONS are allowed. Each PORTION has a simple geometry (CONE or BOX). The TDM’s location and orientation relative to the body (vessel) needs to be specified. - in terms of a set of tanks. Each tank are represented by a mesh of triangular panels described in a STL (STereoLithography) file. The maximum number of ballast tanks is unconstrained. 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 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. 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. \[\Delta \boldsymbol{M}(t)=\begin{bmatrix}1&0&0&0&z&-y\\\\0&1&0&-z&0&x\\\\0&0&1&y&-x&0\\\\0&-z&y&y^2+z^2&-xy&-zx\\\\z&0&-x&-xy&z^2+x^2&-yz\\\\-y&x&0&-zx&-yz&x^2+y^2\end{bmatrix}\cdot m(t)=\begin{bmatrix}\boldsymbol{I}&\boldsymbol{A}^T\\\\\boldsymbol{A}&\boldsymbol{A}\boldsymbol{A}^T\end{bmatrix}\cdot m(t)\] where \(\mathrm {\boldsymbol{I}}\) is a 3x3 unit matrix and \(\mathrm {\displaystyle \boldsymbol{A}=\begin{bmatrix}0&-z&y\\\\z&0&-x\\\\-y&x&0\end{bmatrix}}\) The updated mass matrix will then be \(\mathrm {\boldsymbol{M}(t=t_i)=\boldsymbol{M}_0+\Delta \boldsymbol{M}(t=t_i)}\) where \(\mathrm {\boldsymbol{M}_0}\) is the mass matrix originally specified in the system description file. A general expression of the inertia term in the equation of motion is \[\frac{\partial (\boldsymbol{M}\boldsymbol{\dot {x}})}{\partial t}=\boldsymbol{M}\boldsymbol{\ddot {x}}+\boldsymbol{\dot {x}}\frac{\partial \boldsymbol{M}}{\partial t}\] In the present implementation, the last term is not included, because the change in mass does not necessarily give an impulse to the body. 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 x-axis relative to the vessel should be modified. Figure 1. Example of time-dependent mass. 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 Earth-fixed coordinate system \(\mathrm {V_{m0}}\) is the contribution of the panel to the total volume \(\mathrm {V_{m1}}\) is the contribution of the panel to the total volume center \(\mathrm {V_{m2}}\) is the contribution of the panel to the total second volume moment By summing up the contribution of each "contributing panel", the total fluid volume is computed, as well as the position of the center of volume. The two latter quantities are used to compute the total gravity force and moment, using the density of the fluid contained in the tank. The second volume moments are used to build the mass matrix of the tank, which is transposed to the body-fixed coordinate system, using the generalized Huygens theorem. Computing the position of the fluid free surface One challenge consists in finding the position of the fluid free surface in the tank, given the attitude of the tank and a fluid volume. This is done in SIMO by using an iterative process based on an improved version of the secant method. Permeability factor The permeability factor is used to reduce the maximum amount of fluid a tank can contain. The total volume of a tank is computed by SIMO based on the geometry mesh of the tank. The maximum volume of fluid that can be contained in the tank is then equal to the total volume of the tank times the permeability factor. Damaged tank When the state of tank is "damaged", it is assumed that the tank is severely damaged. The fluid density in the tank is replaced by the sea water density. The amount of fluid in the tank is such that the fluid level corresponds to the calm sea water surface level. For example, if the tank is completely submerged, it will be completely filled with sea water. Limitations In the present SIMO version, it is not possible to change the amount of fluid in the a tank during a time domain simulation. The amount of fluid present in the tank is defined at the beginning of the simulation. 1.3. Hydrodynamic reaction forces Based on first order potential wave theory, the hydrodynamic reaction forces are described by added-mass coefficients, \(\mathrm {\boldsymbol{A}}\), and damping coefficients, \(\mathrm {\boldsymbol{C}}\). Both added-mass and damping are frequency-dependent. They are included in the retardation functions, as described in Solution by convolution integral or the first-order motion transfer functions as described in Separation of motions. 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(x-x_0)}}\) where \(\mathrm {\boldsymbol{K}}\) is the hydrostatic stiffness matrix and \(\mathrm {\boldsymbol{x_0}}\) the stiffness reference position. 1.5.2. Non-linear stiffness model When using this formulation, the hydrostatic force model will return the complete buoyancy force (and not the restoring forces as it is the case when using linear stiffness matrix). The buoyancy force is computed based on the part of the hull geometry volume which is under the calm water line (waves are ignored). The geometry of the hull is described by a STL (STereoLithography) file. The mesh of the geometry is then based on triangular panels. At each time step, SIMO will update the position of the body and proceed to a cut of the mesh as follows: if all the three vertices of the panel are above the calm water line, the panel is disregarded if one vertex is above the calm water line, the panel is cut at the calm water line level. The intersection of the calm water line with the edges of the panel defines two vertices which lay on the calm water line. Two new triangular panels are then defined based on the two vertices under the calm water line and the two vertices on the calm water line. They are included in the "contributing panels". if two vertices are above the calm water line, the panel is cut at the calm water line level. The intersection of the calm water line with the edges of the panel defines two vertices which lay on the calm water line. One new triangular panel is then defined based on the vertex under the calm water line and the two vertices on the calm water line. It is included in the "contributing panels". if all the three vertices are under the calm water line, the panel is included in the "contributing panels". Thanks to this cutting procedure, it is not necessary to refine the mesh about the calm water line. The mesh can be coarse as long as it describes the geometry correctly. For each "contributing" panel, the following integrals are computed analytically: \(\mathrm {V_{m0}=\int_{\Gamma}d\Gamma}\) \(\mathrm {V_{m1}=\int_{\Gamma}\boldsymbol{X}d\Gamma}\) where: \(\mathrm {\Gamma}\) is the integration domain defined by the volume "under" the triangular panel \(\mathrm {\boldsymbol{X}}\) is the position vector in the Earth-fixed coordinate system \(\mathrm {V_{m0}}\) is the contribution of the panel to the total volume \(\mathrm {V_{m1}}\) is the contribution of the panel to the total volume center By summing up the contribution of each "contributing panel", the total submerged volume is computed, as well as the position of the center of buoyancy. The two latter quantities are used to compute the total buoyancy forces and moments which are transformed to the body-fixed coordinate system. 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 \[q_j=C_j(\alpha )v^2\] where: \(j\) degree of freedom \(C\) wind force coefficient for the instantaneous relative direction \(v\) relative wind speed between body and wind \(\alpha\) relative velocity direction in local coordinate system The low frequency body motion is used to calculate the relative wind speed. Optionally, the velocity of the body may be omitted when calculating relative velocity. 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 Earth-fixed coordinate system \(\mathrm {\phi }\) : roll angle of the body \(\mathrm {\theta }\) : pitch angle of the body In that case, the coefficients are interpolated for a given body position and a relative wind direction. The force is computed as in equation Equation (3). 1.7. Wave excitation forces Wave forces may be split into three contributions: first-order forces that oscillate with wave frequencies second-order mean, rapidly varying, and slowly varying wave drift forces. higher-order ringing forces. It is assumed that the wave forces can be formally expressed by the following equation. The two first terms are a Volterra series truncated after the quadratic term, see (Hasselmann, 1966) and (Dalzell, 1976), while the last term is the time-dependent ringing load. \[\begin{array}{l}q(t)=q^{(1)}(t)+q^{(2)}(t)+q^{(R)}(t)\\\\\displaystyle =\frac{1}{2\pi }\int^\infty_{-\infty}h^{(1)}(\tau_1)\zeta(t-\tau_1)\mathrm {d}\tau_1+\frac{1}{4\pi ^2}\int^\infty_{-\infty}\int^\infty_{-\infty}h^{(2)}(\tau_1,\tau_2)\zeta(t-\tau_2)\zeta(t-\tau_1)\mathrm {d}\tau_1\mathrm {d}\tau_2+q^{(R)}(t)\end{array}\] where: ringing moment \(q^{(1)}(t)\) time-dependent first-order wave force \(q^{(2)}(t)\) time-dependent second-order wave force \(q^{(R)}(t)\) time-dependent third-order ringing force or third- and fourth-order ringing moment \(h^{(1)}\) linear impulse response function \(h^{(2)}\) second-order impulse response function It is assumed that the kernels \(\mathrm {h^{(1)}}\) and \(\mathrm {h^{(2)}}\) are smooth, absolutely integrable and possess Fourier transforms as follows \[h^{(1)}(\tau)=\frac{1}{2\pi }\int^\infty_{-\infty}H^{(1)}(\omega )e^{i\omega \tau}\mathrm {d}\omega\] \[H^{(1)}(\omega )=\int^\infty_{-\infty}h^{(1)}(\tau)e^{-i\omega \tau}\] \[h^{(2)}(\tau_1,\tau_2)=\int^\infty_{-\infty}\int^\infty_{-\infty}H^{(2)}(\omega _1,\omega _2)e^{i(\omega _1\tau_1+\omega _2\tau_2)}\mathrm {d}\omega _1\mathrm {d}\omega _2\] \[H^{(2)}(\omega _1\omega _2)=\frac{1}{4\pi ^2}\int^\infty_{-\infty}\int^\infty_{-\infty}h^{(2)}(\tau_1,\tau_2)e^{-i(\omega _1\tau_1+\omega _2\tau_2)}\mathrm {d}\tau_1\mathrm {d}\tau_2\] where: \(H^{(1)}\) first-order transfer function \(H^{(2)}\) second-order transfer function The physical interpretation of first-order transfer functions is well established in linear response analysis. The interpretation of the second-order transfer function is not obvious. A simple example taken from (Dalzell, 1976) shows that the second-order response to a dual harmonic excitation with frequencies \(\mathrm {\omega _i}\) and \(\mathrm {\omega _j}\) consists of a mean value in addition to the contribution from sum frequencies \(\mathrm {(\omega _i+\omega _j)}\) and difference frequencies \(\mathrm {(\omega _i-\omega _j)}\). The magnitude of the second-order forces is relatively small compared to first-order forces, but may be important for structures with natural frequencies below the wave period range. A typical example is resonant heave/roll/pitch motions of TLPs. The second-order forces are usually of importance only if any of the difference or sum frequencies coincide with natural frequencies of the system which are outside the frequency range of the first-order forces. The ringing response of TLPs with high natural heave and roll/pitch periods \((3-5~\mathrm {s})\) or large GBSs in deep water with high fundamental modal period \((3.5-6~\mathrm {s})\) has been shown to introduce extreme dynamic loads causing major concerns for the platform designs. From several studies one has concluded that a sufficiently accurate ringing load model is the most critical part of the modelling of the ringing phenomenon. Ringing is related to response in irregular seas. 1.8. First-order wave excitation forces First order wave forces are described in the frequency domain as a transfer function between wave elevation and force \[\boldsymbol{q}_{WA}^{(1)}(\omega )=\boldsymbol{H}^{(1)}(\omega )\zeta(\omega )\] where: \(\boldsymbol{H}^{(1)}\) complex first order transfer function \(\tilde{\zeta} (\omega)\) complex harmonic wave component 1.9. Second-order wave forces The second-order wave force for unidirectional irregular waves can be expressed in the form \[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_nH_{mn}^{(2+)}e^{i(\omega _m+\omega _n)t}\}+\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}e^{i(\omega _m-\omega _n)t}\}\] where: \(\tilde{\zeta}_m\) complex Fourier component of sea surface elevation with frequency \(\omega_m\) and direction \(\beta\) \(H_{mn}^{(2+)} \equiv H^{(2+)} (\omega_m, \omega_n, \beta, \beta)\) second-order transfer function for the sum-frequency force \(H_{mn}^{(2-)} \equiv H^{(2-)} (\omega_m, \omega_n, \beta, \beta)\) second-order transfer function for the difference-frequency 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 (short-crested) 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 second-order wave forces. This method is outlined in (Newman, 1974). (Faltinsen and Løken, 1979) and (Pinkster, 1981) compared Newman’s simplified approach with the different contributions from the second-order transfer function. Neglecting the sum-frequency force in equation Equation (11), the slowly-varying second-order force, also referred to as the wave drift force, is expressed by \[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}e^{i(\omega _m-\omega _n)t}\}\] where the time average is \[\bar{q}=\frac{1}{2}\mathrm {Re}\{\sum_m\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}\}\] The main contribution to the slowly-varying force is assumed to be associated with frequency combinations for which, \(\mathrm {\omega _m\approx \omega _n}\) or \[|\omega _m-\omega _n|<<\frac{1}{2}(\omega _m+\omega _n)\] Assuming that \[H_{mn}^{(2-)}\approx H_{mm}^{(2-)}\] the slowly-varying force can be approximated by \[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{mn}^{(2-)}e^{i(\omega _m-\omega _n)t}\}\] The double sum is circumvented using the relations \[\begin{array}{lll}\displaystyle (\mathrm {Re}\{\sum_m\tilde{\zeta}_m\sqrt{|H_{mn}^{(2-)}|}e^{i\omega _mt}\})^2&=&\displaystyle \frac{1}{2}\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n\sqrt{|H_{mm}^{(2-)}H_{nn}^{(2-)}|}e^{i(\omega _m+\omega _n)t}\}\\\\\displaystyle &+&\displaystyle \frac{1}{2}\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*\sqrt{|H_{mm}^{(2-)}H_{nn}^{(2-)}|}e^{i(\omega _m-\omega _n)t}\}\end{array}\] and \[\sqrt{|H_{mm}^{(2-)}H_{nn}^{(2-)}|}\approx |H_{mm}^{(2-)}|\] Thus, the slowly varying force may be written \[q(t)=\overline{\big(L^+(t)\big)^2}-\overline{\big(L^-(t)\big)^2}\] where \[L^\pm=\mathrm {Re}\{\sum_m\zeta_m\sqrt{\pm2H_{mm}^{(2-)}}e^{i\omega _mt}\}\] The overline symbolizes low-pass filtering. These sums are to be calculated only for those terms where the argument of the square root is positive in each case. (Kaasen, 1999) has shown that a modification to Newman’s method can be done in order to eliminate the high frequency noise generated by this method. By not using pre-generation of force time series, filtering out this noise must be done simultaneously with the simulation where there is a significant challenge to avoid a phase lag that may affect the low frequency components. Hence, avoiding the high frequency noise is a valuable improvement. By using the difference-frequency force transfer functions in both parts of Equation (10), the force is \[q(t)=\mathrm {Re}\{\sum_m\sum_nH_{mn}^{(2-)}[\tilde{\zeta}_m\tilde{\zeta}_n^*e^{i(\omega _m-\omega _n)t}+\tilde{\zeta}_m\tilde{\zeta}_ne^{i(\omega _m+\omega _n)t}]\}\] As in Equation (13) it is assumed that \(\mathrm {H_{mn}^{(2-)}\approx H_{mm}^{(2-)}}\), and two new complex variables are introduced \[\begin{array}{lll}\displaystyle u(t)=\sum_m\sqrt{H_{mm}^{(2-)}}\mathrm {Re}\{\tilde{\zeta}_me^{i\omega _mt}\}&&\displaystyle v(t)=\sum_m\sqrt{H_{mm}^{(2-)}}\mathrm {Im}\{\tilde{\zeta}_me^{i\omega _mt}\}\end{array}\] where \(\mathrm {H_{mm}^{(2-)}}\) can be positive and negative, and \(\mathrm {v(t)}\) is the Hilbert transform of \(\mathrm {u(t)}\). Then a result for the slowly varying force as in Equation (15), can then be seen from \[q(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*H_{nm}^{(2-)}e^{i(\omega _m-\omega _n)t}\}=\mathrm {Re}\{u(t)^2+v(t)^2\}\] It is possible to pre-calculate wave drift forces for a number of headings before time domain simulation starts, and interpolate between these time series in time domain. It is possible to correct the wave drift force coefficients for the presence of wave-current interaction effects using Aranhas method. This is described in more detail in section Wave-current correction of second order wave drift force and wave drift damping using an extended Aranhas formula. 1.10.1. Second order wave forces in short-crested seas The following method is implemented in order to include the effect of short-crested seas on wave drift forces. The average drift force in irregular waves is \[\bar{q}=\int^\infty_02\bar{H}^{(2)}(\omega )S_\zeta(\omega )\mathrm {d}\omega\] in which the average, or uni-directional drift force coefficient for the direction \(\mathrm {\beta _1}\) is \[\bar{H}^{(2)}(\omega )=\int_{\beta _1-\frac{\pi }{2}}^{\beta _1+\frac{\pi }{2}}H^{(2)}(\beta ,\omega )\theta (\beta -\beta _1)\mathrm {d}\beta\] 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. \[\begin{array}{lllll}\tilde{\zeta}_{kl}(\beta _k,\omega _l)&&k=1,2, \cdots ,n_k && l=1,2, \cdots ,n_l\end{array}\] The surface elevation at the origin is represented by the components \(\mathrm {\tilde{\zeta}_l(\omega _l)}\) \[\tilde{\zeta}_l(\omega _l)=\sum_{k=1}^{n_k}\tilde{\zeta}_{kl}\] The drift force due to the component \(\mathrm {\tilde{\zeta}_{kl}}\) is \[\bar{q}_l=\sum_{k=l}^{n_k}|\tilde{\zeta}_{kl}|^2H^{(2)}(\beta _k,\omega _l)\] An average drift force coefficient for the frequency is thus computed as \[\hat{H}_l^{(2)}=\frac{\displaystyle \sum_{k=l}^{n_k}|\tilde{\zeta}_{kl}|^2H^{(2)}(\beta _k,\omega _l)}{\displaystyle \sum_{k=l}^{n_k}|\tilde{\zeta}_{kl}|^2}\] The expected value of this coefficient is equal to \(\mathrm {\bar{H}^{(2)}(\omega _l)}\), but the discrete sampling over directions gives a scatter similar to the scatter of \(\mathrm {|\tilde{\zeta}_l|}\). Representing the sea by its surface elevation and using the average coefficient \(\mathrm {\hat{H}_l^{(2)}}\), the methods of Hsu and Blenkarn and Newman can be used directly for generating time series of slowly varying wave force. The validity for short-crested seas is unknown and should be subjected to further investigation. The validity also becomes questionable when the frequency dependence of \(\mathrm {H^{(2)}(\omega )}\) is strong, i.e. when \(\mathrm {H^{(2)}(\omega )}\) varies rapidly compared to the wave spectrum \(\mathrm {S_\zeta(\omega )}\). 1.10.2. Second order wave forces calculated with FFT and COS method in short-crested seas Long crested waves The calculation is based on the complex wave Fourier components \(\mathrm {X(\omega )}\) and the real drift force coefficients \(\mathrm {C(\omega )}\) for the actual direction. The drift force is calculated from the complex coefficients \[Z(\omega )=X(\omega )\sqrt{C(\omega )}\] FFT and COSINE simulation method will give the same result. FFT method short-crested seas For each frequency, SIMO will calculate \(\mathrm {Z(\omega )}\) from equivalent coefficients \(\mathrm {X(\omega )}\) and \(\mathrm {C(\omega )}\) \[\begin{array}{l}C(\omega )=\frac{\displaystyle \sum_{i=1}^{n_{dir}}X_i(\omega )X_i^*(\omega )c_i(\omega )}{\displaystyle \sum_{i=1}^{n_{dir}}X_i(\omega )X_i^*(\omega )}\\\\\displaystyle X(\omega )=\sum_{i=1}^{n_{dir}}X_i(\omega )\end{array}\] COS-series short-crested seas For each frequency, SIMO will calculate an equivalent coefficient \(\mathrm {Z(\omega )}\) as \[Z(\omega )=\sum_{i=1}^{n_{dir}}X_i(\omega )\sqrt{c_i(\omega )}\] Depending on the degree of spreading, these two methods will give different results. 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 non-dimensional average wave drift forces with respect to relative velocity. This method uses a frequency-independent 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: \[\mathbf{q}(\mathbf{v}_r) \approx \mathbf{q}(\mathbf{v}_{r,0}) + \mathbf{B}(\mathbf{v}_{r} - \mathbf{v}_{r,0})\] 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, \[\begin{bmatrix} F_1 \\ F_2 \\ F_6 \\ \end{bmatrix} = \begin{bmatrix} B_{11} & B_{12} & B_{16} \\ B_{21} & B_{22} & B_{26} \\ B_{61} & B_{62} & B_{66} \\ \end{bmatrix} \begin{bmatrix} v_{1c} - v_{1b} \\ v_{2c} - v_{2b} \\ -v_{6b} \\ \end{bmatrix}\] 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, \[\begin{bmatrix} F_1 \\ F_2 \\ F_6 \\ \end{bmatrix} = - \begin{bmatrix} B_{11} & B_{12} & B_{16} \\ B_{21} & B_{22} & B_{26} \\ B_{61} & B_{62} & B_{66} \\ \end{bmatrix} \begin{bmatrix} v_{1b} \\ v_{2b} \\ v_{6b} \\ \end{bmatrix}\] The coefficients in the wave drift damping term is time varying and can be written as: \[B_{ij}(t)=\mathrm {Re}\{\sum_m\sum_n\tilde{\zeta}_m\tilde{\zeta}_n^*B_{mn,ij}e^{i(\omega _m-\omega _n)t}\}\] 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. 1.11.2. Simplified method The wave drift damping model is based on the derivatives of non-dimensional average wave drift forces with respect to relative velocity. \[c_{wd}=\frac{\partial \bar{F}_i}{\partial v_{ir}\bar{F}_i}\] where \(\mathrm {\bar{F}_i}\) is the average force for a wave of unit wave amplitude. The coefficients \(\mathrm {c_{wd}}\) normally depend on wave frequency, and can be estimated from model tests or 3D sink-source analyses. Assuming that \[\frac{\partial F(t)}{\partial vF(t)}\approx \frac{\partial \overline{F(t)}}{\partial v\overline{F(t)}}\] and \[\frac{\partial \overline{F(\beta )}}{\partial v\overline{F(\beta )}}\approx \frac{\partial \overline{F(\beta +\Delta \beta )}}{\partial v\overline{F(\beta +\Delta \beta )}}\] the wave drift forces in the time domain simulation are modified according to \[\hat{F}_i(t)=F_i(t)(1+c_{wd}v_{wr}), i=1,2,6\] where \(\mathrm {v_{wr}=v_{wc}-v_{wb}}\) is the difference between current and body velocity in the wave direction. \(\mathrm {c_{wd}}\) is implemented in the program as a function of estimated peak wave period, and is fixed during the time domain simulation. This is a simplified model that reflects both wave drift damping and wave drift - current interaction in the time domain. 1.12. Wave-current 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 off-diagonal 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. Wave-current interaction correction of QTFs using extended Aranha method Consider a stationary floating body subjected to a current speed \(U\) with a going-to direction of \(\theta\) relative to the body x-axis and waves propagating in direction \(\beta\) relative to the body x-axis. The current components in x and y direction then becomes: \[\begin{split} u & =U\cos\theta\\ v & =U\sin\theta \end{split}\] The current components inline and transverse to the wave direction becomes: \[\begin{split} u_w & = u\cos\left(\beta_{i}\right)+v\sin\left(\beta_{i}\right) = U\cos\theta_w\\ v_w & = -u\sin\left(\beta_{i}\right)+v\cos\left(\beta_{i}\right) = U\sin\theta_w \end{split}\] where \(\theta_w = \theta - \beta\). The full QTF with wave-current interaction effects are then taken as, \[Q_{U}\left(\omega_{e,1},\omega_{e,2},\beta_{e,1},\beta_{e,2}\right)=B\left(\omega_{1},\omega_{2},\beta_{1},\beta_{2}\right)Q_{0}\left(\omega_{e,1},\omega_{e,2},\beta_{e,1},\beta_{e,2}\right)\] Here, the frequency of encounter can be written as, \[\omega_{e,i}=\left(1+\frac{u_w}{C_{p,i}}\right)\omega_{i}\] and the direction of encounter can be written as: \[\beta_{e,i}=\beta_{i}+\frac{v_w}{C_{g,i}}\] 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: \[B\left(\omega_{1},\omega_{2},\beta_{1},\beta_{2}\right)=\sqrt{A\left(\omega_{1},\beta_{1}\right)A\left(\omega_{2},\beta_{2}\right)}\] The Aranha scaling factor can be written as, \[A\left(\omega,\beta\right)=1-\frac{\left(\omega\frac{\partial\alpha}{\partial\omega}-2\right)}{C_{g}}u_w\] where, \[\alpha=\frac{C_{g}}{C_{p}}\] The dispersion relation in the finite water depth case reads, \[\omega=\sqrt{kg\tanh\left(kh\right)}\] which gives a phase velocity and group velocity: \[C_{p}=\frac{\omega}{k}=\sqrt{\frac{g\tanh\left(kh\right)}{k}}\] \[C_{g}=\frac{d\omega}{dk}=\alpha C_{p}\] \[\alpha =\frac{C_{g}}{C_{p}}=\frac{1}{2}\left(1+\frac{2kh}{\sinh\left(2kh\right)}\right)=\frac{1}{2}+\frac{kh}{\sinh\left(2kh\right)}\] \[\begin{split} \frac{d\alpha}{d\omega} & =\frac{d\alpha}{dk}\frac{dk}{d\omega}\\ & =\frac{1}{C_{g}}\frac{d\alpha}{dk}\\ & =\frac{h}{C_{g}}\frac{\left(1-2kh\coth\left(2hk\right)\right)}{\sinh\left(2kh\right)}\\ & =\frac{h}{C_{g}\sinh\left(2kh\right)}\left(1-\frac{2kh}{\tanh\left(2kh\right)}\right) \end{split}\] 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: \[\mathbf{F}(\omega,\beta,\mathbf{U}_{r})=(1+2\frac{κ}{C_{g}}(u_{r}\cos\left(\beta\right)+v_{r}\sin\left(\beta\right)))\mathbf{F}(\omega_{e},\beta_{r},\mathbf{0})\] where \(\mathbf{U}_{r}=\begin{bmatrix} u_{r} \\ v_{r} \end{bmatrix} = \begin{bmatrix} u_{0}-u_{b} \\ v_{0}-v_{b} \end{bmatrix}\) is the relative velocity between the current and the body low frequency velocity \(\mathbf{U}_{b}=\begin{bmatrix} u_{b}\\ v_{b} \end{bmatrix}\) Using Taylor expansion at first order, about the current velocity \(\mathbf{U}_{r}=\begin{bmatrix} u_{0}\\ v_{0} \end{bmatrix}\) : \[\mathbf{F}(\omega,\beta,\mathbf{U}_{r})=\mathbf{F}(\omega,\beta,\mathbf{U}_{0})-\frac{\partial\mathbf{F}}{\partial\mathbf{U}_{r}}(\omega,\beta,\mathbf{U}_{0})\mathbf{U}_{b}\] The term \(-\frac{\partial\mathbf{F}}{\partial\mathbf{U}_r}(\omega,\beta,\mathbf{U}_{0})\mathbf{U}_{b}\) represents the wave drift damping that can be interpreted as the change of the drift force due to the low frequency velocity of the vessel. It is possible to calculate the coefficient \(\mathbf{B}(\omega,\beta,\mathbf{U}_0)=\frac{\partial\mathbf{F}(\omega,\beta,\mathbf{U}_r)}{\partial\mathbf{U}_r}_{|\mathbf{U}_r=\mathbf{U}_0}\) by using Equation (53) only for surge, sway and yaw (respectively noted as degrees of freedom no.1,2 and 6): \[\begin{equation} \begin{aligned} \mathbf{B}(\omega,\beta,\mathbf{U}_0) & = \begin{bmatrix} B_{11} & B_{12}\\ B_{21} & B_{22}\\ B_{61} & B_{62} \end{bmatrix} \\ & = \frac{\partial\mathbf{F}(\omega,\beta,\mathbf{U}_r)}{\partial\mathbf{U}_r}_{|\mathbf{U}_r=\mathbf{U}_0} \\ & =P(\mathbf{U}_r)\mathbf{F}(\omega_{e},\beta_{r},\mathbf{0})\left[\begin{array}{cc} \frac{\partial}{\partial u_r} & \frac{\partial}{\partial v_r}\end{array}\right]_{\begin{array}{|c} u_r=u_{0}\\ v_r=v_{0} \end{array}} \end{aligned} \end{equation}\] where \(P(\mathbf{U}_r)=1+2\frac{\kappa}{C_{g}}(u_r\cos\left(\beta\right)+v_r\sin\left(\beta\right))\) After some calculation, we obtain: \[\frac{\partial\mathbf{F}(\omega,\beta,u_r,v_r)}{\partial u_r}_{\begin{array}{|c} u_r=u_{0}\\ v_r=v_{0} \end{array}}=W_{c0}+P(\mathbf{U}_0)W_{u_r}\] and \[\frac{\partial\mathbf{F}(\omega,\beta,u_r,v_r)}{\partial v_r}_{\begin{array}{|c} u_r=u_{0}\\ v_r=v_{0} \end{array}}=W_{s0}+P(\mathbf{U}_0)W_{v_r}\] where: \[W_{c0}=2\frac{\kappa}{C_{g}}\cos\left(\beta\right)\mathbf{F}(\omega_{e_{0}},\beta_{r_{0}},0,0)\] \[W_{s0}=2\frac{\kappa}{C_{g}}\sin\left(\beta\right)\mathbf{F}(\omega_{e_{0}},\beta_{r_{0}},0,0)\] \[P(\mathbf{U}_0)=1+2\frac{\kappa}{C_{g}}(u_{0}\cos\left(\beta\right)+v_{0}\sin\left(\beta\right))\] \[W_{u_r}=\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\omega}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}k\cos\left(\beta\right)-\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\beta}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}\frac{1}{C_{g}}\sin\left(\beta\right)\] \[W_{v_r}=\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\omega}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}k\sin\left(\beta\right)+\frac{\partial\mathbf{F}(\omega,\beta,0,0)}{\partial\beta}_{\begin{array}{|c} \omega=\omega_{e_{0}}\\ \beta=\beta_{r_{0}} \end{array}}\frac{1}{C_{g}}\cos\left(\beta\right)\] where \(\omega_{e_{0}}\) and \(\beta_{r_{0}}\) are respectively the encounter frequency and the encounter heading for the body at rest (i.e. for \(\mathbf{U}_b=\mathbf{0}\)). The wave drift damping coefficients \(\mathbf{B}(\omega,\beta,\mathbf{U}_0)\) can thus be computed from the following inputs: wave drift force coefficients at zero current velocity as function of wave direction and wave frequency (note: when using the QTF as input, the wave drift force coefficients are extracted from the diagonal) wave directions wave frequencies water depth current velocity current direction The obtained wave drift damping coefficients are: 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 wave-induced motion is based on some simplifying assumptions, in particular linearization of the hydrostatic, diffraction and radiation forces about the still water plane and the body’s hydrostatic equilibrium. To partly overcome this deficiency, SIMO offers a way to account for the nonlinearity in the restoring force and the wave force in a quasistatic manner during the simulation. The basic idea is simple and straightforward: The body surface is divided into 4-node panel cells. At each time step, the location of each cell centroid is calculated in earth-fixed coordinates. The wave elevation at this cell centroid is also evaluated. If the cell centroid is above the mean waterline when the body is in its equilibrium position and below the wave surface at any time step during the simulation, the 3D local nonlinear correction in the restoring and wave forces is formed as the product of the panel area and the quasistatic pressure applied on the panel cell (mass density of the water \(\mathrm {\times }\) acceleration of gravity \(\mathrm {\times }\) vertical distance between the wave surface and the cell centroid). The local nonlinear correction is also included if the cell centroid is below the mean waterline when the body is in the equilibrium position, but above the wave surface at any time step during the simulation. The local nonlinear corrections are integrated to give the total nonlinear correction in the 6 degree-of-freedom generalized rigid body forces and moments in the body-fixed coordinate system. Nonlinear correction of wave force and restoring force on body. Quasistatic pressure is computed over the wetted panels, and the resulting time-varying forces and moments are added as corrections to those obtained using frequency-dependent first-order excitation functions and wave drift functions. The panel model need not cover the entire wetted hull, but only the part that will be partly wet or dry, see Nonlinear correction of wave force…. The user specifies a minimum draft and a maximum draft for the panel model, which is read into SIMO from a text file of type GDF, which is used by WAMIT. Using few panels will reduce the computational time. The wave elevation outside the hull can be that of the incoming, undisturbed wave or that corresponding to the diffracted wave based on first-order theory, if available. However, the above-mentioned nonlinear modification has been partially accounted for in the WAMIT second order theory. This part of WAMIT second order force has to be removed in order to avoid double counting. The total second order hydrodynamic force from WAMIT can be written as \[\boldsymbol{F}^{(2)}=\iint_{S_0}p^{(1)}\boldsymbol{n}^{(1)}\mathrm {d}S+\iint_{S_0}p^{(2)}\boldsymbol{n}^{(0)}\mathrm {d}S+\iint_{S_0}p^{(0)}\boldsymbol{n}^{(2)}\mathrm {d}S+\iint_{S-S_0}p^{(1)}\boldsymbol{n}^{(0)}\mathrm {d}S\] \[\boldsymbol{M}^{(2)}=\iint_{S_0}p^{(1)}\boldsymbol{r}\times \boldsymbol{n}^{(1)}\mathrm {d}S+\iint_{S_0}p^{(2)}\boldsymbol{r}\times \boldsymbol{n}^{(0)}\mathrm {d}S+\iint_{S_0}p^{(0)}\boldsymbol{r}\times \boldsymbol{n}^{(2)}\mathrm {d}S+\iint_{S-S_0}p^{(1)}\boldsymbol{r}\times \boldsymbol{n}^{(0)}\mathrm {d}S\] where \(\mathrm {p^{(0)}}\), \(\mathrm {p^{(1)}}\) and \(\mathrm {p^{(2)}}\) are the hydrostatic water pressure, first order and second order hydrodynamic water pressures while \(\mathrm {\boldsymbol{n}^{(0)}}\), \(\mathrm {\boldsymbol{n}^{(1)}}\) and \(\mathrm {\boldsymbol{n}^{(2)}}\) represent the mean, first order and second order normal vectors of the body surface expressed in the earth-fixed coordinate system. \(\mathrm {\boldsymbol{r}}\) is the position vector. \(\mathrm {S_0}\) denotes the mean wetted surface. \(\mathrm {S-S_0}\) is the difference between the instantaneous wetted surface and the mean wetted surface. The last terms in Equation (63) and Equation (64) can be expressed as \[\iint_{S-S_0}p^{(1)}\boldsymbol{n}^{(0)}\mathrm {d}S=\frac{1}{2}\rho g\int_{WL}\boldsymbol{n}^{(0)}[\zeta^{(1)}_r]^2\{1-[n_z^{(0)}]^2\}^{-1/2}\mathrm {d}l\] \[\iint_{S-S_0}p^{(1)}\boldsymbol{r}\times \boldsymbol{n}^{(0)}\mathrm {d}S=\frac{1}{2}\rho g\int_{WL}\boldsymbol{r}\times \boldsymbol{n}^{(0)}[\zeta^{(1)}_r]^2\{1-[n_z^{(0)}]^2\}^{-1/2}\mathrm {d}l\] where \(\mathrm {\rho }\) and \(\mathrm {g}\) are the mass density of water and the acceleration of gravity, respectively. \(\mathrm {WL}\) stands for waterline. \(\mathrm {\zeta_r^{(1)}}\) is the first order relative wave elevation. \(\mathrm {n_z^{(0)}}\) represents the vertical component of the mean normal vector \(\mathrm {\boldsymbol{n}^{(0)}}\). The removal of Equation (65) and Equation (66) is taken care of in the program. If second-order terms are not included in the simulation these terms are not removed. In the SIMO input manual, the nonlinear correction described above is referred to as "Nonlinear Buoyancy Correction" 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 z-velocity may be given. If these TRF’s are included in a simulation, time series of diffracted waves may be generated, but the TRF’s will not affect the simulation results. If a TRF of velocity is specified, time series of both velocity and acceleration will be generated. A body may be in the "wave shadow" of another, larger body. This effect can be included in SIMO by data group "Diffracted wave application" The body will then "feel" the diffracted waves, instead of the original "undisturbed" wave field. There are two ways to specify that a body shall "feel" diffracted waves. 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 (3-dof body), the effect of wave diffraction may be included in data group SMALL BODY HYDRODYNAMIC DATA by referring to a body and a "diffracted wave point". Diffracted waves are not included in other force models. This means that in order to include diffraction effects in these models, e.g. "Second order wave drift coefficients" the coefficients have to be modified. 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: \[q_{CU}^{(k)}(\alpha ,t)=C_1^{(k)}(\alpha )|u(t)|+C_2^{(k)}(\alpha )|u(t)|^2\] \[|u|^2=(v_1-\dot {x}_1)^2+(v_2-\dot {x}_2)^2\] \[\alpha =\arctan\frac{v_2-\dot {x}_2}{v_1-\dot {x}_1}\] where: \(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 Figure 2. Separation of cross-flow due to translation and rotation. Usually, the linear coefficient, \(\mathrm {C_1}\) is not used. The model does not account for the effect from the yaw-induced cross-flow. This cross-flow has normally been included as a separate yaw damping term. The cross-flow based on a pure relative translational velocity and yaw velocity is illustrated in Separation of cross-flow…. Assuming that a quadratic yaw damping term has been included, the viscous yaw moment may be calculated as the sum of the contribution from the translational velocity and the yaw velocity \[q^{(6)}(\alpha )=C_2^{(6)}(\alpha )|u|^2+b_{66}\dot {x}_6|\dot {x}_6|\] where: \(C_2^{(6)}(\alpha)\) Current coefficient in yaw as a function of the incident heading, \(\alpha\) \(b_{66}\) Quadratic damping coefficient in yaw due to a pure yaw velocity \(\dot{x}_6\) Yaw velocity Normally, the current coefficient in yaw, \(\mathrm {C_2^{(6)}(\alpha )}\), also includes the Munk moment which is an effect that can be derived from potential theory. When quadratic current coefficients are included in the simulation model, the Munk moment is excluded from the equations of motion. The Munk moment, denoted \(\mathrm {q_M(\alpha )}\), is given as \[q_M(\alpha )=-\frac{1}{2}|u|^2[A_2-A_1]\sin(2\alpha )\] where: \(A_2\) Added mass in sway \(A_1\) Added mass in surge It is, however, quite clear that the current velocity, the vessel translational velocity and the yaw velocity will influence the distribution of the cross-flow along the hull. Theoretically, the sum of these contributions will result in a linear cross-flow distribution along the hull. Due to the quadratic nature of the viscous hull forces, the forces obtained from the vessel translation and the current and the forces obtained from a yaw-induced cross-flow can’t be separated and then added. Therefore a force model integrating up the sway force and yaw moment along the vessel due to the resulting instantaneous cross-flow has been established as strip model no 1 and strip model no. 2. In order to get simple, analytical equations of the viscous sway force and yaw moment, the drag properties are assumed to be constant along the vessel. For strip model no. 1: \[C_d(x)=\frac{1}{2}\rho c_d(x)T(x)=\frac{C_2^{(2)}(90\deg)}{l}\] For strip model no. 2: \[C_d(x)=\frac{1}{2}\rho c_d(x)T(x)=\frac{C_2^{(2)}(\alpha )}{l\sin^2(\alpha )}\] where: \(c_d(x)\) Non-dimensional drag coefficient of the vessel section \(C_d(x)\) Drag coefficient of the vessel section \(T(x)\) Vessel draught at the vessel section \(l\) Vessel length \(\rho\) Water density The drag force of a strip, \(\mathrm {f_s(x)}\), located at a distance from the centre of motion is then given as \[f_s(x)=C_d(x)(u_{r2}-x\dot {x}_6)^2\] where \[u_{r2}=|\boldsymbol{u}|\sin(\alpha )\] The total sway force, \(\mathrm {q^{(2)}(\alpha )}\), is then given as \[q^{(2)}(\alpha )=\int_{-l_2}^{l_1}C_d(x)(u_{r2}-x\dot {x}_6)^2\mathrm {d}x\] while the yaw moment is given as \[q^{(6)}(\alpha )=q_v^{(6)}(\alpha )+q_M(\alpha )=\int_{-l_2}^{l_1}C_d(x)(u_{r2}-x\dot {x}_6)^2x\mathrm {d}x+q_M(\alpha )\] Generally, the two alternative models give increased damping moments and sway forces compared to the separated motion model. The two-dimensional current coefficient as applied for strip model no.1 may deviate from the current coefficient for the actual heading of the relative velocity, leading to inaccurate quasistatic current sway force. Strip model no. 2, however, gives correct quasistatic current force, but may give numerical problems for small values of \(\mathrm {\sin(\alpha )}\). When using the strip models in the program, the quadratic current coefficients in yaw are not used, and the Munk moment will be included in the equations of motion. 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 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. 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 framework-type seabed structure and four mudmat structures, one in each corner, for analysis of dynamics in the splash zone or during seabed penetration. The framework can be modelled as a 6-dof body and the mudmats as 3-dof bodies with position-dependent coefficients. imaginary bodies attached to a floating vessel, in order to include non-linear deviations from force coefficients calculated by linear wave diffraction theory programs. As an alternative model, a small 3-dof body can be rigidly connected to a large body by a set of springs having sufficiently high stiffness to avoid excitation of internal resonance. However, this will need extremely short time steps in the simulation. 1.16.2. Element definition Both element types may have any orientation. The main body must be of BODY TYPE 1, i.e. 6-dof body. Small-body theory is used to calculate forces both for the slender element and the concentrated fixed element. Stiff connection implies that all the "sub-body" (slender element or fixed, concentrated element) forces are calculated and transferred to the "main" body. Each slender element specified by the user is divided into NSTRIP strips with equal length. 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 right-handed. For XG and XB, see General Assumptions and Notation. Figure 3. Definition of local strip coordinate system. Figure 3 illustrates a single strip, introducing the local strip coordinate system. The XS(1)-axis is directed along the longitudinal axis of the strip, and hence also the element, the XS(2)-axis is placed perpendicular to the XS(1)-axis in the local strip xy-plane and the XS(3)-axis is then defined from the definition of orthogonal right-handed coordinate systems. 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 \[-\Delta \boldsymbol{M}(t) = \begin{bmatrix} 1&0&0&0&z&-y\\ 0&1&0&-z&0&x\\ 0&0&1&y&-x&0\\ 0&-z&y&y^2+z^2&-xy&-zx\\ z&0&-x&-xy&z^2+x^2&-yz\\ -y&x&0&-zx&-yz&x^2+y^2 \end{bmatrix} \cdot mL = \begin{bmatrix} \boldsymbol{I}&\boldsymbol{A}^T\\ \boldsymbol{A}&\boldsymbol{A}\boldsymbol{A}^T \end{bmatrix}\cdot mL\] where \(\mathrm {\boldsymbol{I}}\) is a 3x3 unit matrix and \[\boldsymbol{A}=\begin{bmatrix}0&-z&y\\\\z&0&-x\\\\-y&x&0\end{bmatrix}\] The mandatory specified BODY MASS DATA for the main body can be zero (but not necessarily), and will be added to the contributions from the slender elements. Contribution from added mass The contribution from the added mass will be an extension to the above, since the added mass of each strip is different in the 3 directions. Assume the following parameters: \(\boldsymbol{V}\) Linear acceleration vector (main body) \(\boldsymbol{W}\) Rotation acceleration vector (main body) \(\boldsymbol{R}\) Radius vector to the strip \(\boldsymbol{M}_a\) Added mass matrix for the strip \[\begin{array}{lllllll}\boldsymbol{V}=\begin{bmatrix}e\\f\\g\end{bmatrix}&&\boldsymbol{W}=\begin{bmatrix}p\\q\\r\end{bmatrix}&&\boldsymbol{R}=\begin{bmatrix}x\\y\\z\end{bmatrix}&&\boldsymbol{M}_a=\begin{bmatrix}a_x&0&0\\0&a_y&0\\0&0&a_z\end{bmatrix}\end{array}\] The local acceleration of the strip will be \[\boldsymbol{V}_L=\boldsymbol{V}+\boldsymbol{W}\times \boldsymbol{R}\] The local force and moment can be expressed by \[\begin{array}{l}\boldsymbol{F}_L=\boldsymbol{M}_a\cdot \boldsymbol{V}+\boldsymbol{M}_a\cdot (\boldsymbol{W}\times \boldsymbol{R})\\\\\boldsymbol{M}_L=\boldsymbol{R}\times \boldsymbol{F}_L=\boldsymbol{R}\times (\boldsymbol{M}_a\cdot \boldsymbol{V})+\boldsymbol{R}\times \big(\boldsymbol{M}_a\cdot (\boldsymbol{W}\times \boldsymbol{R})\big)\end{array}\] Further expansion yields \[\begin{array}{l}\boldsymbol{A}\cdot \boldsymbol{V}=\boldsymbol{R}\times (\boldsymbol{M}_a\cdot \boldsymbol{V})=\begin{bmatrix}-a_yzf+a_zyg\\a_xze-a_zxg\\-a_xye+a_yxf\end{bmatrix}=\begin{bmatrix}0&-a_yz&a_zy\\a_xz&0&-a_zx\\-a_xy&a_yx&0\end{bmatrix}\cdot \begin{bmatrix}e\\f\\g\end{bmatrix}\\\\\boldsymbol{B}\cdot \boldsymbol{W}=\boldsymbol{R}\times \big(\boldsymbol{M}_a\cdot (\boldsymbol{W}\times \boldsymbol{R})\big)=\begin{bmatrix}a_zy^2+a_yz^2&-a_zxy&-a_yxz\\-a_zxy&a_xz^2+a_zx^2&-a_xyz\\-a_yxz&-a_xyz&a_yx^2+a_xy^2\end{bmatrix}\cdot \begin{bmatrix}p\\q\\r\end{bmatrix}\end{array}\] The resulting added mass contribution is \[\Delta \boldsymbol{M}_a=\begin{bmatrix}\boldsymbol{M}_a&\boldsymbol{A}^T\\\boldsymbol{A}&\boldsymbol{B}\end{bmatrix}\] The contribution from all strips will be added to the body mass matrix. Depth-dependent hydrodynamic coefficient corrections may be specified for both slender and fixed elements. Added mass data for the main body may be zero (but not necessarily), and will be added to the contributions from the slender elements. 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 z-direction and is written \[\boldsymbol{F}_{\mathrm {grav},G}=\begin{bmatrix}0\\0\\-mg\mathrm {d}S\end{bmatrix}\] Buoyancy forces The buoyancy force acts in the global z-direction, through the center of buoyancy. The load vector for a strip, expressed in the global coordinate system, can be written, \[\boldsymbol{F}_{\mathrm {buoy},G}=\begin{bmatrix}0\\0\\\rho V_sg\mathrm {d}S\end{bmatrix}\] Wave forces The wave force on a strip acts through the center of buoyancy. The load vector for a strip, expressed in the local strip coordinate system can be written \[\begin{equation} \boldsymbol{F}_{W,S} = \rho V_S\boldsymbol{a}_S + \boldsymbol{m}_h \circ \boldsymbol{a}_S + \boldsymbol{F}_{q,S} + \boldsymbol{C}_l \circ \boldsymbol{U}_{rel,S} \end{equation}\] where the quadratic drag force is computed as, \[\begin{equation} \boldsymbol{F}_{q,S} = \left[\begin{array}{c} C_{qx}U_{rel,sx}\left|U_{rel,sx}\right|\\ C_{qy}U_{rel,sy}\sqrt{U_{rel,sy}^{2}+U_{rel,sz}^{2}}\\ C_{qz}U_{rel,sz}\sqrt{U_{rel,sy}^{2}+U_{rel,sz}^{2}} \end{array}\right] \end{equation}\] and, \[\begin{equation} \boldsymbol{U}_{rel,S} = \boldsymbol{U}_{S}+\boldsymbol{v}_{S}-\boldsymbol{\dot{X}}_{S} \end{equation}\] and we have that: \(V_S\) is submerged volume per length unit, calculated to \(z = 0\) \(\boldsymbol{a}_S =\begin{bmatrix} a_{sx} & a_{sy} & a_{sz} \end{bmatrix}^T\) is water particle acceleration components in local strip coordinate system \(\boldsymbol{U}_{rel,S} = \begin{bmatrix} U_{rel,sx} & U_{rel,sy} & U_{rel,sz} \end{bmatrix}^T\) is the relative water particle velocity components in local strip coordinate system (relative to strip velocity) \(\boldsymbol{v}_S =\begin{bmatrix} v_{sx} & v_{sy} & v_{sz} \end{bmatrix}^T\) is wave induced water particle velocity components in local strip coordinate system \(\boldsymbol{\dot{X}}_S =\begin{bmatrix} \dot{X}_S & \dot{Y}_S & \dot{Z}_S \end{bmatrix}^T\) is strip velocity components in local strip coordinate system \(\boldsymbol{U}_S =\begin{bmatrix} U_{sx} & U_{sy} & U_{sz} \end{bmatrix}^T\) is current flow velocity in local strip coordinate system and \(\qquad \qquad \quad \circ\) denotes elementwise multiplication The first and second terms in equation Equation (87) contains the Froude-Krylov and diffraction force. The third term is the quadratic drag term of the Morison formula. The fourth term represents linear drag. Slamming forces The slamming force on a moving strip, can be related to the change in (added) mass with time. Expressed in local strip coordinates, the slamming force can be expressed by \[\boldsymbol{F}_{S,S} =-\frac{\partial \boldsymbol{m}_s}{\partial t}\boldsymbol{U}_{rel,S} =-\frac{\partial \boldsymbol{m}_s}{\partial h}\frac{\partial h}{\partial t}\boldsymbol{U}_{rel,S} =-\frac{\partial }{\partial h}\begin{bmatrix}m_{h,x(h)}&0&0\\0&m_{h,y(h)}&0\\0&0&m_{h,z(h)}\end{bmatrix}\frac{\partial h}{\partial t}\boldsymbol{U}_{rel,S}\] where: \(h\) distance between instantaneous surface elevation and strip origin in global z-direction 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. Table 1. Summary of force components Force component IFOADD \(\mathrm {=1}\) IFOADD \(\mathrm {=0}\) Inertia: \(\mathrm {F=(m+m_h)\ddot {x}}\) \(\mathrm {a}\) to \(\mathrm {z=\zeta}\) \(\mathrm {a}\) to \(\mathrm {z=\zeta}\) Froude-Krylov: \(\mathrm {F=(\rho V+m_h)\ddot {\zeta}}\) longitudinal \(\mathrm {a,V}\) to \(\mathrm {z=0}\) \(\mathrm {a,V}\) to \(\mathrm {z=0}\) transverse \(\mathrm {a,V}\) to \(\mathrm {z=\zeta}\) \(\mathrm {a,V}\) to \(\mathrm {z=0}\) Slamming: \(\displaystyle F=v_r\frac{\mathrm {d}m_h}{\mathrm {d}t}\) \(\mathrm {a}\) to \(\mathrm {z=\zeta}\) \(\mathrm {a}\) to \(\mathrm {z=\zeta}\) Drag: \(\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}\) Froude-Krylov adjustment: \(\mathrm {F=\rho gV}\) \(\mathrm {V;z=0}\) to \(\mathrm {\zeta}\) \(\mathrm {V;z=0}\) to \(\mathrm {\zeta}\) Buoyancy \(\mathrm {F=\rho gV}\) \(\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 Froude-Krylof force, for elements crossing the surface. 1.16.7. Wind force on slender element The aerodynamic force on a slender element is calculated strip-wise. The element is divided into a number of slices, as indicated in Wind force on slender element…. The fineness of the division is chosen by the user. It is assumed that the cross sectional geometry is constant and does not vary along the X-axis. Cartesian axes XYZ are defined as shown in Wind force on slender element…. Let wind load coefficients \(\mathrm {C_x^{wi}}\), \(\mathrm {C_y^{wi}}\) and \(\mathrm {C_z^{wi}}\) be associated with the three axial directions. Let \(\mathrm {v_x^{wi}}\), \(\mathrm {v_y^{wi}}\) and \(\mathrm {v_z^{wi}}\) be the components of relative wind velocity at the centre of one slice. Here relative wind velocity means the hypothetical wind speed at the location of the slice centre (i.e. the wind velocity that would be if there was no slender element) minus the velocity of the slice centre. The three force components on the slice are then modelled as \[\begin{array}{l}\Delta F_x^{wi}=C_x^{wi}\:|v_x^{wi}|\:v_x^{wi}\Delta x\\\\\Delta F_y^{wi}=C_{xy}^{wi}\:v_{yz}^{wi}\:v_y^{wi}\Delta x\\\\\Delta F_z^{wi}=C_z^{wi}\:v_{yz}^{wi}\:v_z^{wi}\Delta x\end{array}\] where \(\mathrm {\Delta x}\) is the length of the slice and \(\mathrm {v_{yz}^{wi}}\) is the component of the relative wind velocity in the YZ plane. \[v_{yz}^{wi}=\sqrt{(v_y^{wi})^2+(v_z^{wi})^2}\] The forces from all the strips are added to give the total force on the slender element. The corresponding components of force moment are obtained by multiplication by the respective moment arms from the defined point of reference, which may not be located as shown in Wind force on slender element…. Wind force on slender element. Decomposition of relative wind velocity vector shown. the coefficients \(\mathrm {C_x^{wi}}\), \(\mathrm {C_y^{wi}}\) and \(\mathrm {C_z^{wi}}\) are distributed coefficients with unit \([\mathrm {kN~\frac{(s/m)^2}{m}}]\) or equivalent. The element is allowed to penetrate the water surface. In this case the aerodynamic force is computed on the dry part of the element. If a slice is partly above the water and partly below, the computation of relative wind velocity, force and moment is based on the mid-point of the part of the slice that is above the water surface. As for hydrodynamic force, the calculation will use the still water level or the instantaneous surface level, dependent on the parameter IFOADD, see Body Data Specification in the SIMO User Guide. 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. Small-body hydrodynamic forces The small volume body is described as a point with mass and volume. The motion is described by up to 3 degrees of freedom, all translations. The hydrodynamic forces on the object are found according to Morison’s formula \[q_h=\rho V'c_m\dot {u}_r+u_r\big(c_L+|u_r|c_Q\big)\] where: \(V'\) submerged volume of the body \(c_m\) added mass coefficient \(u_r, \dot{u}_r\) relative velocity and acceleration between body and water particles \(c_L, c_Q\) linear and quadratic drag at velocity = 1 According to (Lehn and Øritsland, 1987) and (Sarpkaya and Isaacson, 1981), the dynamic equilibrium equation can be written \[q=\frac{\partial (m\dot {x})}{\partial t}=\frac{\partial m}{\partial t}\dot {x}+m\frac{\partial \dot {x}}{\partial t}\] \[\frac{\partial m}{\partial t}\dot {x}=\frac{\partial m}{\partial z}\frac{\partial z}{\partial t}\dot {x}=\frac{\partial m_h}{\partial z}\frac{\partial z}{\partial t}\dot {x}\] \[m\frac{\partial \dot {x}}{\partial t}=m\ddot {x}\] where: \(m\) mass and hydrodynamic mass of the body \(m_h\) hydrodynamic mass of the body \(z\) depth of the object relative to instantaneous sea surface \(\rho\) density of water The dynamic equilibrium equation is written \[(m_0+m_{h,i})\ddot {x}_i=q_i+q_{w,i}+q_{s,i}\] where: \('\) denotes position-dependent parameters/coefficients \(i\) denotes degree of freedom (\(i=1,2,3\)) \(x,\dot{x},\ddot{x}\) body position, velocity and acceleration \(m_0\) (dry) mass of the body \(m_h\) hydrodynamic mass of the body \(q_i\) hydrodynamic force in addition to \(m_{h,i} \ddot{x}_i\) \(q_w\) partially submerged body weight, \(q_w = m_0g - \rho V'g = m_0g - \rho V_0 c_vg\) \(\rho\) density of water \(q_s\) other static or time-dependent forces such as coupling forces \(V_0\) fully-submerged volume \(g\) acceleration of gravity The hydrodynamic forces, \(\mathrm {q_i}\), are found according to the following equation \[q_i=\big(\rho V'+m_{h,i}\big)\dot {u}_i+\big(c_L+|u_r|\:c_Q\big)u_{r,i}+\frac{\partial m_{h,i}}{\partial z'}\big(\dot {x}_3-\dot {x}_{3,W}\big)u_{r,i}\] where: \(c_L, c_Q\) linear and quadratic drag for unit velocity \(u_{r,i}\) relative velocity between body and water particles in \(i\)’th direction \(\dot{u}_i\) wave particle acceleration in \(i\)’th direction \(u_r\) relative velocity, \(\displaystyle u_r = \sqrt{\sum_{i=1}^3 u_{r,i}^2}\) \(z'\) vertical displacement between object and sea surface \(\dot{x}_{3,W}\) velocity of particle on sea surface in vertical direction 1.19. Soil penetration forces The soil penetration models apply for both skirted and piled structures. The soil penetration model is related to a small-volume body, which means that all rotations are excluded for the model. There are three conditions which are considered in the modelling of soil resistance to seabed penetration. These are: The structure causes no pressure build-up such that the soil resistance against vertical movement is from vertical friction forces and resistance on the skirt/pile tip area. 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 depth-dependent hydrodynamic coefficients to account for any proximity effects from the sea bed on added mass and damping properties. The modelling of soil reaction forces includes: Soil buoyancy force (if the soil can be regarded as a fluid) Soil friction and bearing force For structures with closed compartments, the modelling also includes: Pressure forces in undisturbed soil, including suction. Pressure buildup in a fractured compartment caused by hydrodynamic forces related to the underside of the closed compartments. These forces may give additional fractures. 1.19.1. Soil buoyancy force The buoyancy due to higher density of the soil compared to water is included in the model as \[F_b=(\rho _s-\rho _w)A_shg\] where \(\mathrm {\rho _s}\) is mass density of soil, \(\mathrm {\rho _w}\) is density of water, \(\mathrm {A_s}\) is cross-section area of penetrating part of the structure, and \(\mathrm {h}\) is penetration. 1.19.2. Pressure buildup in a closed compartment The volume flow is calculated as \[q=-c\sqrt{p_c}\] Where \(\mathrm {c}\) is the flow coefficient and \(\mathrm {p_c}\) is the net pressure inside the cavity, including any suction. \(\mathrm {V}\) is the total volume flow since the cavity became sealed and is calculated as \[V=\int_{t_s}^tq\mathrm {d}t\] where \(\mathrm {t_s}\) is the time since sealing occurred. The compressibility of water (and soil) is considered. Once the skirt reaches the bottom of a new failure (or the sea bed, if no failure occurs at landing) the closed compartment is sealed, and the pressure build-up starts. The static force due to pressure inside the cavity can be approximated as \[F=k(\frac{V}{A}-\delta )\] where: cavity \(\delta\) penetration depth, down from the sealing level \(k\) effective stiffness due to compressibility of water and soil inside the cavity \(A\) section area of cavity The stiffness, \(\mathrm {k}\), can be approximated by \[k=\frac{EA}{L}\] where: taken as the depth of full penetration) \(E\) elasticity module of water \((=2150~\mathrm{MPa})\) \(L\) height of cavity (constant value used during the entire penetration, taken as the depth of full penetration) When the skirt tip reaches the level of the undisturbed soil, the structure is arrested until the pressure inside the skirt exceeds the soil capacity against failure. This capacity depends on penetration depth. The water entry and exit of the closed compartment is a dynamic model which sets requirements to the time step. Subdivisions of time steps should possibly be used. The soil penetration model should only be used with the modified Euler integration method. 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: \[\big(B_2(z)-B_{20}\big)\:|\dot {z}|\dot {z}+\big(B_1(z)-B_{10}\big)\dot {x}\] Inertia: \[\rho V\big(C_a(z)-C_{a0}\big)\ddot {z}\] Where \(\mathrm {B_1(z)}\), \(\mathrm {B_2(z)}\) and \(\mathrm {C_a(z)}\) are depth-dependent hydrodynamic coefficients and \(\mathrm {B_{10}}\), \(\mathrm {B_{20}}\) and \(\mathrm {C_{a0}}\) are values valid for no sea bed interference. Each single fracture will open up a ventilation channel along a part of the skirt periphery. The width and depth of the channels at fracture must be assessed by the user, and the input parameters (soil forces and hydrodynamic forces) given accordingly. Depending on the above penetration forces and the hydrodynamic parameters valid for ventilation underneath the skirt, the damping of motion can increase sharply with depth. 1.19.4. Soil friction force The soil reaction force is calculated as: Downward motion: \[F_d=F_{d0}(z)+F_1(z)\] Upward motion: \[F_u=F_{u0}(z)-F_1(z)\] where resistance, depth-dependent ), depth-dependent skirt wall or pile friction \(F_{d0}\) Bearing force for DOWNward motion, typically representing skirt tip resistance, depth-dependent \(F_{u0}\) Bearing force for UPward motion, usually taken as a fraction of \(F_{d0}\), depth-dependent \(F_1\) Friction force for both DOWNward and UPward, typically representing skirt wall or pile friction This force is directed opposite to the velocity of the penetrating structure, or if the velocity is zero, opposite to the net force (in this case not exceeding the net force). Force models Hydrodynamic interaction between bodies