Contact Formulation 1. General The following features are modelled by use of contact formulation: seafloor contact soil contact roller contact tensioner contact The seafloor contact formulation includes contact normal to the seafloor modelled by springs and in-plane contact modelled by a combination of springs and frictional forces. The soil contact formulation enables modelling monopile foundations of marine structures by springs and in-plane contact modelled by springs connected to beam element nodes. The roller and tensioner contact formulations enable modelling of pipelines during laying operations. This includes contact forces between the pipe and rollers on the lay barge/stinger and applied tension from a tensioner (Holthe et al. 1995, Ormberg et al. 1995). 2. Seafloor Contact The seafloor may either be horizontal or be defined using depths at each point in a regular 2D grid. In a nonlinear static or time domain analysis, the stiffness matrix will normally be re-established at each step and iteration. Nodes may therefore obtain or lose seafloor contact, or may slide along the seafloor. The interaction forces are introduced into the unbalanced force vector during equilibrium iterations. In a linear analysis in either time or frequency domain, the stiffness matrix is established at the final static equilibrium position and used thereafter. In-plane spring-friction is linearized to the spring stiffness, regardless of whether sliding occurred in the static analysis. No equilibrium corrections or iterations are made during the subsequent dynamic analysis. Interaction forces are therefore not calculated. Effects of nodes coming in contact with or lifting off the seabed or starting/stopping sliding are therefore not included. One consequence is that a node that did not have contact at the end of the static analysis, will not come into contact with the seabed - even if it is lowered below the seafloor. The seafloor contact is handled separately in the direction normal to the seafloor and in the seafloor plane. In the seafloor plane, the contact is currently handled separately along the projection of the corresponding beam or bar element on the seafloor and transverse to the projection. These three directions are respectively denoted normal, axial and lateral. 2.1. Contact normal to the seafloor Contact normal to the seafloor may either be modelled with bi-linear springs or with the nonlinear consolidated riser-soil model. 2.1.1. Bi-linear spring contact The riser-seabed interaction is described by a stiffness parameter giving the interaction force per indentation depth per length of riser. For each element having seafloor contact at one or two nodes, the corresponding spring stiffnesses are added to the system stiffness matrix. The value of the spring stiffness is calculated from the normal stiffness parameter and the unstressed element length. Interaction forces will be calculated from the spring values and the indentation depths. Indentation depths are found from the water depth at the location, \(\mathrm {D}\), the external contact radius if applicable and the nodal displacements. For an element lying on the seabed - at depth \(\mathrm {D}\), the spring forces are in equilibrium with the element’s submerged weight. The indentation may then be found as \[\delta =\delta _0=\frac{w_s}{k}\] where \(\mathrm {w_s}\) is the submerged weight per length riser \([\mathrm {F/L}]\) and \(\mathrm {k}\) is the vertical stiffness parameter \([\mathrm {F/L^2}]\). In general, the indentation depth is found as \[\delta =-(D+Z-R_e)\] where \(\mathrm {Z}\) is the node’s vertical coordinate - negative below the mean water level and \(\mathrm {R_e}\) is the external contact radius, if applicable. All nodes with \[Z-R_e\leq-D\] will be in contact with the sea floor. Note that the external contact radius, \(\mathrm {R_e}\) , is set to zero in the original seafloor contact formulation. The nodal interaction force is subsequently found as \[F_n=\frac{1}{2}kL_e\delta\] where \(\mathrm {L_e}\) is the element length. The nodal interaction force is added to the internal force vector of the system. The corresponding nodal stiffness contribution to the system stiffness matrix is calculated according to \[k_n=\frac{1}{2}kL_e\] Damping may be handled similarly. 2.1.2. Consolidated riser-soil model The Consolidated riser-soil model is described in detail in (Voie et al. 2014). The contact model includes nonlinear stiffness as the pipe is pressed into the seafloor and a separate unloading curve as the pipe is lifted. Suction forces are activated as the pipe is lifted off the seafloor. At a new contact cycle, the soil deformations from the previous cycle are retained so that contact will occur at a slightly greater depth. The memory of seafloor contact follows each seafloor contact element, so that a point that once has been in contact with the seafloor will thereafter not have contact with undisturbed seafloor – even if the new contact is at some distance from the old contact. The model is thus best suited to cases where a given node with contact remains within a limited area. The main parameters to the model are geotechnical soil parameters; soil submerged weight, soil shear strength at seafloor, gradient of soil shear strength with depth, Poisson ratio and shear modulus. A single node seafloor contact element is added at each end of each beam or bar element with potential seafloor contact. The contact element’s properties will depend in the seafloor contact component and on the orientation and contact radius of the element that it is attached to. The projection of the beam / bar element on the seafloor defines the local axial direction, which together with the normal vector to the seafloor determines the in-plane lateral direction. A riser-soil contact element becomes active when its external contact radius is at or below the seafloor. The contact element’s external contact radius is set to the external contact radius of the beam / bar element that it is attached to. This external contact radius must have a positive value. The riser-soil contact normal to the seafloor is not influenced by the in-plane spring-friction contact. This means that the soil capacity is not reduced due to non-normal contact. The in-plane friction forces will, however, depend on the normal contact force. In the present version, the riser-soil contact element memory is activated together with the activation of the in-plane seafloor contact; i.e. the static load "FRIC". 2.2. In-plane contact In-plane contact with the seafloor is modelled as a combination of initial springs and friction forces due to sliding. Behaviour in axial and lateral directions is currently assumed to be completely independent of each other. Input to the model consists of four parameters: \(\mathrm {k_a\quad }\) : Axial spring stiffness parameter \([\mathrm {F/L^2}]\) \(\mathrm {k_l\quad }\) : Lateral spring stiffness parameter \([\mathrm {F/L^2}]\) \(\mathrm {\mu_a\quad }\) : Axial friction coefficient \(\mathrm {[-]}\) \(\mathrm {\mu_l\quad }\) : Lateral friction coefficient \(\mathrm {[-]}\) Initial springs are introduced when a node is in contact with the seafloor. Springs are removed if contact is lost and reinstated if new contact arises. The in-plane seafloor contact is in general modelled as shown in Figure 1. Figure 1. In-plane seafloor contact. For the in-plane contact, the reference distance, \(\mathrm {d}\), is the sum of incremental axial/lateral displacements since the corresponding initial spring for this node was introduced. The interaction force is calculated from the reference distance and the spring value when sliding has not occurred. If sliding has occurred, the interaction force is calculated from the normal contact force \(\mathrm {F_n}\) (denoted \(\mathrm {F_z}\) inthe figure) and the friction coefficient. The incremental stiffnesses are defined as \[\begin{array}{l}\frac{1}{2}k_iL_e\quad \mathrm {for}\quad \frac{1}{2}k_iL_e\delta _i\leq\mu_iF_n\\\\0.0\quad \mathrm {for}\quad \frac{1}{2}k_iL_e\delta _i>\mu_iF_n\end{array}\] and the interaction force, \(\mathrm {F_i}\), as \[\begin{array}{l}F_i=\frac{1}{2}k_iL_e\delta _i\quad \mathrm {for}\quad \frac{1}{2}k_iL_e\delta _i\leq\mu_iF_n\\\\F_i=\mu_iF_n\quad \mathrm {for}\quad \frac{1}{2}k_iL_e\delta _i>\mu_iF_n\end{array}\] where the index \(\mathrm {i}\) denotes either axial or lateral direction. Note that sliding can start due to two different changes from one step to another; either an increase of the reference distance or a decrease of the normal contact force. Optionally, the lateral friction may be applied at the lower edge of the element, thus introducing a torsional moment. 3. Soil Contact Contact between soil and embedded structures is defined through a soil layer profile. The profile defines the sub-surface below a mudline as a vertical slice of soil layers where each layer is defined by relevant properties such as specific weight, strength and stiffness, together with specific soil characteristics. The soil characteristics are usually established empirically and may be re-used at multiple sites and depths. The two main types of soil are sand and clay. The embedded structure gains added stiffness from soil contact with respect to lateral displacement. Depending on the modelling methodology there may also be added stiffness with respect to vertical displacements or shear deformation (cross-section rotation). The added stiffness is modelled through nonlinear-elastic springs acting on element nodes. The material law of the springs are known as soil reaction curves, the shape of which are defined by the soil layer profile properties, soil characteristics and pile geometry properties such as embedded length, depth below the mudline and pile diameter. The currently implemented soil contact is based on PISA methodology which is intended for modelling of monopile foundations with large diameter-to-embedded-length ratio. For details on the PISA methodology see, e.g., (Byrne et al. 2017, Byrne et al. 2019 and Burd et al. 2020). 4. Roller Contact The main intention of the roller contact formulation is to enable the user to model the contact between pipe and rollers on a lay barge/stinger during laying operations. The contact is modelled by bilinear or nonlinear spring and a bilinear dash pot damper. Possible frictional forces are not included. The contact force is assumed to act normal to the pipe and the roller. It is treated as a discrete element load acting on the pipe, while the contact force acting on the roller is transferred as a nodal load to the stinger. During nonlinear analysis, the computation of the contact force is straightforward and is based on relative positions of the pipe and the roller and the roller stiffness characteristics. In addition, the incremental stiffness relation of the contact force contributes to the system stiffness. The incremental stiffness relation consists of a term associated with the roller spring stiffness and terms corresponding to the change in the contact force direction and point of attack. In linearized dynamic analysis, the contact force is represented by the linearized stiffness at the static equilibrium position. In the following, a consistent formulation of the incremental stiffness relation for the contact force between a pipe element and a roller is developed. Figure 2 defines the geometry of the problem. Figure 2. Problem geometry It is assumed that the contact force between the pipe element and the roller will be normal to both the element and to the roller axis. Possible friction forces are neglected. A unit vector parallel to the contact force is given by \[\boldsymbol{\tilde{n}}=\frac{\Delta \boldsymbol{X}_r\times \Delta \boldsymbol{X}_e}{|\Delta \boldsymbol{X}_r\times \Delta \boldsymbol{X}_e|}\] where \(\mathrm {\begin{array}{l}\Delta \boldsymbol{X}_r=\boldsymbol{X}_{r2}-\boldsymbol{X}_{r1}\\\\\Delta \boldsymbol{X}_e=\boldsymbol{X}_{e2}-\boldsymbol{X}_{e1}\end{array}}\) It is necessary to define the appropriate side of the element relative to the roller. The correct side is defined by a vector, called \(\mathrm {\boldsymbol{d}}\), such that \[\begin{array}{l}\mathrm {if}\quad \boldsymbol{\tilde{n}}^T\boldsymbol{d}>0\quad \mathrm {then}\quad m=1\\\\\mathrm {if}\quad \boldsymbol{\tilde{n}}^T\boldsymbol{d} < 0 \quad \mathrm {then}\quad m=-1\end{array}\] The vector \(\mathrm {\boldsymbol{d}}\) is directed ``from the roller to the element''. That is, the unit direction of a compression contact force acting on the pipe element is \[\boldsymbol{n}=m\boldsymbol{\tilde{n}}\] A contact situation is shown in Figure 3. The roller length, \(\mathrm {l_r}\), and the element length, \(\mathrm {l_e}\), are found by \[\begin{array}{l}l_r=|\Delta \boldsymbol{X}_r|\\\\l_e=|\Delta \boldsymbol{X}_e|\end{array}\] Figure 3. The contact situation. The two contact points, \(\mathrm {\boldsymbol{X}_r}\) and \(\mathrm {\boldsymbol{X}_e}\), are found by setting \[\boldsymbol{X}_e-\boldsymbol{X}_r=\boldsymbol{X}_{e1}+\zeta_e\Delta \boldsymbol{X}_e-\boldsymbol{X}_{r1}-\zeta_r\Delta \boldsymbol{X}_r=v\boldsymbol{n}\] Here, \(\mathrm {v}\) is the distance (may be positive or negative) from \(\mathrm {\boldsymbol{X}_r}\) to \(\mathrm {\boldsymbol{X}_e}\). The above equation may be written as: \[\begin{bmatrix}\boldsymbol{n}_r,&-\boldsymbol{n}_e,&\boldsymbol{n}\end{bmatrix}\begin{bmatrix}l_r\zeta_r\\\\l_e\zeta_e\\\\v\end{bmatrix}=\boldsymbol{X}_{e1}-\boldsymbol{X}_{r1}\] The contact points are now found from inverting Equation (13), which gives \[\begin{bmatrix}l_r\zeta_r\\\\l_e\zeta_e\\\\v\end{bmatrix}=\begin{bmatrix}\displaystyle \frac{m}{\sin(\varphi)}(\boldsymbol{n}_e\times \boldsymbol{n})^T\\\\\displaystyle \frac{m}{\sin(\varphi)}(\boldsymbol{n}_r\times \boldsymbol{n})^T\\\\\boldsymbol{n}^T\end{bmatrix}(\boldsymbol{X}_{e1}-\boldsymbol{X}_{r1})\] A clear restriction is that both \(\mathrm {\zeta_r}\) and \(\mathrm {\zeta_e}\) are within the range \(\mathrm {(0,1)}\). The angle \(\mathrm {\varphi}\) is the angle between the roller and the element and is found by \[\begin{array}{l}\cos(\varphi)=\boldsymbol{n}^T_r\boldsymbol{n}_e\\\\\sin(\varphi)=\sqrt{1-\cos^2(\varphi)}\end{array}\] The magnitude of the contact force (positive as tension) is given as a function of the distance \(\mathrm {v}\), see Figure 4. To include damping, the contact force is also a function of the velocity of the distance \(\mathrm {v}\), \(\mathrm {S=S(v,v_{,t})}\) where \(v_{,t}=\displaystyle \frac{\mathrm {d}v}{\mathrm {d}t}\) (time derivative). Figure 4. Contact force vs contact distance. When the distance is greater than the sum of the roller and pipe external radii, the contact force should be zero (no contact). The contact force vector, see Figure 3, is now found as \[\boldsymbol{S}=S\boldsymbol{n}\] where \(\mathrm {\boldsymbol{S}}\) may be considered as the force vector acting on a nonlinear contact spring, defined in Figure 4. The equilibrium forces acting at the roller and element end points are \[\begin{array}{c}\begin{bmatrix}\boldsymbol{S}_{r1}\\\\\boldsymbol{S}_{r2}\\\\\boldsymbol{S}_{e1}\\\\\boldsymbol{S}_{e2}\end{bmatrix}=\begin{bmatrix}-(1-\zeta_r)\boldsymbol{I}\\\\-\zeta_r\boldsymbol{I}\\\\-(1-\zeta_e)\boldsymbol{I}\\\\-\zeta_e\boldsymbol{I}\end{bmatrix}\boldsymbol{S}\\\\\mathrm {or}\\\\\boldsymbol{S}_{\mathrm {tot}}=\boldsymbol{PS}\end{array}\] The consistent incremental stiffness relation is derived from the above force expression by differentiation of \(\mathrm {\zeta_r}\), \(\mathrm {\zeta_e}\), \(\mathrm {\boldsymbol{n}}\) and \(\mathrm {S}\). Note that a dot ( . ) means differentiation with respect to load level (or time). \[\boldsymbol{\dot {S}}_{\mathrm {tot}}=\boldsymbol{\dot {P}S}+\boldsymbol{P}(\boldsymbol{\dot {n}}S+\boldsymbol{n}\dot {S})\] First, \(\mathrm {\boldsymbol{\dot {P}S}}\) is considered. \[\boldsymbol{\dot {P}S}=\begin{bmatrix}\dot {\zeta}_r\boldsymbol{S}\\\\-\dot {\zeta}_r\boldsymbol{S}\\\\-\dot {\zeta}_e\boldsymbol{S}\\\\\dot {\zeta}_e\boldsymbol{S}\end{bmatrix}\] Differentiation of Equation (13) gives \[\begin{bmatrix}\Delta \boldsymbol{\dot {X}}_r,&-\Delta \boldsymbol{\dot {X}}_e,&\boldsymbol{\dot {n}}\end{bmatrix}\begin{bmatrix}\zeta_r\\\\\zeta_e\\\\v\end{bmatrix}+\begin{bmatrix}\boldsymbol{n}_r,&-\boldsymbol{n}_e,&\boldsymbol{n}\end{bmatrix}\begin{bmatrix}l_r\dot {\zeta}_r\\\\l_e\dot {\zeta}_e\\\\v\end{bmatrix}=\boldsymbol{\dot {X}}_{e1}-\boldsymbol{\dot {X}}_{r1}\] By making use of Equation (14) the following expression is obtained, \[\begin{bmatrix}l_r\dot {\zeta}_r\\\\l_e\dot {\zeta}_e\\\\v\end{bmatrix}=\begin{bmatrix}\displaystyle \frac{m}{\sin(\varphi)}(\boldsymbol{n}_e\times \boldsymbol{n})^T\\\\\displaystyle \frac{m}{\sin(\varphi)}(\boldsymbol{n}_r\times \boldsymbol{n})^T\\\\\boldsymbol{n}^T\end{bmatrix}(\boldsymbol{P}^T\dot {X}_{\mathrm {tot}}-v\boldsymbol{\dot {n}})\] where \[\boldsymbol{\dot {X}}_{\mathrm {tot}}=\begin{bmatrix}\boldsymbol{\dot {X}}_{r1}\\\\\boldsymbol{\dot {X}}_{r2}\\\\\boldsymbol{\dot {X}}_{e1}\\\\\boldsymbol{\dot {X}}_{e2}\end{bmatrix}\] The expression for \(\mathrm {\boldsymbol{\dot {n}}}\) is found by differentiation of Equation (8), and after some manipulation the result is \[\begin{array}{c}\boldsymbol{\dot {n}}=\displaystyle \frac{m}{\sin(\varphi)}\begin{bmatrix}\displaystyle \frac{1}{l_r}\boldsymbol{N}^T_e,&-\displaystyle \frac{1}{l_r}\boldsymbol{N}^T_e,&-\displaystyle \frac{1}{l_e}\boldsymbol{N}^T_r,&\displaystyle \frac{1}{l_e}\boldsymbol{N}^T_r\end{bmatrix}\boldsymbol{\dot {X}}_{\mathrm {tot}}\\\\\mathrm {or}\quad \boldsymbol{\dot {n}}=\boldsymbol{R}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}\end{array}\] where \[\begin{array}{l}\boldsymbol{N}_e=\boldsymbol{n}(\boldsymbol{n}_e\times \boldsymbol{n})^T\\\\\boldsymbol{N}_r=\boldsymbol{n}(\boldsymbol{n}_r\times \boldsymbol{n})^T\end{array}\] Equation (21) now gives: \[\begin{array}{c}\displaystyle \dot {\zeta_r}=\frac{m}{l_r\sin(\varphi)}(\boldsymbol{n}_e\times \boldsymbol{n})^T\boldsymbol{P}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}-\frac{v}{\sin^2(\varphi)}\begin{bmatrix}\displaystyle \frac{1}{l^2_r}\boldsymbol{n}^T,&\displaystyle -\frac{1}{l^2_r}\boldsymbol{n}^T,&\displaystyle -\frac{1}{l_rl_e}\cos(\varphi)\boldsymbol{n}^T,&\displaystyle \frac{1}{l_rl_e}\cos(\varphi)\boldsymbol{n}^T\end{bmatrix}\boldsymbol{\dot {X}}_{\mathrm {tot}}\\\\\displaystyle \dot {\zeta_e}=\frac{m}{l_e\sin(\varphi)}(\boldsymbol{n}_r\times \boldsymbol{n})^T\boldsymbol{P}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}-\frac{v}{\sin^2(\varphi)}\begin{bmatrix}\displaystyle \frac{1}{l_rl_e}\cos(\varphi)\boldsymbol{n}^T,&\displaystyle -\frac{1}{l_rl_e}\cos(\varphi)\boldsymbol{n}^T,&\displaystyle -\frac{1}{l^2_e}\boldsymbol{n}^T,&\displaystyle \frac{1}{l^2_e}\boldsymbol{n}^T\end{bmatrix}\boldsymbol{\dot {X}}_{\mathrm {tot}}\\\\\dot {v}=\boldsymbol{n}^T\boldsymbol{P}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}\\\\\boldsymbol{n}^T\boldsymbol{\dot {n}}=\boldsymbol{n}^T\boldsymbol{R}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}=\boldsymbol{0}\\\\\end{array}\] And, finally, the expression for \(\mathrm {\boldsymbol{\dot {P}S}}\) is \[\begin{array}{c}\boldsymbol{\dot {P}S}=\{\displaystyle \frac{mS}{\sin(\varphi)}\begin{bmatrix}\displaystyle \frac{1}{l_r}\boldsymbol{N}_e\\\\-\displaystyle \frac{1}{l_r}\boldsymbol{N}_e\\\\-\displaystyle \frac{1}{l_e}\boldsymbol{N}_r\\\\\displaystyle \frac{1}{l_e}\boldsymbol{N}_r\end{bmatrix}\boldsymbol{P}^T-\frac{Sv}{\sin^2(\varphi)}\begin{bmatrix}\displaystyle \frac{1}{l^2_r}\boldsymbol{N},&\displaystyle -\frac{1}{l^2_r}\boldsymbol{N},&\displaystyle -\frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N},&\displaystyle \frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N}\\\\\displaystyle -\frac{1}{l^2_r}\boldsymbol{N},&\displaystyle \frac{1}{l^2_r}\boldsymbol{N},&\displaystyle \frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N},&\displaystyle -\frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N}\\\\\displaystyle -\frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N},&\displaystyle \frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N},&\displaystyle \frac{1}{l^2_e}\boldsymbol{N},&\displaystyle -\frac{1}{l^2_e}\boldsymbol{N}\\\\\displaystyle \frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N},&\displaystyle -\frac{\cos(\varphi)}{l_rl_e}\boldsymbol{N},&\displaystyle -\frac{1}{l^2_e}\boldsymbol{N},&\displaystyle \frac{1}{l^2_e}\boldsymbol{N}\end{bmatrix}\}\boldsymbol{\dot {X}}_{\mathrm {tot}}\\\\\mathrm {or}\\\\\boldsymbol{\dot {P}S}=(S\boldsymbol{RP}^T\displaystyle -\frac{Sv}{\sin^2(\varphi)}\boldsymbol{A})\boldsymbol{\dot {X}}_{\mathrm {tot}}\\\\\end{array}\] Note that matrix \(\mathrm {\boldsymbol{A}}\) is symmetric, and \(\mathrm {\boldsymbol{N}}\) is defined in Equation (33). The expression for \(\mathrm {\boldsymbol{P}(\boldsymbol{\dot {n}}S+\boldsymbol{n}\dot {S})}\), Equation (18), will not be found. \[\boldsymbol{P}\boldsymbol{\dot {n}}S=S\boldsymbol{PR}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}\] To evaluate \(\mathrm {\boldsymbol{P\dot {n}}S}\) it must be considered that the contact force is a function of both \(\mathrm {v}\) and \(\mathrm {v_{,t}}\) which implies \[S=S(v)+\bar{c}v_{,t}\] where \(\mathrm {\bar{c}}\) is a constant damping coefficient. Differentiation yields \[\dot {S}=k\dot {v}+\bar{c}\dot {v}_{,t}\] By time differentiation of \(\mathrm {\dot {v}}\), Equation (25), it is found that \[\dot {v}_{,t}=\boldsymbol{n}^T_{,t}\boldsymbol{P}^T\boldsymbol{\dot {X}}_{\mathrm {tot}}+\boldsymbol{n}^T\boldsymbol{P}^T_{,t}\boldsymbol{\dot {X}}_{\mathrm {tot}}+\boldsymbol{n}^T_{,t}\boldsymbol{P}^T\boldsymbol{\dot {X}}_{\mathrm {tot,t}}\] or, by using the relations \(\boldsymbol{n}^T_{,t}=\boldsymbol{X}_{\mathrm {tot},t}\boldsymbol{R}\) from Equation (23) and \(\boldsymbol{n}^T\boldsymbol{P}^T_{,t}=\boldsymbol{X}^T_{\mathrm {tot},t}(\boldsymbol{PR}^T-\displaystyle \frac{v}{\sin^2(\varphi)}\boldsymbol{A})\) from Equation (26), it is found that \[\dot {v}_{,t}=\boldsymbol{X}^T_{\mathrm {tot},t}\begin{bmatrix}\boldsymbol{RP}^T+\boldsymbol{PR}^T\displaystyle -\frac{v}{\sin^2(\varphi)}\boldsymbol{A}\end{bmatrix}\boldsymbol{\dot {X}}_{\mathrm {tot}}+\boldsymbol{n}^T\boldsymbol{P}^T\boldsymbol{X}_{\mathrm {tot},t}\] By combining the above results, the final result is \[\begin{array}{l} \boldsymbol{\dot {S}}_{\mathrm {tot}}= \begin{bmatrix}k\boldsymbol{PNP}^T+S(\boldsymbol{RP}^T+\boldsymbol{PR}^T\displaystyle -\frac{v}{\sin^2(\varphi)}\boldsymbol{A}) \end{bmatrix}\boldsymbol{\dot {X}}_{\mathrm {tot}}\\ + \bar{c}\boldsymbol{PnX}^T_{\mathrm {tot},t}(\boldsymbol{RP}^T+\boldsymbol{PR}^T\displaystyle -\frac{v}{\sin^2(\varphi)}\boldsymbol{A})\boldsymbol{\dot {X}}_{\mathrm {tot}}\\ + \bar{c}\boldsymbol{PNP}^T\boldsymbol{\dot {X}}_{\mathrm {tot},t} \end{array}\] where \[\boldsymbol{N}=\boldsymbol{nn}^T\] The above incremental relation may be rewritten as \[\boldsymbol{\dot {S}}_{\mathrm {tot}}=\boldsymbol{K}_{Ik}\boldsymbol{\dot {X}}_{\mathrm {tot}}+\boldsymbol{K}_{Iv}\boldsymbol{\dot {X}}_{\mathrm {tot}}+\boldsymbol{C}_I\boldsymbol{\dot {X}}_{\mathrm {tot},t}\] where \(\mathrm {\boldsymbol{K}_{Ik}\quad }\) : incremental stiffness matrix from material and geometry stiffness \(\mathrm {\boldsymbol{K}_{Iv}\quad }\) : incremental stiffness matrix from damping (unsymmetric) \(\mathrm {\boldsymbol{C}_I\quad \:\:\,}\) : incremental damping matrix 5. Tensioner Contact The main intention of the tensioner contact formulation is to enable the user to model applied tension from a tensioner acting on a pipe during laying operations. The force from the tensioner acting on the pipe is directed along the pipe and should be within an upper and lower limit. Due to possible dynamic behaviour at the lay barge and/or the pipe, the tensioner force may increase/decrease, but it will not exceed the upper limit nor decrease below the lower limit. Between the upper and lower limit, the tensioner is described as a linear spring having a relatively large stiffness. During static analyses, a constant tensioner force directed along the pipe is applied. The force should be within the upper and the lower limit. In nonlinear dynamic analysis, the tensioner is modelled by the actual tensioner force and incremental stiffness based on the tensioner spring stiffness. During linearized analysis, the tensioner force is represented by the linearized stiffness at static equilibrium. In this case the applied tensioner force may pass the limits. In the following, the expressions for tensioner force and incremental stiffness are developed. See Figure 5. The tensioner is defined to be fixed to the stinger at a specified stinger point, \(\mathrm {P}\), and to follow the translations and rotations of that point. The intersection point, \(\mathrm {C}\), between the tensioner and the beam is uniquely defined by the coordinate of \(\mathrm {P}\), the normal direction of the mid-plane through the tensioner at \(\mathrm {P}\) and the straight line between the beam nodal points \(\mathrm {A}\) and \(\mathrm {B}\). The distance from point \(\mathrm {A}\) to the intersection point \(\mathrm {C}\) is found to be expressed as \[\xi \boldsymbol{L}_{AB}=\frac{\boldsymbol{n}^T_p(\boldsymbol{X}_p-\boldsymbol{X}_A)}{\boldsymbol{n}^T_p\boldsymbol{n}}\] We now define the interaction force, \(\mathrm {F_D}\) , between the stinger and the beam element to be directed along \(\mathrm {-\boldsymbol{n}}\). The corresponding equilibrium forces directed along \(\mathrm {+\boldsymbol{n}}\) at nodes \(\mathrm {A}\) and \(\mathrm {B}\) are respectively: \[\begin{array}{l}\begin{bmatrix}F_A\\\\F_B\end{bmatrix}=\begin{bmatrix}1-\xi \\\\\xi \end{bmatrix}F_D\\\\F_A+F_B=F_D\end{array}\] Figure 5. The tensioner model. Similarly, the incremental displacement along \(\mathrm {\boldsymbol{n}}\) at point \(\mathrm {D}\) is expressed by linear combinations of the incremental nodal displacements as \[\dot {U}_D=\begin{bmatrix}1-\xi ,&\xi \end{bmatrix}\begin{bmatrix}\dot {U}_A\\\\\dot {U}_B\end{bmatrix}\] The incremental nodal displacements along \(\mathrm {\boldsymbol{n}}\) may be expressed by the incremental global displacements as follows. \[\begin{bmatrix}\dot {U}_A\\\\\dot {U}_B\end{bmatrix}=\begin{bmatrix}\boldsymbol{n}^T&\boldsymbol{0}\\\\\boldsymbol{0}&\boldsymbol{n}^T\end{bmatrix}\begin{bmatrix}\boldsymbol{\dot {U}}_A\\\\\boldsymbol{\dot {U}}_B\end{bmatrix}\] In the global system the nodal forces at \(\mathrm {A}\) and \(\mathrm {B}\) become \[\begin{bmatrix}\boldsymbol{F}_A\\\\\boldsymbol{F}_B\end{bmatrix}=\begin{bmatrix}\boldsymbol{n}&\boldsymbol{0}\\\\\boldsymbol{0}&\boldsymbol{n}\end{bmatrix}\begin{bmatrix}F_A\\\\F_B\end{bmatrix}\] or \[\begin{bmatrix}\boldsymbol{F}_A\\\\\boldsymbol{F}_B\end{bmatrix}=\begin{bmatrix}\boldsymbol{n}&\boldsymbol{0}\\\\\boldsymbol{0}&\boldsymbol{n}\end{bmatrix}\begin{bmatrix}1-\xi \\\\\xi \end{bmatrix}F_D\] An approximation of an incremental version of the above equilibrium equation may be written \[\begin{bmatrix}\boldsymbol{\dot {F}}_A\\\\\boldsymbol{\dot {F}}_B\end{bmatrix}=\begin{bmatrix}\boldsymbol{n}&\boldsymbol{0}\\\\\boldsymbol{0}&\boldsymbol{n}\end{bmatrix}\begin{bmatrix}1-\xi \\\\\xi \end{bmatrix}\dot {F}_D\] where only the change in \(\mathrm {F_D}\) is considered. Further, \[\dot {F}_D=k\dot {U}_{CD}\cong k\dot {U}_D\] where \(\mathrm {k}\) is the slope of the force-displacement curve in Figure 6 (bilinear). Note that the change in the location of point \(\mathrm {C}\) is not considered in the incremental stiffness. Figure 6. Tensioner force model. Again, neglecting the change in \(\mathrm {\xi }\), an approximation to the incremental stiffness matrix of the tensioner is obtained as \[\begin{bmatrix} \boldsymbol{\dot {F}}_A\\\\\boldsymbol{\dot {F}}_B \end{bmatrix}=k \begin{bmatrix} \boldsymbol{n}&\boldsymbol{0}\\\\\boldsymbol{0}&\boldsymbol{n} \end{bmatrix} \begin{bmatrix} 1-\xi \\\\\xi \end{bmatrix} \begin{bmatrix}1-\xi ,&\xi \end{bmatrix}\begin{bmatrix}\boldsymbol{n}^T&\boldsymbol{0}\\\\\boldsymbol{0}&\boldsymbol{n}^T\end{bmatrix}\begin{bmatrix}\boldsymbol{\dot {U}}_A\\\\\boldsymbol{\dot {U}}_B \end{bmatrix}\] Or, after multiplication \[\begin{bmatrix} \boldsymbol{\dot {F}}_A\\\\\boldsymbol{\dot {F}}_B \end{bmatrix}=k \begin{bmatrix} \boldsymbol{nn}^T(1-\xi )^2,&\boldsymbol{nn}^T(1-\xi )\xi \\\\\boldsymbol{nn}^T(1-\xi )\xi ,&\boldsymbol{nn}^T\xi ^2 \end{bmatrix} \begin{bmatrix} \boldsymbol{\dot {U}}_A\\\\\boldsymbol{\dot {U}}_B \end{bmatrix}\] Note that the value of \(\mathrm {k}\) and the force, \(\mathrm {F_D}\), from the tensioner on the beam are found from the current configuration (updated in the equilibrium iterations) relative to the last equilibrium configuration. 6. References Holte, K.H. and Ormberg, H. (1995): ICPS Component Modelling, PINTIES Report ZA-E 70 012, Trondheim. Ormberg, H. and Holthe, K.H. (1995): ICPS Development, PINTIES Report ZA-E 70 013, Trondheim. Voie, P. E., Skeie, G. and Bergan-Haavik, J. (2014): Importance Rating of Riser-Soil Interaction Effects OMAE2014-24179, ASME 2014 33rd International Conference on Ocean, Offshore and Arctic Engineering. San Francisco, 2014 Byrne, B. W., et al. (2017): PISA: new design methods for offshore wind turbine monopiles Offshore Site Investigation Geotechnics 8th International Conference Proceeding. Vol. 142. No. 161. Society for Underwater Technology, 2017. Byrne, B. W., et al. (2019): PISA design methods for offshore wind turbine monopiles Offshore Technology Conference. Offshore Technology Conference, 2019. Burd, H. J., et al. (2020): PISA design model for monopiles for offshore wind turbines: application to a marine sand, Géotechnique 70.11: 1048-1066 [https://doi.org/10.1680/jgeot.18.P.277]