Dynamic Time Domain Analysis

1. General

The dynamic equilibrium of a spatially discretized finite element system model can in general be expressed as

\[\boldsymbol{R}^I(\boldsymbol{r},\boldsymbol{\ddot {r}},\boldsymbol{t})+\boldsymbol{R}^D(\boldsymbol{r},\boldsymbol{\dot {r}},\boldsymbol{t})+\boldsymbol{R}^S(\boldsymbol{r},\boldsymbol{t})=\boldsymbol{R}^E(\boldsymbol{r},\boldsymbol{\dot {r}},\boldsymbol{t})\]

where

  • \(\mathrm {\boldsymbol{R}^I\qquad\,}\) - inertia force vector

  • \(\mathrm {\boldsymbol{R}^D\qquad}\) - damping force vector

  • \(\mathrm {\boldsymbol{R}^S\qquad\,}\) - internal structural reaction force vector

  • \(\mathrm {\boldsymbol{R}^E\qquad}\) - external force vector

  • \(\mathrm {\boldsymbol{r}}\), \(\mathrm {\boldsymbol{\dot {r}}}\), \(\mathrm {\boldsymbol{\ddot {r}}\quad }\)- structural displacement, velocity and acceleration vectors

This is a nonlinear system of differential equations due to the displacement dependencies in the inertia and damping forces, and the coupling between the external load vector and structural displacement and velocity. In addition, there is a nonlinear relationship between internal forces and displacements as discussed for the static situation, see Static Finite Element Analysis.

All force vectors are established by assembly of element contributions and specified discrete nodal forces. A description of the internal reaction force vector is given in Static Finite Element Analysis, the other terms in Equation (1) are further detailed in the following.

The external force vector accounts for:

1.1. Inertia force

The inertia force vector can be expressed as

\[\boldsymbol{R}^I(\boldsymbol{r},\boldsymbol{\ddot {r}},\boldsymbol{t})=\Big[\boldsymbol{M}^S+\boldsymbol{M}^F(\boldsymbol{r})+\boldsymbol{M}^H(\boldsymbol{r})\Big]\boldsymbol{\ddot {r}}\]

where

  • \(\mathrm {\boldsymbol{M}^S\qquad\:\:}\) - structural mass matrix

  • \(\mathrm {\boldsymbol{M}^F(\boldsymbol{r})\quad }\) - mass matrix accounting for internal fluid flow. Will be displacement- and time-dependent if slug flow is considered

  • \(\mathrm {\boldsymbol{M}^H(\boldsymbol{r})\quad }\) - displacement-dependent hydrodynamic mass matrix accounting for the structural acceleration terms in the Morison equation as added mass contributions in local directions, see Hydrodynamic Load Models for Submerged Elements.

For the general beam cross section CRS7, an additional intertia force vector representing the Coriolis centripetal load is included. It is calculated based on the rotational diogonal terms of the local element mass matrix and the rotational speeds around local axes as described in Equation (3).

\[ boldsymbol{S}^C= \begin{bmatrix}S_{x1}^c\\S_{y1}^c\\S_{z1}^c\\S_{\theta x1}^c\\S_{\theta y1}^c\\S_{\theta z1}^c\\S_{x2}^c\\S_{y2}^c\\S_{z2}^c\\S_{\theta x2}^c\\S_{\theta y2}^c\\S_{\theta z2}^c\end{bmatrix}= \begin{bmatrix} 0\\ 0\\ 0\\ (m_{6,6}-m_{5,5})\dot \theta _{y1}\dot \theta _{z1}\\ (m_{4,4}-m_{6,6})\dot \theta _{x1}\dot \theta _{z1}\\ (m_{5,5}-m_{4,4})\dot \theta _{x1}\dot \theta _{y1}\\ 0\\ 0\\ 0\\ (m_{12,12}-m_{11,11})\dot \theta _{y2}\dot \theta _{z2}\\ (m_{10,10}-m_{12,12})\dot \theta _{x2}\dot \theta _{z2}\\ (m_{11,11}-m_{10,10})\dot \theta _{x2}\dot \theta _{y2}\\ \end{bmatrix}\]

The local element Coriolis load vectors are included in the system inertia force vector, Equation (2).

1.2. Damping force

The damping force vector is expressed as

\[\boldsymbol{R}^D(\boldsymbol{r},\boldsymbol{\dot {r}})=\Big[\boldsymbol{C}^S(\boldsymbol{r})+\boldsymbol{C}^H(\boldsymbol{r})+\boldsymbol{C}^D(\boldsymbol{r},\boldsymbol{\dot {r}})\Big]\boldsymbol{\dot {r}}\]

where

  • \(\mathrm {\boldsymbol{C}^S(\boldsymbol{r})\quad \quad \:}\) - internal structural damping matrix

  • \(\mathrm {\boldsymbol{C}^H(\boldsymbol{r})\quad \quad }\) - hydrodynamic damping matrix accounting for diffraction effects for floating, partly submerged elements (irrelevant in other situations)

  • \(\mathrm {\boldsymbol{C}^D(\boldsymbol{r},\boldsymbol{\dot {r}})\quad }\) - matrix of specified discrete dashpot dampers which may be displacement- (and velocity-) dependent

Internal structural damping must be included in the dynamic analysis to account for energy dissipation in the structure itself. The physical nature of the structural damping is strongly dependent on the cross sectional properties of the structure, e.g. bonded and unbonded flexible pipes, umbilicals, anchor lines composed of wire and/or chain segments and steel risers.

Measurements of structural damping of flexible pipes have shown that different structural damping levels can be related to axial, torsional and bending deformations. A rather low linear viscous type of damping is typically related to axial and torsional deformations while an amplitude-dependent hysteretic type of damping is related to large bending deformations. Several structural damping models are therefore optionally available to enable description of a realistic structural damping level in dynamic analysis. Further details are given in Section 6.

The most important nonlinear effects that must be considered in dynamic analyses of slender marine structures are:

  • Geometric stiffness, i.e. contribution from axial force to transverse stiffness.

  • Nonlinear material properties.

  • Hydrodynamic loading according to the generalized Morison equation expressed by relative velocities.

  • Integration of loading to actual surface elevation.

  • Contact problems, e.g. bottom contact, riser collision, vessel/riser, pipe/stinger and contact.

A non-Gaussian response must therefore in general be foreseen for systems exposed to irregular loading due to wave and vessel motions. The non-Gaussian response characteristics are, however, strongly system - and excitation dependent (Sødahl 1991, 1992).

