This module only:

Hydrodynamic Load on Partly Submerged, Floating Elements

This section will give an introduction to the hydrodynamic load models for floating partly submerged elements, which are relevant for analysis of floating flexible fish farm systems, floating offshore loading hoses, etc. The presentation is to a large extent based on Ormberg (1991).

A floating body of small volume may, during a wave cycle, accomplish large motions from its mean position. Thus, the hydrodynamic loads should be computed at the instantaneous position rather than at the mean position, which means that the assumptions within linear potential theory are violated. The applied hydrodynamic load model, described by Aarsnes et al. (1988), is an engineering approach to a solution of the problem. It involves hydrodynamic coefficients deduced by linear theory. In addition, an adequate description of the wave kinematics in the wave zone is required. Application of the model including comparison with model test results is described by Ormberg (1991).

The hydrodynamic forces are calculated based on two-dimensional strip theory, which requires the structural members to be long and slender. Hydrodynamic interaction between structural parts is neglected. The wave-induced excitation forces (Froude-Krylov and diffraction forces) are computed by a long wavelength approximation which involves added mass and potential damping of the actual cross section together with the wave kinematics. The viscous loads are computed using the drag term in Morison’s equation.

In the following, a rigid structural member (section) with constant cross sectional dimensions, is considered. The section with local (co-rotating) coordinate system \(\mathrm {(x,y,z)}\) , is illustrated in Figure 1. The figure also defines the space-fixed (global) coordinate system with X- and Y-axes in the still water plane and with the Z-axis pointing upwards.

Figure 1. Structural member divided into sub-elements, for hydrodynamic load calculation

As indicated in Figure 1 , the section is divided into sub-elements for hydrodynamic load calculation. The forces acting on a sub-element are calculated based on the instantaneous orientation and position relative to the free surface. The water particle kinematics, which are assumed to be constant within the sub-element, are calculated at the buoyancy center of the sub-element. Correspondingly, the load per unit length is constant over the length of the sub-element and attacks through the buoyancy center. The total force normal to the section is found by simply adding the force contribution on each sub-element.

For the computer program implementation, a table of the two-dimensional added mass and damping coefficients for different levels of immersion at the actual wave frequency is pre-computed by use of the Frank close fit method, Frank (1967).

In the following, the hydrodynamic forces are understood as the total forces acting on the member from the surrounding fluid. This implies that buoyancy and current forces are also included in the hydrodynamic forces. Thus, the total hydrodynamic forces are given as

\[\boldsymbol{F^H}=\boldsymbol{F^{\mathrm {Pot}}}+\boldsymbol{F^D}=\boldsymbol{F^{FK}}+\boldsymbol{F^S}+\boldsymbol{F^R}+\boldsymbol{F^D}\]

where:

\(\mathrm {\boldsymbol{F^{Pot}}}\) : the potential flow contribution to hydrodynamic forces

\(\mathrm {\boldsymbol{F^{FK}}}\) : Froude-Krylov forces including buoyancy forces

\(\mathrm {\boldsymbol{F^S}}\) : diffraction forces

\(\mathrm {\boldsymbol{F^R}}\) : added mass and damping forces

\(\mathrm {\boldsymbol{F^D}}\) : the viscous forces or drag forces

1. Potential Flow Contribution

1.1. Loads in longitudinal direction

The long and slender shape assumption implies that the forces on the end surfaces are dominated by the Froude-Krylov contribution. Hence, other contributions are neglected. The importance of including the pressure force on end plates is thoroughly discussed by Hooft (1972).

The end pressure forces are calculated based on dynamic Froude-Krylov pressure and the hydrostatic pressure

\[p_{\textrm{tot}} = p_{\textrm{dyn}} + p_{\textrm{st}} = p_{\textrm{dyn}} - \rho gZ\]

where \(\mathrm {g}\) is the acceleration due to gravity and \(\mathrm {\rho }\) is the density of water.

Thus, the total potential contribution to load in local x-direction of the structural member is found by the vector sum of the end pressure forces, which formally is written

\[\boldsymbol{F_x^{FK}}=(\int_{A_w^a}\!{p_{\mathrm {tot}}^{(a)}}\textrm{d}{A}-\int_{A_w^b}\!{p_{\mathrm {tot}}^{(b)}}\textrm{d}{A})\boldsymbol{i_I}\]

