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 inplane contact modelled by a combination of springs and frictional forces.
The soil contact formulation enables modelling monopile foundations of marine structures by springs and inplane 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 reestablished 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. Inplane springfriction 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 bilinear springs or with the nonlinear consolidated risersoil model.
2.1.1. Bilinear spring contact
The riserseabed 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
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
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
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
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
Damping may be handled similarly.
2.1.2. Consolidated risersoil model
The Consolidated risersoil 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 inplane lateral direction.
A risersoil 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 risersoil contact normal to the seafloor is not influenced by the inplane springfriction contact. This means that the soil capacity is not reduced due to nonnormal contact. The inplane friction forces will, however, depend on the normal contact force.
In the present version, the risersoil contact element memory is activated together with the activation of the inplane seafloor contact; i.e. the static load "FRIC".
2.2. Inplane contact
Inplane 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 inplane seafloor contact is in general modelled as shown in Figure 1.
For the inplane 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
and the interaction force, \(\mathrm {F_i}\), as
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 subsurface 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 reused 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 (crosssection rotation).
The added stiffness is modelled through nonlinearelastic 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 diametertoembeddedlength 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.
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
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
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
A contact situation is shown in Figure 3. The roller length, \(\mathrm {l_r}\), and the element length, \(\mathrm {l_e}\), are found by
The two contact points, \(\mathrm {\boldsymbol{X}_r}\) and \(\mathrm {\boldsymbol{X}_e}\), are found by setting
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:
The contact points are now found from inverting Equation (13), which gives
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
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).
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
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
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).
First, \(\mathrm {\boldsymbol{\dot {P}S}}\) is considered.
Differentiation of Equation (13) gives
By making use of Equation (14) the following expression is obtained,
where
The expression for \(\mathrm {\boldsymbol{\dot {n}}}\) is found by differentiation of Equation (8), and after some manipulation the result is
where
Equation (21) now gives:
And, finally, the expression for \(\mathrm {\boldsymbol{\dot {P}S}}\) is
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.
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
where \(\mathrm {\bar{c}}\) is a constant damping coefficient.
Differentiation yields
By time differentiation of \(\mathrm {\dot {v}}\), Equation (25), it is found that
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
By combining the above results, the final result is
where
The above incremental relation may be rewritten as
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 midplane 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
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:
Similarly, the incremental displacement along \(\mathrm {\boldsymbol{n}}\) at point \(\mathrm {D}\) is expressed by linear combinations of the incremental nodal displacements as
The incremental nodal displacements along \(\mathrm {\boldsymbol{n}}\) may be expressed by the incremental global displacements as follows.
In the global system the nodal forces at \(\mathrm {A}\) and \(\mathrm {B}\) become
or
An approximation of an incremental version of the above equilibrium equation may be written
where only the change in \(\mathrm {F_D}\) is considered. Further,
where \(\mathrm {k}\) is the slope of the forcedisplacement curve in Figure 6 (bilinear). Note that the change in the location of point \(\mathrm {C}\) is not considered in the incremental stiffness.
Again, neglecting the change in \(\mathrm {\xi }\), an approximation to the incremental stiffness matrix of the tensioner is obtained as
Or, after multiplication
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 ZAE 70 012, Trondheim.
Ormberg, H. and Holthe, K.H. (1995): ICPS Development, PINTIES Report ZAE 70 013, Trondheim.
Voie, P. E., Skeie, G. and BerganHaavik, J. (2014): Importance Rating of RiserSoil Interaction Effects OMAE201424179, 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: 10481066 [https://doi.org/10.1680/jgeot.18.P.277]