The distinguishing feature among available analysis methods is how nonlinearities are treated. The following finite element analysis techniques are available:

1.2.1. Nonlinear time domain analysis

Step-by-step numerical integration of the incremental dynamic equilibrium equations, with a Newton-Raphson type of equilibrium iteration at each time step. This approach allows for a proper treatment of all the described nonlinearities. Nonlinear dynamic analysis is, however, rather time consuming due to repeated assembly of system matrices (mass, damping and stiffness) and triangularisation during the iteration process at each time step.

1.2.2. Linearized time domain analysis

Step-by-step numerical integration of the dynamic equations found by linearization of mass, damping and stiffness matrices at the static equilibrium position. This means that the system matrices are kept constant throughout the analysis. Nonlinear hydrodynamic loading is, however, included. The linearized time domain analysis gives a significant reduction in computation time when compared to the nonlinear approach and is an attractive alternative for systems where the hydrodynamic loading is the dominant nonlinear effect.

The nonlinear and linearized time domain analyses are further detailed in Section 3 and Section 4, respectively.

2. Methods for Numerical Time Integration

The step-by-step numerical integration of the dynamic equilibrium equations is based on the well-known Newmark \(\mathrm {\beta }\)-family including the Wilson \(\mathrm {\theta }\)-method considering a constant time step throughout the analysis. Both methods can be applied for nonlinear as well as linearized analysis.

These methods apply the following relations between displacements, velocity and acceleration vectors at time \(\mathrm {t}\) and \(\mathrm {t+\Delta \tau}\) :

\[\begin{array}{l}\displaystyle \boldsymbol{\dot {r}}_{t+\Delta \tau}=\boldsymbol{\dot {r}}_t+(1-\gamma )\boldsymbol{\ddot {r}}_t\Delta \tau+\gamma \boldsymbol{\ddot {r}}_{t+\Delta \tau}\Delta \tau\\\\\displaystyle \boldsymbol{r}_{t+\Delta \tau}=\boldsymbol{r}_t+\boldsymbol{\dot {r}}_t\Delta \tau+(\frac{1}{2}-\beta )\boldsymbol{\ddot {r}}_t(\Delta \tau)^2+\beta \boldsymbol{\ddot {r}}_{t+\Delta \tau}(\Delta \tau)^2\end{array}\]

where \(\mathrm {\Delta \tau=\theta \Delta t}\), \(\mathrm {\theta \geq1.0}\).

\(\mathrm {\gamma }\), \(\mathrm {\beta }\) and \(\mathrm {\theta }\) are parameters in the integration methods defining the functional change in displacement, velocity and acceleration vectors over the time step \(\mathrm {\Delta t}\).

The Newmark \(\mathrm {\beta }\)-family is found for \(\mathrm {\theta =1}\). Values of \(\mathrm {\gamma }\), \(\mathrm {\beta }\) and \(\mathrm {\theta }\) corresponding to well-known time integration techniques are summarized in Table 1.

Table 1. Integration methods covered by the Newmark \(\mathrm {\beta }\)-family including Wilson \(\mathrm {\theta }\)-method
Method \(\mathrm {\gamma }\) \(\mathrm {\beta }\) \(\mathrm {\theta }\)

Central difference method

1/2

0

1

Fox-Goodwin method

1/2

1/12

1

Linear acceleration method

1/2

1/6

1

Constant average acceleration method

1/2

1/4

1

Wilson \(\mathrm {\theta }\)-method

1/2

1/6

> 1

The value of \(\mathrm {\gamma }\) determines the artificial (numerical) damping of the method:

  • \(\mathrm {\gamma >1/2}\) : \(\mathrm {\quad }\) positive numerical damping

  • \(\mathrm {\gamma <1/2}\) : \(\mathrm {\quad }\) negative numerical damping

  • \(\mathrm {\gamma =1/2}\) : \(\mathrm {\quad }\) no numerical damping .

\(\mathrm {\gamma =1/2}\) is normally used to obtain second order accuracy.

The integration methods are implicit for \(\mathrm {\beta >0}\), which is assumed in the following. Equation (5) can be rearranged to express the relations between incremental displacement, velocity and acceleration vectors over the time interval \(\mathrm {\Delta \tau}\):

\[\begin{array}{l}\displaystyle \Delta \boldsymbol{\dot {r}}_t=\boldsymbol{\dot {r}}_{t+\Delta \tau}-\boldsymbol{\dot {r}}_t=\frac{\gamma }{\beta \Delta \tau}\Delta \boldsymbol{r}_t-\frac{\gamma }{\beta }\boldsymbol{\dot {r}}_t-(\frac{\gamma }{2\beta }-1)\Delta \tau\boldsymbol{\ddot {r}}_t\\\\\displaystyle \Delta \boldsymbol{\ddot {r}}_t=\boldsymbol{\ddot {r}}_{t+\Delta \tau}-\boldsymbol{\ddot {r}}_t=\frac{1}{\beta (\Delta \tau)^2}\Delta \boldsymbol{r}_t-\frac{1}{\beta \Delta \tau}\boldsymbol{\dot {r}}_t-\frac{1}{2\beta }\boldsymbol{\ddot {r}}_t\end{array}\]

The Wilson \(\mathrm {\theta }\)-method interpolates the acceleration at time \(\mathrm {t+\Delta t}\) between the values at times \(\mathrm {t}\) and \(\mathrm {t+\Delta \tau}\) according to

\[\boldsymbol{\ddot {r}}_{t+\Delta t}=\boldsymbol{\ddot {r}}_t+\frac{1}{\theta }(\boldsymbol{\ddot {r}}_{t+\Delta \tau}-\boldsymbol{\ddot {r}}_t)\]

The corresponding displacement and velocity vectors at time \(\mathrm {t+\Delta t}\), \(\mathrm {\boldsymbol{r}_{t+\Delta t}}\) and \(\mathrm {\boldsymbol{\dot {r}}_{t+\Delta t}}\), are found according to Equation (5) where \(\mathrm {\Delta \tau}\) is substituted by \(\mathrm {\Delta t}\).

The velocity and acceleration vectors corresponding to prescribed displacements are calculated applying Equation (6).

Further details regarding the numerical time integration methods are given in standard textbooks on structural dynamics, see for instance Langen and Sigbjørnsson (1979) and Clough and Penzien (1975).

3. Nonlinear Dynamic Analysis

This section will give an introduction to the method of analysis used in nonlinear dynamic analysis. The presentation mainly follows the approach outlined by Remseth (1978) and Langen and Sigbjørnsson (1978).