where \(\mathrm {A_w}\) is the instantaneous immersed area of the end surface. The indexes (a) and (b) indicate the end surface number, see Figure 2 For practical purposes it is assumed that the following simplified formula is sufficient

\[\boldsymbol{F_x^{FK}}=(p_{\textrm{tot}}^{(a)}A_w^{(a)}-p_{\textrm{tot}}^{(b)}A_w^{(b)})\boldsymbol{i_I}\]

where \(\mathrm {p_\textrm{tot}}\) is calculated at the wetted area center.

Figure 2. End pressure forces

1.2. Froude-Krylov (FK) forces:

The Froude-Krylov forces including the buoyancy forces acting normal to a sub-element may according to Kaplan and Hue (1959) and Dixon et al. (1979), be computed by

\[\Delta \boldsymbol{F_n^{FK}}=\boldsymbol{f_n^{FK}}\Delta x=(\boldsymbol{f_n^{\mathrm {Buoy}}}+\rho A_S\boldsymbol{\dot u_n})\Delta x\]

where \(\mathrm {\boldsymbol{f_n^{\textrm{Buoy}}}}\) is the buoyancy force calculated based on the instantaneous immersed cross section area. \(\mathrm {\Delta x}\) is the length of the sub-element. \(\mathrm {\boldsymbol{\dot u_n}}\) is the water particle acceleration vector component normal to the sub-element, i.e.

\[\boldsymbol{\dot u_n}=\dot u_y\boldsymbol{i_2}+\dot u_z\boldsymbol{i_3}\]

The buoyancy force referred to the global system is simply calculated by

\[\boldsymbol{\tilde{f}^{\mathrm {\,Buoy}}}=\rho gA_S\boldsymbol{I_3}\]

where \(\mathrm {A_S}\) is the immersed cross section area of the sub-element.

Transforming the buoyancy load to local system, neglecting the axial contribution to avoid double representation of this term, the FK-force referred to local system is written

\[\Delta \boldsymbol{F_n^{FK}}=\boldsymbol{f_n^{FK}}\Delta x=(\rho gA_ST_{23}+\rho A_S\dot u_y)\Delta x\boldsymbol{i_2}+(\rho gA_ST_{33}+\rho A_S\dot u_z)\Delta x\boldsymbol{i_3}\]

where \(\mathrm {\boldsymbol{T_{ij}}}\) is the transformation matrix between the global and local coordinate systems, i.e.

\[\boldsymbol{i_i}=\boldsymbol{T_{ij}I_j}\]

1.3. Diffraction forces:

According to the long wave assumption, the diffraction forces may be expressed in terms of two-dimensional added mass and damping

\[\Delta \boldsymbol{F_n^S}=\boldsymbol{f_n^S}\Delta x=(A_{22}^{(2D)}\boldsymbol{\dot u_y}+B_{22}^{(2D)}\boldsymbol{u_y})\Delta x\boldsymbol{i_2}+(A_{33}^{(2D)}\boldsymbol{\dot u_z}+B_{33}^{(2D)}\boldsymbol{u_z})\Delta x\boldsymbol{i_3}\]

where \(A^{(2D)}=A^{(2D)}(Z_\mathrm {rel})\) and \(B^{(2D)}=B^{(2D)}(Z_{\mathrm {rel}})\) are position-dependent added mass and damping per unit length, respectively. \(Z_\mathrm {rel}\) is the vertical position relative the free surface.

1.4. Added mass and damping forces:

The added mass and damping forces including moment around local x-axis are calculated according to

\[\begin{bmatrix}\Delta F_y^R\\\Delta F_z^R\\\Delta F_{\theta x}^R\end{bmatrix}=(\begin{bmatrix}A_{22}^{(2D)}&0&A_{24}^{(2D)}\\0&A_{33}^{(2D)}&0\\A_{42}^{(2D)}&0&A_{44}^{(2D)}\end{bmatrix}\begin{bmatrix}\ddot v_y\\\ddot v_z\\\ddot v_{\theta x}\end{bmatrix}+\begin{bmatrix}B_{22}^{(2D)}&0&B_{24}^{(2D)}\\0&B_{33}^{(2D)}&0\\B_{42}^{(2D)}&0&B_{44}^{(2D)}\end{bmatrix}\begin{bmatrix}\dot v_y\\\dot v_z\\\dot v_{\theta x}\end{bmatrix})\]

where \(\mathrm {\boldsymbol{\dot v}}\) and \(\mathrm {\boldsymbol{\ddot v}}\) are the structural velocity and acceleration, respectively.