The incremental form of the dynamic equilibrium equation, Equation (1), is obtained by considering dynamic equilibrium at two configurations a short time interval \(\mathrm {\Delta t}\) apart:

\[(\boldsymbol{R}^I_{t+\Delta t}-\boldsymbol{R}^I_t)+(\boldsymbol{R}^D_{t+\Delta t}-\boldsymbol{R}^D_t)+(\boldsymbol{R}^S_{t+\Delta t}-\boldsymbol{R}^S_t)=(\boldsymbol{R}^E_{t+\Delta t}-\boldsymbol{R}^E_t)\]

Equation Equation (8) states that the increment in external loading is balanced by increments in inertia-, damping- and structural reaction forces over the time interval \(\mathrm {\Delta t}\).

For numerical solution, the nonlinear incremental equation of motion is linearized by introducing the tangential mass-, damping- and stiffness matrices at the start of the increment. The linearized incremental equation of motion can be expressed as

\[\boldsymbol{M}_t\Delta \boldsymbol{\ddot {r}}_t+\boldsymbol{C}_t\Delta \boldsymbol{\dot {r}}_t+\boldsymbol{K}_t\Delta \boldsymbol{r}_t=\Delta \boldsymbol{R}^E_t\]

where \(\mathrm {\boldsymbol{M}_t}\), \(\mathrm {\boldsymbol{C}_t}\) and \(\mathrm {\boldsymbol{K}_t}\) denote the tangential mass-, damping- and stiffness matrices computed at time \(\mathrm {t}\). \(\mathrm {\Delta \boldsymbol{r}_t}\), \(\mathrm {\Delta \boldsymbol{\dot {r}}_t}\), \(\mathrm {\Delta \boldsymbol{\ddot {r}}_t}\) and \(\mathrm {\Delta \boldsymbol{R}^E_t}\) are the incremental displacement, velocity, acceleration and external force vectors, respectively, i.e.

\[\begin{array}{l}\Delta \boldsymbol{r}_t=\boldsymbol{r}_{t+\Delta t}-\boldsymbol{r}_t\\\\\displaystyle \Delta \boldsymbol{\dot {r}}_t=\boldsymbol{\dot {r}}_{t+\Delta t}-\boldsymbol{\dot {r}}_t\\\\\Delta \boldsymbol{\ddot {r}}_t=\boldsymbol{\ddot {r}}_{t+\Delta t}-\boldsymbol{\ddot {r}}_t\\\\\Delta \boldsymbol{R}^E_t=\boldsymbol{R}^E_{t+\Delta t}-\boldsymbol{R}^E_t\end{array}\]

Due to nonlinearities, the dynamic equilibrium equation is not satisfied at the end of the time step. The residual force vector, due to change in mass, damping and stiffness over the time step is given by

\[\delta \boldsymbol{R}_t=\boldsymbol{R}^E_{t+\Delta t}-(\boldsymbol{R}^I_{t+\Delta t}+\boldsymbol{R}^D_{t+\Delta t}+\boldsymbol{R}^S_{t+\Delta t})\]

To prevent error accumulation, the residual force vector is added to the incremental equilibrium equation at the next time step. Thus, the incremental equation of motion including equilibrium correction is written as

\[\boldsymbol{M}_t\Delta \boldsymbol{\ddot {r}}_t+\boldsymbol{C}_t\Delta \boldsymbol{\dot {r}}_t+\boldsymbol{K}_t\Delta \boldsymbol{r}_t=\boldsymbol{R}^E_{t+\Delta t}-(\boldsymbol{R}^I_t+\boldsymbol{R}^D_t+\boldsymbol{R}^S_t)\]

By rewriting Equation (12) for dynamic equilibrium at time \(\mathrm {t+\Delta \tau}\) and inserting Equation (6), the incremental equation expressed by the incremental displacement vector over the time interval \(\mathrm {\Delta \tau+\theta \Delta t}\), is written

\[\boldsymbol{\hat{K}}_t\Delta \boldsymbol{r}_t=\Delta \boldsymbol{\hat{R}}_t\]

where the effective stiffness, \(\mathrm {\boldsymbol{\hat{K}}_t}\), and effective incremental load vector, \(\mathrm {\Delta \boldsymbol{\hat{R}}_t}\), are defined by the relations

\[\boldsymbol{\hat{K}}_t=\frac{1}{\displaystyle \beta (\Delta \tau)^2}\boldsymbol{M}_t+\frac{\gamma }{\beta \Delta \tau}\boldsymbol{C}_t+\boldsymbol{K}_t\]
\[\Delta \boldsymbol{\hat{R}}_t=\boldsymbol{R}^E_{t+\Delta \tau}-(\boldsymbol{R}^I_t+\boldsymbol{R}^D_t+\boldsymbol{R}^S_t)+(\frac{1}{\beta \Delta \tau}\boldsymbol{\dot {r}}_t+\frac{1}{2\beta }\boldsymbol{\ddot {r}}_t)+\boldsymbol{C}_t(\frac{\gamma }{\beta }\boldsymbol{\dot {r}}_t+(\frac{\gamma }{2\beta }-1)\Delta \tau\boldsymbol{\ddot {r}}_t)\]

To obtain dynamic equilibrium at the end of the time step, an iterative approach similar to the static equilibrium iteration is applied. Equilibrium iterations are performed according the equation

\[\Big[{^{i-1}}\boldsymbol{R}^I_{t+\Delta \tau} +{^i}\Delta \boldsymbol{R}^I_{t+\Delta \tau}\Big] +\Big[{^{i-1}}\boldsymbol{R}^D_{t+\Delta \tau} +{^i}\Delta \boldsymbol{R}^D_{t+\Delta \tau}\Big] +\Big[{^{i-1}}\boldsymbol{R}^S_{t+\Delta \tau} +{^i}\Delta \boldsymbol{R}^S_{t+\Delta \tau}\Big] ={^{i-1}}\boldsymbol{R}^E_{t+\Delta \tau}\]

The sums in brackets represent the total inertia-, damping- and structural reaction force vectors, respectively, after iteration number \(\mathrm {i}\).

Introducing the displacement correction vector \(\mathrm {\boldsymbol{^i\Delta _r}}\) the improved incremented displacement vector after iteration number \(\mathrm {i}\) is given by

\[^i\Delta \boldsymbol{r}_t={^{i-1}}\Delta \boldsymbol{r}_t+\boldsymbol{{^i}\Delta _r}\]

and the total displacement vector reads

\[^i\boldsymbol{r}_{t+\Delta \tau}={^{i-1}}\boldsymbol{r}_{t+\Delta \tau}+\boldsymbol{{^i}\Delta _r}\]

The corresponding velocity and acceleration vectors are given by

\[\begin{array}{l}\displaystyle ^i\boldsymbol{\dot {r}}_{t+\Delta \tau}=\,^{i-1}\boldsymbol{\dot {r}}_{t+\Delta \tau}+\frac{\gamma }{\beta \Delta \tau}\:\,\boldsymbol{^i\Delta _r}\\\\\displaystyle ^i\boldsymbol{\ddot {r}}_{t+\Delta \tau}=\,^{i-1}\boldsymbol{\ddot {r}}_{t+\Delta \tau}+\frac{\gamma }{\beta (\Delta \tau)^2}\:\,\boldsymbol{^i\Delta _{r}}\\\\\end{array}\]

The dynamic equilibrium at \(\mathrm {t+\Delta \tau}\), Equation (16), using the relation in Equation (6) is rewritten in the form

\[\begin{array}{l}\displaystyle (\frac{1}{\beta (\Delta \tau)^2}\:\,^{i-1}\boldsymbol{M}_{t+\Delta \tau}+\frac{\gamma }{\beta \Delta \tau}\:\,^{i-1}\boldsymbol{C}_{t+\Delta \tau}+\:\,^{i-1}\boldsymbol{K}_{t+\Delta \tau})\,\boldsymbol{^i\Delta _r}\\\\\displaystyle =\boldsymbol{R}^I_{t+\Delta \tau}+(\,^{i-1}\boldsymbol{R}^I_{t+\Delta \tau}+\,^{i-1}\boldsymbol{R}^D_{t+\Delta \tau}+\,^{i-1}\boldsymbol{R}^S_{t+\Delta \tau})\end{array}\]

The structural reaction force vector, \(\mathrm {^{i-1}\boldsymbol{R}^S_{t+\Delta \tau}}\), is calculated from the state of stress in the elements after iteration \(\mathrm {i-1}\). The inertia and damping force vectors after iterations \(\mathrm {i-1}\), are calculated according to

\[\begin{array}{l}^{i-1}\boldsymbol{R}^I_{t+\Delta \tau}=\boldsymbol{M}_{t+\Delta \tau}\,^{i-1}\boldsymbol{\ddot {r}}_{t+\Delta \tau}\\\\^{i-1}\boldsymbol{R}^D_{t+\Delta \tau}=\boldsymbol{C}_{t+\Delta \tau}\,^{i-1}\boldsymbol{\dot {r}}_{t+\Delta \tau}\end{array}\]

In Equation (20) it is assumed that the tangential mass, damping and stiffness matrices are recalculated at each iteration cycle, which will give a Newton-Raphson iteration procedure. A modified Newton-Raphson procedure can also be included by keeping the matrices constant over several iteration cycles.

Contributions to the external load vector from prescribed displacements due to vessel motions are applied at each time step. There are, however, no contributions from prescribed displacements in the external load vector during the equilibrium iteration.

The criterion used for termination of the iteration is a modified Euclidian displacement norm similar to the norm used in static finite element analysis, see Incremental Equilibrium Iterations.

The iteration procedure will give a stable numerical solution provided that the time step is small enough to give a good starting point for the iteration and allow for treatment of incremental rotations as vectorial quantities. This is similar to load step requirements as discussed for the static equilibrium iteration in Incremental Equilibrium Iterations.

Note that the equilibrium is defined at time \(\mathrm {t+\Delta \tau}\) and thus not satisfied at the end of the time step \(\mathrm {t+\Delta t}\) for \(\mathrm {\theta >1}\).

It should also be noted that hydrodynamic loading initially is based on the position and velocity vectors at time \(\mathrm {t}\). Hydrodynamic loading is, however, based on updated position and velocity vectors at time \(\mathrm {t+\Delta t}\) during the equilibrium iteration.

The sequence of calculations for the nonlinear integration procedure including equilibrium iteration is summarized in the following:

A. Initial calculation

  1. Establish integration constants based on the integration parameters \(\mathrm {\beta }\), \(\mathrm {\gamma }\) and \(\mathrm {\theta }\).

  2. Establish initial conditions \(\\\boldsymbol{r}_0=\boldsymbol{r}_{\mathrm {static}}\\\boldsymbol{\dot {r}}_0=0\\\boldsymbol{\ddot {r}}_0=0\)

B. Each time step

  1. Calculate the effective stiffness matrix, \(\mathrm {\boldsymbol{\hat{K}}_t}\), according to Equation (14).

  2. Calculate the effective load vector, \(\mathrm {\Delta \boldsymbol{\hat{R}}_t}\), according to Equation (15).

  3. Compute the incremental displacement, \(\mathrm {\Delta \boldsymbol{r}_t}\), according to Equation (13).

  4. Calculate velocity,\(\mathrm {\boldsymbol{\dot {r}}_{t+\Delta \tau}}\), and acceleration, \(\mathrm {\boldsymbol{\ddot {r}}_{t+\Delta \tau}}\) according to Equation (6).

  5. Perform the equilibrium iteration

    1. Set: \(\mathrm {\\\\^0\Delta \boldsymbol{r}_t=\Delta \boldsymbol{r}_t\\\\^0\boldsymbol{\dot {r}}_{t+\Delta \tau}=\boldsymbol{\dot {r}}_{t+\Delta \tau}\\\\^0\boldsymbol{\ddot {r}}_{t+\Delta \tau}=\boldsymbol{\ddot {r}}_{t+\Delta \tau}\\\\^0\boldsymbol{r}_{t+\Delta \tau}=\boldsymbol{r}_{t}+{^0}\Delta \boldsymbol{r}_t\\}\)

    2. Establish the effective stiffness matrix, \(\mathrm {\boldsymbol{\hat{K}}}\), based on the tangential mass, damping and stiffness matrices at iteration \(\mathrm {i-1}\). Left hand side of Equation (20).

    3. Calculate the effective residual load vector. Right hand side of Equation (20), where inertia and damping contributions are given by Equation (21).

    4. Compute the additional displacement increments, \(\mathrm {\boldsymbol{^i\Delta _r}}\) Equation (20).

    5. Calculate the improved displacement increment, \(\mathrm {^i\Delta \boldsymbol{r}_t}\) Equation (17).

    6. Calculate velocity, \(\mathrm {^i\boldsymbol{\dot {r}}_{t+\Delta \tau}}\), and acceleration, \(\mathrm {^i\boldsymbol{\ddot {r}}_{t+\Delta \tau}}\), according to Equation (19) and displacement vector, \(\mathrm {^i\boldsymbol{r}_{t+\Delta \tau}}\), according to Equation (18).

    7. Test for convergence. If the test is not fulfilled go to b.