2. Simplified Approach

A further simplification of the hydrodynamic problem is included by introducing constant added mass coefficients, i.e. neglecting the depth dependency. This means that the added mass is taken as proportional to the submerged cross sectional area. Accordingly (Eq. 7.99) and (Eq. (7.100) simplify to

\[\Delta \boldsymbol{F_n^S}=\boldsymbol{f_n^S}\Delta x=\rho A_SC_{my}\dot u_y\Delta x\boldsymbol{i_2}+\rho A_SC_{mz}\dot u_z\Delta x\boldsymbol{i_3}\]

and

\[\Delta \boldsymbol{F_n^R}=\boldsymbol{f_n^R}\Delta x=-\rho A_SC_{my}\ddot v_y\Delta x\boldsymbol{i_2}-\rho A_SC_{mz}\ddot v_z\Delta x\boldsymbol{i_3}\]

where \(\mathrm {A_S}\) is submerged cross section area. \(\mathrm {C_{my}}\) and \(\mathrm {C_{mz}}\) are constant added mass coefficients referred to completely submerged cross section.

3. Viscous Forces

The viscous forces are computed using the drag force term of Morison’s equation. The cross flow principle first introduced by Hoerner (1965), is used to determine the forces in transverse direction of the sub-elements. The drag coefficients are assumed to be constant and given for a submerged cross section.

The viscous forces or drag forces are calculated based on relative velocity (\(\mathrm {\boldsymbol{u_r}}\)) in the local system expressed by

\[\boldsymbol{u_r}=\boldsymbol{u_c}+\boldsymbol{u_w}-\boldsymbol{\dot v}=u_{rx}\boldsymbol{i_1}+u_{ry}\boldsymbol{i_2}+u_{rz}\boldsymbol{i_3}\]

where \(\mathrm {u_c}\) is the current velocity, \(\mathrm {u_w}\) is the wave particle velocity and \(\mathrm {\dot v}\) is the structural velocity.

3.1. Loads in longitudinal direction

A friction force contribution is included in the axial direction

\[\Delta \boldsymbol{F^D_x}=\boldsymbol{f_x^D}\Delta x=\frac{1}{2}\rho C_{Dx}L_W|u_{rx}|u_{rx}\Delta x\boldsymbol{i_1}\]

where \(\mathrm {C_{Dx}}\) is the skin friction coefficient and \(\mathrm {L_W}\) is the instantaneous wetted part of the cross section circumference.

3.2. Transverse loads

From Equation (14), the transverse relative velocity vector is written

\[\boldsymbol{u_{rn}}=u_{ry}\boldsymbol{i_2}+u_{rz}\boldsymbol{i_3}\]

The transverse drag force is calculated according to

\[\Delta \boldsymbol{F_n^D}=\boldsymbol{f_n^D}\Delta x=\frac{1}{2}\rho C_{Dy}h_{\mathrm {rel}}|\boldsymbol{u_{rn}}|u_{ry}\Delta x\boldsymbol{i_2}+\frac{1}{2}\rho C_{Dz}b_\mathrm {rel}|\boldsymbol{u}_{rn}|u_{rz}\Delta x\boldsymbol{i_3}\]

where

\(\begin{array}{l}\displaystyle h_\mathrm {rel}=\frac{A_S}{A}h\\\\\displaystyle b_\mathrm {rel}\,=\frac{A_S}{A}b\end{array}\)

\(\mathrm {C_{Dy}}\) and \(\mathrm {C_{Dz}}\) are drag coefficients in local y and z directions respectively, \(\mathrm {A_S}\) is the instantaneous submerged cross section area and \(\mathrm {A}\) is the cross section area. \(\mathrm {b}\) and \(\mathrm {h}\) are characteristic width (y-direction) and height (z-direction) of the section respectively.

4. Moment around Local x-axis

Assuming the transverse forces attack at the center of buoyancy, see Figure 3, the moment around the local x-axis due to wave excitation forces and drag forces is approximated by

\[\Delta \boldsymbol{F_{\theta x}^h}=\boldsymbol{f_{\theta x}^h}\Delta x=(\boldsymbol{r_B}\times \boldsymbol{f_n^h})\Delta x\]

where \(\mathrm {\boldsymbol{r_B}=y_B\boldsymbol{i_2}+z_B\boldsymbol{i_3}}\) is the position vector of the buoyancy center and the transverse hydrodynamic excitation forces are given by

\[\boldsymbol{f_n^h}=\boldsymbol{f_n^{FK}}+\boldsymbol{f_n^S}+\boldsymbol{f_n^D}\]

Figure 3. Excitation forces acting at the buoyancy center

5. Discussion

Computation of added mass and potential damping according to the Frank close fit method requires that the structural member is horizontal, i.e. the local xy-plane is parallel with the still water plane, and that the cross section is symmetric about its local x-z plane. During time domain simulation, the added mass and damping are taken as direction-dependent coefficients connected to local (structural) axis directions, while strictly they are valid for global directions. Hence, introduced as a practical limit, the instantaneous pitch and roll angles have to be within \(\mathrm {\pm30^{\circ}}\) inclination to the surface.

The loads due to slamming or water impact are not included. These may be included as the time rate of change of momentum associated with the immersed portion of the structure. Calculating the vertical impact forces on horizontal members for the water entry case is well documented by e.g. Greenhow (1987), Kaplan and Silbert (1976), Faltinsen (1977). However, in this more general case regarding structure orientation, water exit as well as entry, further investigation is needed. Hence, this term is neglected.

The long wave length assumption is used in the calculation of potential flow acting on the sub-element, which normally requires the wave length / diameter ratio to be greater than 5. It is further assumed that the local angular displacements are not too large relative to the free surface. This assumption is due to that the computation of hydrodynamic coefficients are precalculated for a horizontal structural member.

Experiments indicate that the drag coefficients for a cross section in the free surface zone are strongly dependent on the submergence, see e.g. Øritsland (1986). These effects are not included. The dependency on Reynolds number, Keulegan-Carpenter number and roughness should be taken into account when specifying the constant drag coefficients as input.

6. Implementation in Finite Element Formulation

Different formulations might be used for the numerical implementation of hydrodynamic loads. These consider the order of the polynomial used for the displacement functions and the representation of the distribution of the hydrodynamic element loads. In the current version of the program, a consistent formulation is adopted for the numerical implementation of hydrodynamic loads, i.e. using the same displacement functions as for the element deformation description. The contributions from added mass and damping forces are included in the inertia and damping load vectors of the system, while the other hydrodynamical loads contribute to the external load vector.

Analogous to hydrodynamic strip theory, the load distribution is taken care of by dividing the element into a specified number of sub-elements for load calculation, see figure Figure 4. Within a sub-element, the hydrodynamic loads are assumed to be constant. The corresponding hydrodynamic mass- and damping matrices and nodal load vector are computed by numerical integration simply by adding up the contribution from each sub-element. This is certainly not the most accurate method for numerical integration, but is found to be adequate considering the uncertainties in the hydrodynamic load computation.

Based on this formulation, the element length may be chosen with regard to structural deformations of the actual problem, while the lengths of the sub-elements are chosen with regard to hydrodynamic load distribution independently of the element length.

The hydrodynamic loads may in general be described as non-conservative loads i.e. dependent on position, direction and/or loaded area. By linearizing the incremental dynamic equation, tangential mass, damping and stiffness matrices are introduced, see chapters Finite Element Formulation and Dynamic Time Domain Analysis. However, the linearization also introduces an additional stiffness due to the displacement dependent loads. This stiffness is called load correction stiffness, are described in section Incremental Equilibrium Iterations. The inclusion of the load correction stiffness only has a consequence for the rate of convergence in a numerical solution procedure. Provided uniqueness of the solution, it is the equilibrium equation alone that governs the final solution. A thorough discussion regarding type of loads and the consequence of implementing the corresponding stiffness matrix is presented by Mathisen (1990).

In the current version of the program, a load correction stiffness matrix associated with the water plane stiffness (buoyancy springs) is included.

In the following, the formulation of hydrodynamic loads are described as developed for beam elements. This formulation applies for the bar elements by use of corresponding displacement functions and by omitting rotational degrees of freedom.

Figure 4. Element divided into sub-elements for hydrodynamic load calculation

6.1. Hydrodynamic mass and damping matrices

The hydrodynamic (added) mass and damping matrices are established based on the hydrodynamic mass and damping calculated along the secant length of the deformed element.

The components of these matrices are computed according to

\[\begin{array}{l}\begin{aligned}m_{vv}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_v^T}(x_\textrm{i})\Delta m_{vv}^h(x_i)\boldsymbol{N_v}(x_i)\\m_{ww}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_w^T}(x_\textrm{i})\Delta m_{ww}^h(x_i)\boldsymbol{N_w}(x_i)\\m_{{\theta _x}{\theta _x}}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_{\theta _x}^T}(x_\textrm{i})\Delta m_{{\theta _x}{\theta _x}}^h(x_i)\boldsymbol{N_{\theta _x}}(x_i)\\m_{v{\theta _x}}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_v^T}(x_\textrm{i})\Delta m_{v{\theta _x}}^h(x_i)\boldsymbol{N_{\theta _x}}(x_i)\\\\\end{aligned}\end{array}\]