C. Conditions at time \(\mathrm {t+\Delta t}\) ( irrelevant for \(\mathrm {\theta =1}\) )

  • \(\mathrm {\boldsymbol{\ddot {r}}_{t+\Delta t}}\) according to Equation (7).

  • \(\mathrm {\boldsymbol{\dot {r}}_{t+\Delta t}}\) and \(\mathrm {\boldsymbol{r}_{t+\Delta t}}\) according to Equation (5) where \(\mathrm {\Delta \tau=\Delta t}\).

4. Linearized Analysis

The linear time domain analysis is based on linearization of the dynamic equilibrium equations with regard to stiffness, damping and inertia forces at the static equilibrium position according to Equation (22).

\[\boldsymbol{M}\boldsymbol{\ddot {r}}+\boldsymbol{C}\boldsymbol{\dot {r}}+\boldsymbol{K}\boldsymbol{r}=\boldsymbol{R}^E(\boldsymbol{r},\boldsymbol{\dot {r}},\boldsymbol{t})\]

where

  • \(\mathrm {\boldsymbol{M}}\), \(\mathrm {\boldsymbol{C}}\), \(\mathrm {\boldsymbol{K}\qquad\qquad\quad }\) - tangential mass, damping and stiffness matrices evaluated at the static equilibrium position

  • \(\boldsymbol{r}=\boldsymbol{r}_\mathrm {tot}-\boldsymbol{r}_\mathrm {stat}\qquad\:\:\) - dynamic displacement vector relative to the static position

  • \(R^E=R_{tot}^E-R_{stat}^E\) - dynamic load vector expressing the difference between the total external load vector and the external load vector at static equilibrium

By rewriting Equation (22) at time \(\mathrm {t+\Delta \tau}\) and inserting Equation (6), the dynamic equation for displacement over the time interval \(\mathrm {\Delta \tau=\theta \Delta t}\) is written

\[\boldsymbol{\hat{K}}\boldsymbol{r}_{t+\Delta \tau}=\boldsymbol{\hat{R}}_{t+\Delta \tau}\]

where the effective stiffness matrix, \(\mathrm {\boldsymbol{\hat{K}}}\), and the effective dynamic load vector, \(\mathrm {\boldsymbol{\hat{R}}}\), are defined by the relations

\[\boldsymbol{\hat{K}}=\frac{1}{\beta (\Delta \tau)^2}\boldsymbol{M}+\frac{\gamma }{\beta \Delta \tau}\boldsymbol{C}+\boldsymbol{K}\]
\[\begin{array}{l}\displaystyle \boldsymbol{\hat{R}}_{t+\Delta \tau}=\boldsymbol{R}^E_{t+\Delta \tau}+\boldsymbol{M}(\frac{1}{\beta \Delta \tau^2}\boldsymbol{r}_t+\frac{1}{\beta \Delta \tau}\boldsymbol{\dot {r}}_t+(\frac{1}{2\beta }-1)\boldsymbol{\ddot {r}}_t)\\\\\displaystyle +\:\boldsymbol{C}(\frac{\gamma }{\beta \Delta \tau}\boldsymbol{r}_t+(\frac{\gamma }{\beta }-1)\boldsymbol{\dot {r}}_t+(\frac{\gamma }{2\beta }-1)\Delta \tau\boldsymbol{\ddot {r}}_t)\end{array}\]

The advantage of this formulation is that the effective stiffness matrix is constant throughout the integration. The effective stiffness matrix is therefore triangulated prior to the integration, and the dynamic displacement vector is found by a simple back substitution at each time step.

The external loads may be nonlinear. The external load vector is computed at the static equilibrium position. This means that normally the only nonlinearity in the external loads is the drag load expressed by quadratic relative velocity.

The coupling between external loads and structural velocity can be accounted for by an iteration procedure. In most cases, however, sufficient accuracy is obtained by using the structural velocity from the previous step.

Alternatively, the drag load calculation may be based on the following approximation for the structural velocity:

\[\boldsymbol{\dot {r}}_{t+\Delta \tau}=\boldsymbol{\dot {r}}_t+\Delta \tau\cdot \boldsymbol{\ddot {r}}_t\]

Load iteration+

During the load iteration, the external loads are updated due to improved estimates of structural velocities. All other contributions to the effective load vector are kept constant.

The convergence criterion is related to the change in magnitude of dynamic load components during the iteration according to:

\[\frac{\mathrm {max}\Big|\,^{i-1}_j\boldsymbol{\hat{R}}+\,^i_j\boldsymbol{\hat{R}}\Big|}{\mathrm {max}\Big|\,^i_k\boldsymbol{\hat{R}}\Big|}<\varepsilon\]

where the upper left indices, \(\mathrm {^i}\) and \(\mathrm {^{i-1}}\), indicate the iteration number, and the lower left indices, \(\mathrm {_j}\) and \(\mathrm {_k}\), indicate the degree of freedom. Only force components are considered, moments are not included.

The test, Equation (27), expresses that the maximum correction of a force component during the iteration is compared to the maximum force component at the same time instant. As indicated they are not necessarily identified at the same degree of freedom.

The sequence of calculations for the linearized integration procedure including load iterations is summarized in the following:

A. Initial calculation

  1. Establish integration constants based on the integration parameters \(\mathrm {\beta }\), \(\mathrm {\gamma }\) and \(\mathrm {\theta }\). For details, see Langen and Sigbjørnsson (1979).

  2. Establish initial conditions \(\mathrm {\\\boldsymbol{r}_0=0\\\boldsymbol{\dot {r}}_0=0\\\boldsymbol{\ddot {r}}_0=0}\)

  3. Calculate the effective stiffness matrix, \(\mathrm {\boldsymbol{\hat{K}}}\), based on the static equilibrium position according to Equation (24).

  4. Triangularize the effective stiffness matrix.

B. Each time step

  1. Calculate the effective dynamic load vector, \(\mathrm {\Delta \boldsymbol{\hat{R}}_{t+\Delta \tau}}\), according to Equation (25).

  2. Compute the dynamic displacement, \(\mathrm {\Delta \boldsymbol{r}_{t+\Delta \tau}}\), according to Equation (23).

  3. Calculate velocity, \(\mathrm {\Delta \boldsymbol{\dot {r}}_{t+\Delta \tau}}\), and acceleration, \(\mathrm {\Delta \boldsymbol{\ddot {r}}_{t+\Delta \tau}}\), based on Equation (6).

  4. Iterate on external loads

    1. Set: \(\mathrm {\\\\^0\boldsymbol{\hat{R}}_{t+\Delta \tau}=\Delta \boldsymbol{\hat{R}}_{t+\Delta \tau}\\\\^0\boldsymbol{r}_{t+\Delta \tau}=\boldsymbol{r}_{t+\Delta \tau}\\\\^0\boldsymbol{\dot {r}}_{t+\Delta \tau}=\boldsymbol{\dot {r}}_{t+\Delta \tau}\\\\^0\boldsymbol{\ddot {r}}_{t+\Delta \tau}=\boldsymbol{\ddot {r}}_{t+\Delta \tau}\\}\)

    2. Calculate the effective load vector, \(\mathrm {^i\boldsymbol{\hat{R}}_{t+\Delta \tau}}\), where the contribution from external loads is based on velocity estimates from iteration number \(\mathrm {i-1}\). All other contributions are kept constant.

    3. Compute the dynamic displacement, \(\mathrm {^i\boldsymbol{r}_{t+\Delta \tau}}\), Equation (23).

    4. Calculate velocity, \(\mathrm {^i\boldsymbol{\dot {r}}_{t+\Delta \tau}}\), and acceleration, \(\mathrm {^i\boldsymbol{\ddot {r}}_{t+\Delta \tau}}\) based on Equation (6).

    5. Test for convergence. If the test is not fulfilled go to b.

C. Conditions at time \(\mathrm {t+\Delta \tau}\) ( irrelevant for \(\mathrm {\theta =1}\) )

Calculate displacements, velocities and acceleration vectors at \(\mathrm {t+\Delta t}\). Calculate \(\mathrm {\boldsymbol{\ddot {r}}_{t+\Delta t}}\) according to Equation (7), and calculate \(\mathrm {\boldsymbol{\dot {r}}_{t+\Delta t}}\) and \(\mathrm {\boldsymbol{r}_{t+\Delta t}}\) according to Equation (5) where \(\mathrm {\Delta \tau=\Delta t}\).

5. Forced Vessel Rotations

Vessel rotations are irrelevant when a moment-free vessel attachment is used (i.e. bar elements or ball-joint connector used at vessel termination). Correct modelling of vessel rotations is, however, crucial for description of moments and curvatures close to the termination when the vessel attachment is a clamped beam element.

Vessel displacements and rotations are according to Forced Vessel Motion. >> described as linear responses to waves. This means that vessel rotations are treated as vector quantities in the vessel motion analysis.

The large rotation formulation applied in the finite element formulations is linearized at the static equilibrium position in linearized dynamic analysis. In the approach, it is therefore consistent to apply vessel rotations directly as prescribed rotations, i.e. 3D rotations are treated as a vector.

Special considerations are, however, required to obtain consistency between the linear vessel response model and the co-rotated finite element formulations applied in nonlinear dynamic analysis.

Nodal rotations are in the nonlinear finite element formulation treated as an Euler rotation sequence in the order \(\mathrm {\theta _1}\), \(\mathrm {\theta _2}\), \(\mathrm {\theta _3}\). Following this convention, the orientation of the nodal coordinate system relative to the initial position is described by a rotation matrix \(\mathrm {\boldsymbol{T}(\theta _1,\theta _2,\theta _3)}\) where the elements of \(\mathrm {\boldsymbol{T}}\) are given in Updating of the Rotation Matrix.

Adopting the same conventions for vessel rotations, the relation between vessel rotations at time \(\mathrm {t_k}\) and \(\mathrm {t_{k+1}}\) can be established from the equation for \(\boldsymbol{\bar{T}}_{ij}^{(n+1)}\):

\[\boldsymbol{T}(\theta ^{k+1}_1,\theta ^{k+1}_2,\theta ^{k+1}_3)=\boldsymbol{T}(\theta ^k_1,\theta ^k_2,\theta ^k_3)\boldsymbol{T}(\Delta \theta ^k_1,\Delta \theta ^k_2,\Delta \theta ^k_3)\]

where

  • \(\mathrm {\theta ^k_1}\), \(\mathrm {\theta ^k_2}\), \(\mathrm {\theta ^k_3\quad \quad \quad \quad }\) - Vessel rotations at time \(\mathrm {t_k}\)

  • \(\mathrm {\theta ^{k+1}_1}\), \(\mathrm {\theta ^{k+1}_2}\), \(\mathrm {\theta ^{k+1}_3\quad \:}\) - Vessel rotations at time \(\mathrm {t_{k+1}}\)

  • \(\mathrm {\Delta \theta ^k_1}\), \(\mathrm {\Delta \theta ^k_2}\), \(\mathrm {\Delta \theta ^k_3\:\:\quad }\) - Incremental vessel rotations from time \(\mathrm {t_k}\) to \(\mathrm {t_{k+1}}\)

Equation (28) can therefore be used to establish the incremental vessel rotations from time \(\mathrm {t_k}\) to \(\mathrm {t_{k+1}}\) that are consistent with the finite element formulation. Taking advantage of the orthogonality of the rotation matrices, the incremental rotation matrix can be expressed as

\[\boldsymbol{T}(\Delta \theta ^k_1,\Delta \theta ^k_2,\Delta \theta ^k_3)=\boldsymbol{T}^T(\theta ^k_1,\theta ^k_2,\theta ^k_3)\boldsymbol{T}(\theta ^{k+1}_1,\theta ^{k+1}_2,\theta ^{k+1}_3)\]

Using the equation for \(\boldsymbol{\tilde{T}}_{kj}\), explicit expressions for the incremental rotations are given by

\[\Delta \theta ^k_1=\sin^{-1}\bigg[-\frac{T_{32}}{\cos(\Delta \theta ^k_2\:)}\bigg]\]
\[\Delta \theta ^k_2=\sin^{-1}(T_{31})\]
\[\Delta \theta ^k_3=\sin^{-1}\bigg[-\frac{T_{21}}{\cos(\Delta \theta ^k_2\:)}\bigg]\]

where \(\mathrm {T_{ij}}\) denotes the elements of the incremental rotation matrix computed by Equation (29).

The described formulation is crucial for correct representation of simultaneous forced rotations around 2 or 3 axes in nonlinear dynamic analysis for commonly used support vessels. A prescribed rotation will, however, be identical for linear and nonlinear dynamic analyses involving forced rotation about one axis only, e.g. in-plane excitation.