\[\begin{array}{l}\begin{aligned}c_{vv}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_v^T}(x_\textrm{i})\Delta b_{vv}^h(x_i)\boldsymbol{N_v}(x_i)\\c_{ww}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_w^T}(x_\textrm{i})\Delta b_{ww}^h(x_i)\boldsymbol{N_w}(x_i)\\c_{{\theta _x}{\theta _x}}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_{\theta _x}^T}(x_\textrm{i})\Delta b_{{\theta _x}{\theta _x}}^h(x_i)\boldsymbol{N_{\theta _x}}(x_i)\\c_{v{\theta _x}}^h&=\sum_{\textrm{i=1}}^{\textrm{nsub}}\,\boldsymbol{N_v^T}(x_\textrm{i})\Delta b_{v{\theta _x}}^h(x_i)\boldsymbol{N_{\theta _x}}(x_i)\\\end{aligned}\end{array}\]

where \(\mathrm {x_i}\) is taken at the buoyancy center of the sub-element. \(\mathrm {\Delta m_{\alpha \gamma }^h}\) and \(\mathrm {\Delta b_{\alpha \gamma }}\) are added mass and damping for the sub-element. These terms are calculated according to the applied hydrodynamic load model.

If the added mass and damping forces are established according to Equation (8) the added mass for the sub-element is computed by