6. Structural Damping Models

6.1. Global Rayleigh Damping Model

In the global Rayleigh damping formation, the tangential damping matrix is established as a linear combination of the global tangential mass- and stiffness matrices:

\[\boldsymbol{C}=\alpha _1\boldsymbol{M}_t+\alpha _2\boldsymbol{K}_t\]

\(\mathrm {\alpha _1}\) and \(\mathrm {\alpha _2}\) are denoted the mass- and stiffness-proportional damping coefficients, respectively. The motivation for introducing this approach is mainly due to computational convenience regarding implementation and specification of the structural damping level.

The most important property of the global Rayleigh damping is that the structural damping matrix as defined by Equation (33) is orthogonal with respect to the eigenvectors. The orthogonality can be used to express the modal damping for a linear dynamic system as a function of the damping coefficients

\[\lambda_i=\frac{1}{2}\bigg[\frac{\alpha _1}{\omega _i}+\alpha _2\omega _i\bigg]\]

where \(\mathrm {\lambda_i}\) is the modal damping ratio relative to critical and \(\mathrm {\omega _i}\) is the eigenfrequency. Equation (34) can also be used as a reference for specification of Rayleigh damping for nonlinear dynamic systems.

It should be noted that the damping coefficients \(\mathrm {\alpha _1}\) and \(\mathrm {\alpha _2}\) apply to all global degrees of freedom. Hence, it is impossible to specify different damping levels related to axial, bending and torsional deformations. Application of the global Rayleigh damping model is therefore restricted to description of a realistic overall damping level for the structure.

The global Rayleigh damping formulation has been widely used to model structural damping of fixed structures. Estimation of modal damping ratios at two eigenfrequencies can be used to adapt the empirical damping model to an actual structure.

The global Rayleigh damping model must, however, be used more carefully for compliant structures undergoing large rigid body displacements. The mass-proportional damping is usually omitted for such structures to avoid unphysical structural damping due to rigid body motions (i.e. \(\mathrm {\alpha _1=0}\)) (Bech and Skallerud, 1992).

The stiffness-proportional damping, \(\mathrm {\alpha _2\boldsymbol{K}}\), will give a modal damping ratio, \(\mathrm {\lambda_i=\displaystyle \frac{\alpha _2\omega _i}{2}}\), which is seen to be proportional to the eigenfrequency. In practical applications, \(\mathrm {\alpha _2}\) is selected to give a realistic energy dissipation at the peak period of the loading, i.e. peak period of wave spectrum or support vessel motion spectrum in heave.

The Rayleigh damping model can be applied in nonlinear as well as linearized dynamic analysis. An updated or constant structural damping matrix can alternatively be applied in nonlinear dynamic analyses. These alternatives are detailed in the following:

6.1.1. Updated Rayleigh damping in nonlinear dynamic analysis

The updated structural damping matrix is computed according to Equation (33) using the instantaneous tangential mass and stiffness matrices.

This approach will give a nonlinear displacement-dependent structural damping matrix. The relative importance of the additional nonlinearities will increase with the damping coefficients. The nonlinear time domain analysis tends to become more numerically sensitive with increasing structural damping when an updated damping matrix is used.