\[\begin{array}{l}\begin{aligned}&\Delta m_{vv}^h&=A_{22}^{(2D)}\Delta x_0\\&\Delta m_{ww}^h&=A_{33}^{(2D)}\Delta x_0\\&\Delta m_{\theta _x\theta _x}^h&=A_{44}^{(2D)}\Delta x_0\\&\Delta m_{v\theta _x}^h&=A_{24}^{(2D)}\Delta x_0\\\end{aligned}\end{array}\]

and the damping computed by

\[\begin{array}{l}\begin{aligned}&\Delta b_{vv}^h&=B_{22}^{(2D)}\Delta x_0\\&\Delta b_{ww}^h&=B_{33}^{(2D)}\Delta x_0\\&\Delta b_{\theta _x\theta _x}^h&=B_{44}^{(2D)}\Delta x_0\\&\Delta b_{v\theta _x}^h&=B_{24}^{(2D)}\Delta x_0\\\end{aligned}\end{array}\]

where \(\mathrm {A_{ij}^{(2D)}}\) and \(\mathrm {B_{ij}^{(2D)}}\) are 2-dimensional added mass and damping per unit length dependent on the instantaneous immersion. \(\mathrm {\Delta x_0}\) is the length of the sub-element in the initial condition.

If the simplified approach is applied for the hydrodynamic load calculations, Equation (13), no hydrodynamic damping is considered. The added mass for the sub-element is established according to

\[\begin{array}{l}\begin{aligned}&\Delta m_{vv}^h&=\rho A_SC_{my}\Delta x_0\\&\Delta m_{ww}^h&=\rho A_SC_{mz}\Delta x_0\\\end{aligned}\end{array}\]

where \(\mathrm {C_{my}}\) and \(\mathrm {C_{mz}}\) are directional dependent added mass coefficients, and \(\mathrm {A_S}\) is submerged cross section area.

6.2. External hydrodynamic load vector

The hydrodynamic loads are computed along the secant length of the deformed element.

Considering the transverse load contributions the hydrodynamic load on each sub-element is written

\[\boldsymbol{f_n^h}=\boldsymbol{f_n^{FK}}+\boldsymbol{f_n^S}+\boldsymbol{f_n^D}\]

where

\(\mathrm {\boldsymbol{f_n^{FK}}}\) : Froude-Krylov force including buoyancy, computed according to Equation (8).

\(\mathrm {\boldsymbol{f_n^S}}\) : diffraction force computed according to (Eq. 7.99) or Equation (12).

\(\mathrm {\boldsymbol{f_n^D}}\) : drag force computed according to Equation (17).

The external moment around local x-axis is approximated according to

\[\boldsymbol{f_{\theta _x}^h}=\boldsymbol{r^B}\times \boldsymbol{f_n^h}\]

where \(\mathrm {\boldsymbol{r^B}}\) is the vector from the principal axis to the buoyancy center, see Figure 4. In the axial direction, \(\mathrm {F_x^{FK}}\) and \(\mathrm {f_x^D}\) are computed according to Equation (4) and Equation (15), respectively.

The element end forces, \(\mathrm {F_x^{FK}}\) , Froude-Krylov forces including hydrostatic pressure, are assumed to be constant over the element length. Hence, the cross sectional axial force (stress resultant) is expressed in terms of effective axial force, see section Overview of Load Effects. The effect axial force is used directly in the expression for the geometric stiffness.

Based on these considerations, the components of the external hydrodynamic load vector is computed according to

\[\begin{array}{l}\displaystyle \boldsymbol{S_u^h}=\frac{1}{2}F_x^{FK}\begin{bmatrix}1\\1\end{bmatrix}\sum_{\mathrm {i}=1}^{\textrm{nsub}}\,\boldsymbol{N_u^T}(x_\mathrm {i})f_x^D(x_\mathrm {i})\Delta x_0\\\\\displaystyle \boldsymbol{S_v^h}=\sum_{\mathrm {i}=1}^{\textrm{nsub}}\,\boldsymbol{N_v^T}(x_\mathrm {i})f_y^h(x_\mathrm {i})\Delta x_0\\\\\displaystyle \boldsymbol{S_w^h}=\sum_{\mathrm {i}=1}^{\textrm{nsub}}\,\boldsymbol{N_w^T}(x_\mathrm {i})f_z^h(x_\mathrm {i})\Delta x_0\\\\\displaystyle \boldsymbol{S_{\theta _x}^h}=\sum_{\mathrm {i}=1}^{\textrm{nsub}}\,\boldsymbol{N_{\theta _x}^T}(x_\mathrm {i})f_{\theta _x}^h(x_\mathrm {i})\Delta x_0\\\\\displaystyle \boldsymbol{S_{\theta _y}^h}=0\\\\\displaystyle \boldsymbol{S_{\theta _z}^h}=0\end{array}\]

where \(\mathrm {N}\) denotes the actual displacement function. \(\mathrm {N}\) and \(\mathrm {f}\) are taken at the buoyancy center \(\mathrm {x_i}\) of each sub-element. \(\mathrm {\Delta x_0}\) is the initial length of the sub-element.

6.3. Hydrodynamic stiffness matrix

The stiffness associated with water plane stiffness (buoyancy springs) is included in the stiffness matrix. The stiffness terms are limited to contributions due to translational displacement in global Z-direction and torsional rotation around local x-axis. For the computation of the hydrodynamic stiffness matrix, linear displacement functions are applied.

Thus, the element hydrodynamic stiffness matrix with regard to Z-displacement, \(\mathrm {\boldsymbol{\tilde{k}_{ww}^{E}}}\) , is calculated directly in global system according to

\[\boldsymbol{\tilde{k}_{ww}^{E}}=\sum_{\mathrm {i}=1}^{\textrm{nsub}}\,\boldsymbol{N^T}(x_\mathrm {i})\frac{\partial \,\tilde{f}_Z^{\mathrm {\,Buoy}}(x_\mathrm {i})}{\partial Z}\boldsymbol{N}(x_\mathrm {i})\Delta x_0\]

where \(\mathrm {\frac{\partial \,\tilde{f}_Z^{~\textrm{Buoy}}}{\partial Z}}\) is the water plane stiffness of the sub-element based on the water plane area in the global XY-plane

The local torsional hydrodynamic stiffness matrix is calculated according to

\[\boldsymbol{k_{\theta _x\theta _x}^E}=\sum_{\mathrm {i}=1}^{\textrm{nsub}}\,\boldsymbol{N^T}(x_\mathrm {i})(\frac{\partial f^\textrm{Buoy}_{\theta _x}(x_\mathrm {i})}{\partial \theta _x})\boldsymbol{N}(x_\mathrm {i})\Delta x_0\]

Note that this formulation results in a symmetric hydrodynamic stiffness matrix.