Instability problems have also been observed when using an updated stiffness-proportional damping, i.e. \(\mathrm {\alpha _1=0}\), for so-called low tension problems (i.e. analysis of compliant slender structures involving zero or close to zero axial force in parts of the structure). This is because the contribution from the geometric stiffness is insignificant for an axial force close to zero, resulting in a very low structural damping related to transverse displacements. Hence, the structural damping is `lost' when it is most needed to obtain a stable numerical solution.

From this discussion it can be concluded that use of updated Rayleigh damping mainly is applicable for description of an overall moderate damping level for numerically well-behaved systems.

6.1.2. Constant Rayleigh damping in nonlinear dynamic analysis

In this approach, the local Rayleigh damping matrices for all elements evaluated at the static equilibrium position are kept constant throughout the analysis. The local Rayleigh element damping matrix, \(\mathrm {\boldsymbol{c}}\), is in accordance with Equation (35) defined as a linear combination of the local mass, \(\mathrm {\boldsymbol{m}}\), and stiffness, \(\mathrm {\boldsymbol{k}}\), element matrices:

\[\boldsymbol{c}=\alpha _1\boldsymbol{m}+\alpha _2\boldsymbol{k}\]

The global Rayleigh structural damping is then at each time step found by assembly of the local damping matrices transformed to global system according to the instantaneous element orientation. The described approach gives a more robust formulation, and the shortcomings experienced when using an updated Rayleigh damping are to a large extent eliminated.

6.2. Separated Rayleigh Damping Model

A separated Rayleigh damping model which allows for specification of separate mass- and stiffness-proportional damping coefficients for tension, bending and torsion has been proposed by Bech et al. (1992).

In this approach, the separated coefficients are applied to the relevant degrees of freedom in the local mass and stiffness matrices before transformation to global degrees of freedom and assembly. The element damping matrix can be expressed as:

\[\boldsymbol{c}=\alpha _{1t}\boldsymbol{m}_t+\alpha _{1t0}\boldsymbol{m}_{t0}+\alpha _{1b}\boldsymbol{m}_b+\alpha _{2t}\boldsymbol{k}_t+\alpha _{2t0}\boldsymbol{k}_{t0}+\alpha _{2b}\boldsymbol{k}_b\]

where subscripts \(\mathrm {_t}\), \(\mathrm {_{t0}}\) and \(\mathrm {_b}\) refer to tension, torsion and bending contributions, respectively, and the matrices, \(\mathrm {\boldsymbol{c}}\), \(\mathrm {\boldsymbol{m}}\) and \(\mathrm {\boldsymbol{k}}\), are local element matrices. For instance, \(\mathrm {\boldsymbol{k}_b}\) includes all bending deformation terms in the local element stiffness matrix.

Separate viscous damping levels for tension, bending and torsion can be directly introduced in the numerical model by a proper selection of mass and stiffness coefficients for tension, bending and torsion, respectively. Use of mass-proportional damping in torsion should still be avoided to prevent damping due to rigid body motion.

This approach will give a global damping matrix that in general is not orthogonal with respect to the eigenvectors. It has, however, been experienced that Equation (33) still can be used to specify the damping independently for tension, torsion and bending with sufficient accuracy (Hanson et al., 1994).

The separated model can be used in nonlinear as well as linearized dynamic analysis. An updated or constant structural damping matrix can be used in nonlinear analyses as described in Section 6.1.

6.3. Combined Use of Rayleigh Damping and Hysteretic Models

Typical cross-sectional properties of unbonded flexible pipes can be summarized as:

  1. Different structural damping levels can be related to axial, torsional and bending deformations.

  2. Axial and torsional stiffness can be regarded as constant.

  3. Low viscous damping is related to axial and torsional deformations.

  4. Significant hysteretic bending moment/curvature behaviour for large amplitude bending vibrations.

  5. Low viscous damping and high (constant) bending stiffness for small amplitude bending vibrations.

A realistic modelling of such cross-sectional behaviour can be obtained by combining the global or separated Rayleigh damping model with a hysteretic moment/curvature description.

The hysteretic model will give an amplitude dependent damping and bending stiffness for large amplitude bending vibrations, and a high constant bending stiffness and no damping for small amplitude vibrations.

The Rayleigh model accounts for the viscous structural damping related to axial, torsional and small amplitude bending vibrations.

7. Guidance to Dynamic Analysis

7.1. Nonlinear Dynamic Analysis

7.1.1. Applications

The nonlinear dynamic analysis should be applied when at least one of the following nonlinear phenomena is of significance for the dynamic response under consideration:

  • Nonlinear material properties.

  • Time-dependent geometric stiffness variation caused by significant tension and geometry variations. This is especially important in analyses of low tension problems involving snapping and possible compression in parts of the structure. Can also be important for prediction of angular deflections, curvatures, and moments close to supports in case of large transverse displacements.

  • Time-dependent location of seafloor touchdown point. This can be of special importance for analysis of anchor lines exposed to first and second order vessel motion excitation. It can also be important for prediction of curvature, moments, and transverse displacements in touchdown area for flexible risers, umbilicals and pipelines under installation.

  • Elastic contact problems such as vessel/riser impact, collision between risers and pipeline/stinger interaction.

  • Integration to actual surface elevation. Nonlinear analysis must be considered if variation in added mass is of significance.

  • The external load is in nonlinear analysis always applied at the instantaneous structural position. Nonlinear analysis should therefore also be considered when using wave kinematics computed at instantaneous structural position to obtain full consistency in the numerical model.

  • Structures undergoing large dynamic rotations in 3D space. Nonlinear analysis is required to give an adequate representation of moments and curvatures.

7.1.2. Nonlinear analysis is mandatory in the following situations:

  • Structures including partly submerged floating elements. Examples of such structures are floating flexible fish farm systems and floating/partly floating offshore loading hoses.

  • Release/rupture analysis. Such analyses are typically carried out for anchor lines and flexible riser systems.

  • Dynamic analyses including modelling of internal slug flow.

7.1.3. Numerical solution procedures

Use of constant average acceleration over each time step (\(\mathrm {\gamma =1/2}\), \(\mathrm {\beta =1/4}\), \(\mathrm {\theta =1}\)) is recommended for nonlinear dynamic analysis. It is further recommended to use true Newton-Raphson iteration with convergence criteria similar to the static situation, see Guidance to Static Analysis.

The constant Rayleigh damping formulation will contribute to a more numerically well-behaved system and should therefore be applied whenever possible.

The time step size required to obtain a stable numerical solution is strongly system- and excitation dependent. The time step should in general be sufficiently small to give an adequate representation of the external loading and to give stable integration of all eigenperiods of the discretized system model.

For structural analyses of moderately nonlinear systems it has been experienced that using 70-200 time steps per average load period normally will be sufficient. A finer time resolution is required for systems with more pronounced nonlinear characteristics.

Examples of more numerically sensitive analyses requiring a rather small time step and narrow convergence criterion are:

  • Systems undergoing large transverse seafloor excursions.

  • Low tension problems including snapping and compression.

  • Nonlinear cross sectional modelling using bending moment/curvature hysteresis.

  • Transient release/rupture analysis.

7.2. Linearized Dynamic Analysis

7.2.1. Application

The linearized dynamic analysis is a numerically robust and efficient formulation which can be applied with sufficient accuracy for a wide range of systems. The linearized analysis is in particular attractive in cases with moderate structural nonlinearity, where the external hydrodynamic loading is the main source of nonlinearity. Examples of such structures are standard flexible riser configurations and anchor lines in normal operating conditions.

The computational time consumption in linearized analysis is typically 1/10 of nonlinear analysis. The linearized approach should therefore be used whenever possible to reduce the computational efforts. Nonlinear verification analysis of a few, critical cases should, however, always be considered when using linearized analysis.

7.2.2. Numerical solution procedures

Use of constant acceleration over each time step (\(\mathrm {\gamma =1/2}\), \(\mathrm {\beta =1/4}\), \(\mathrm {\theta =1}\)) or the Wilson \(\mathrm {\theta }\)-method (\(\mathrm {\gamma =1/2}\), \(\mathrm {\beta =1/6}\), \(\mathrm {\theta =1.4}\)) is recommended in linearized analysis.

The external loading normally is well predicted using the structural velocity from previous time step. Load iteration is therefore in general not required.

The number of time steps per average load period is typically 50-150 for standard linearized analysis.

8. References

Bech, A. and Skallerud, B. (1992): Structural Damping in Flexible Pipes: Comparisons between Dynamic Tests and Numerical Simulations, 2nd International Offshore and Polar Engineering Conference, San Francisco.

Bech, A., Skallerud, B. and Sødahl, N. (1992): Structural Damping in Design Analysis of Flexible Risers, Marinflex 92, London.

Clough, W. and Penzien, J. (1975): Dynamics of Structures, McGraw-Hill.

Hanson, T.D., Otteren, A. and Sødahl, N. (1994): Response Calculation Using an Enhanced Model for Structural Damping in Flexible Risers Compared with Full Scale Measurements, Proceedings of the international conference on hydroelasticity in Marine Technology, Trondheim.

Langen, I. and Sigbjørnsson, R. (1979): Dynamisk analyse av konstruksjoner, Tapir, Trondheim.

Remseth, S.N. (1978): Nonlinear Static and Dynamic Analysis of Space Structures, Dr.ing. Thesis, Report No. 78-2, Division of Structural Mechanics, The Norwegian Institute of Technology, Trondheim.

Sødahl, N. (1991): Methods for Design and Analysis of Flexible Riser Systems, Dr.ing. Thesis, Div. of Marine Structures, The Norwegian Institute of Technology, Trondheim.

Sødahl, N. and Larsen, C.M. (1992): Methods for Prediction of Extreme Responses of Flexible Risers, 2nd International Offshore and Polar Engineering Conference, San Francisco.