Finite Element Formulation

1. Introduction

The structural analysis part of the RIFLEX program is based on finite element modelling. This chapter presents the theoretical basis for the finite element formulation. The current formulation is with few exceptions equal to the one used in FENRIS, which is a large scale, general purpose program for nonlinear finite element analysis. A detailed description of the theoretical basis for FENRIS has been given by Bergan et al. (1984, 1985), Mollestad (1983) and Engseth (1984).

The most important features included in the basic finite element formulation are:

  • Unlimited rotations and translations in 3D space.

  • Beam and bar element based on small strain theory.

  • Use of beam as well as bar elements in system modelling.

  • Connector element for modelling of swivels, hinges, etc.

  • Description of nonlinear material properties.

  • Stiffness contribution from material properties as well as geometric stiffness, i.e. contribution of axial force to the transverse stiffness.

  • General element assembly allowing for description of arbitrary system topology, varying cross-sectional properties, arbitrary boundary conditions at supports, prescribed support displacements and specified external forces.

The described features are crucial for adequate modelling and analysis of slender marine structures in general. Application of the basic finite element formulation in static and dynamic analysis are detailed in Static Finite Element Analysis and Dynamic Time Domain Analysis, respectively.

A short review of basic continuum mechanics theory is given in Section 2. Using generalized strain/stress measures allows the description of large cross-sectional deformations. Implementation of the finite element method including description of translation and rotations in 3D space is described in Section 3. All element formulations are based on the small strain approximation which is reasonable for application to slender marine structures.

2. Fundamental Continuum Mechanics Theory

2.1. Description of Motion

The theoretical foundation is based on well established principles of continuum mechanics. Details of the formulation may be found in the textbook by Malvern (1969).

The Lagrangian description is used to describe the motion of material particles. This motion is referred to a fixed global Cartesian coordinate frame defined by the base vectors \(\mathrm {\boldsymbol{I}_i}\) (see Figure 1). Each particle is identified by its position vector \(\mathrm {\boldsymbol{X}}\) in an arbitrarily chosen reference configuration, which may be the initial configuration \(\mathrm {C_0}\) as shown in Figure 1. The motion of the particle can be expressed as

\[\boldsymbol{x}=\boldsymbol{x}(\boldsymbol{X},t)\]

where \(\mathrm {\boldsymbol{x}}\) denotes the position of the particle at the time \(\mathrm {t}\). The displacement vector \(\mathrm {\boldsymbol{u}}\) is defined through

\[\boldsymbol{x}=\boldsymbol{X}+\boldsymbol{u}.\]

In a nonlinear analysis it is necessary to consider the initial configuration \(\mathrm {C_{0}}\) of the body, the deformed configuration \(\mathrm {C_n}\) at a given time \(\mathrm {t}\) and a new incremental configuration \(\mathrm {C_{n+1}}\) for time \(\mathrm {t+\Delta t}\). This is illustrated in Figure 2. The new configuration \(\mathrm {C_{n+1}}\) is close to \(\mathrm {C_n}\) since \(\mathrm {\Delta t}\) is assumed to be small.

A formulation where the strains in \(\mathrm {C_n}\), \(\mathrm {C_{n+1}}\), and all other configurations are referred to the initial configuration \(\mathrm {C_0}\), is usually termed a total Lagrangian formulation. However, in some cases it is convenient to relate the motion of material particles within a body to a local, rectangular coordinate frame which translates and rotates along with the average motion of the body. The total motion is then found by combining the motion of the local position vector and the motion of the local reference system. An example of this is the concept of co-rotated ghost reference, introduced by Bergan (1985). This formulation assumes that the initial configuration translates and rotates as a rigid body so that it at any time is located close to the actual deformed configuration. In Figure 2 this is illustrated by \(\mathrm {C_n}\) and \(\mathrm {C_{0n}}\), where the original base vectors \(\mathrm {\boldsymbol{I}_i}\) have rotated along with \(\mathrm {C_{0n}}\) and are denoted \(\mathrm {\boldsymbol{\tilde{I}}_i}\)

RIFLEX TheoryManual 42 v0 23
Figure 1. Motion of a material particle in space
RIFLEX TheoryManual 42 v0 24
Figure 2. Various configurations of a body.

In RIFLEX the bar element is formulated using a total Lagrangian description, while the beam element formulation uses a co-rotated ghost reference description. This will be described in more detail in Section 3.4 and in Section 3.

2.2. Definitions of Strain and Stresses

For the Lagrangian formulation the strains are measured in terms of the Green strain tensor \(\mathrm {\boldsymbol{E}}\). If \(\mathrm {{C}_0}\) is used as reference configuration, this strain tensor is defined by

\[ \mathrm {d}S^2_n-\mathrm {d}S^2_0=2\mathrm {d}\boldsymbol{X}\cdot \boldsymbol{E}\cdot \mathrm {d}\boldsymbol{X}\]

where \(\mathrm {d}S_0\) and \(\mathrm {d}S_n\) denote the length of the line segment PQ before and after deformation, (see Figure 2). If \(\mathrm {C_{0n}}\) is used as reference configuration, the definition of the Green strain tensor reads

\[ \mathrm {d}S^2_n-\mathrm {d}S^2_{0n}=2\mathrm {d}\boldsymbol{\tilde{X}}\cdot \boldsymbol{\tilde{E}}\cdot \mathrm {d}\boldsymbol{\tilde{X}}\]

The strain tensors may also be expressed as

\[\boldsymbol{E}=E_{ij}\boldsymbol{I}_{i}\boldsymbol{I}_j\]
\[\boldsymbol{\tilde{E}}=\tilde{E}_{ij}\boldsymbol{I}_{i}\boldsymbol{I}_j=\overset{\approx }{E}_{ij}\boldsymbol{\tilde{I}}_{i}\boldsymbol{\tilde{I}}_j\]

where \(\mathrm {E_{ij}}\) and \(\mathrm {\tilde{E}_{ij}}\) denote the components of\(\mathrm {\boldsymbol{E}}\) and\(\mathrm {\boldsymbol{\tilde{E}}}\) in the original coordinate frame \(\mathrm {\boldsymbol{I}_j}\). \(\mathrm {\tilde{E}_{ij}}\) denotes the components of\(\mathrm {\boldsymbol{\tilde{E}}}\) in the co-rotated coordinate frame \(\mathrm {\boldsymbol{\tilde{I}}_j}\). The strain tensor \(\mathrm {\boldsymbol{\tilde{E}}}\) is generally different from \(\mathrm {\boldsymbol{E}}\). However, it can be shown that the following identity holds between the components

\[E_{ij}=\overset{\approx }{E}_{ij}\]

Equation (7) signifies that the strain components remain unchanged when both the reference body and the base vectors are rotated in the same way. No tensor transformations of the components are necessary. This result is very useful and makes this formulation especially attractive.

The rectangular components of \(\mathrm {E}\) referred to \(\mathrm {\boldsymbol{I_i}}\) and \(\mathrm {\boldsymbol{I_j}}\), may be expressed as

\[E_{ij}=\frac{1}{2}(\frac{\partial u_i}{\partial X_j}+\frac{\partial u_j}{\partial X_i}+\frac{\partial u_m}{\partial X_i}\frac{\partial u_m}{\partial X_j})\]

where the components of the displacement vector\(\mathrm {\boldsymbol{u}}\) have been introduced. Note that the differentiation in Equation (8) is performed with respect to the material coordinates \(\boldsymbol{\mathrm {X_i}}\). It may also be seen that\(\mathrm {\boldsymbol{E}}\) is a symmetric tensor consisting of both linear and quadratic terms. The quadratic contributions are only of importance for large displacements and stability problems.

An appropriate stress measure is also needed in the analysis. The symmetric Piola-Kirchhoff stress tensor \(\mathrm {\boldsymbol{S}}\) is always used in conjunction with Green strain. A detailed description of this stress measure is omitted here, but may be found in Malvern (1969).

The symmetric Piola-Kirchhoff stress tensor referred to\(\mathrm {C_0}\) may be expressed as

\[\boldsymbol{S}=S_{ij}\boldsymbol{I}_{i}\boldsymbol{I}_j\]

Similar to what was done for the strains, this stress tensor may also be referred to a co-rotated ghost reference \(\mathrm {C_{0n}}\). In this case the stress tensor may be written as

\[\boldsymbol{\tilde{S}}=\tilde{S}_{ij}\boldsymbol{I}_{i}\boldsymbol{I}_j=\overset{\approx }{S}_{ij}\boldsymbol{\tilde{I}}_{i}\boldsymbol{\tilde{I}}_j\]

As for the strains, it may be shown that

\[S_{ij}=\overset{\approx }{S}_{ij}\]

which means that no transformation of the stress components is necessary when the co-rotated ghost reference formulation is used.

2.3. Virtual Work Equations

The virtual work equation may be used to express equilibrium for a finite body. Using Green strains and Piola-Kirchhoff stresses, this equation may be written

\[\int_{V_0}\boldsymbol{S}:\delta \boldsymbol{E}\mathrm {d}V_0=\int_{A_0}\boldsymbol{t}_0\cdot \delta \boldsymbol{u}\mathrm {d}A_0+\int_{V_0}\boldsymbol{f}_0\cdot \delta \boldsymbol{u}\mathrm {d}V_0\]

where \(\mathrm {\delta }\) indicates virtual quantities, and \(\mathrm {A_0}\) and \(\mathrm {V_0}\) express the surface and volume of the initial reference configuration. The surface traction\(\mathrm {\boldsymbol{t}_0}\) and body forces\(\mathrm {\boldsymbol{f}_0}\) are referred to a unit surface and a unit volume in the initial reference state.

Most solution procedures for nonlinear problems employ a linearized incremental form of the equilibrium equations based on two neighbouring equilibrium configurations.\(\mathrm {C_n}\) and \(\mathrm {C_{n+1}}\) in Figure 2 are examples of this. The incremental equations are best expressed using an incremental form of the virtual work principle.

\[\int_{V_0}(\boldsymbol{S}:\delta \Delta \boldsymbol{E}+\Delta \boldsymbol{S}:\delta \boldsymbol{E})\mathrm {d}V_0=\int_{A_0}\Delta \boldsymbol{t}_0\cdot \delta \boldsymbol{u}\mathrm {d}A_0+\int_{V_0}\Delta \boldsymbol{f}_0\cdot \delta \boldsymbol{u}\mathrm {d}V_0\]

where the \(\mathrm {\Delta }\) is used to denote finite but small increments between\(\mathrm {C_n}\) and \(\mathrm {C_{n+1}}\).

The first term in Equation (13) which depends on the current state of stress \(\mathrm {\boldsymbol{S}}\), gives rise to the geometric stiffness matrix. The second term which depends on the incremental material law, is the basis for the material stiffness matrix. Note that both of these stiffness matrices are symmetric.

Equations Equation (12) and Equation (13) are also valid when a co-rotated ghost reference system is used. Merely replace \(\mathrm {{E}_{ij}}\) and \(\mathrm {{S}_{ij}}\) by \(\mathrm {\overset{\approx }{E}_{ij}}\) and \(\mathrm {\overset{\approx }{S}_{ij}}\).

2.4. Dynamic Equilibrium Equations

The dynamic equilibrium equation expressed in terms of virtual work can be written, see Remseth (1978)

\[\int_{V_0}\boldsymbol{S}:\delta \boldsymbol{E}\mathrm {d}V_0+\int_{V_0}\rho _0\boldsymbol{\ddot {u}}\cdot \delta \boldsymbol{u}\mathrm {d}V_0+\int_{V_0}\tilde{c}\boldsymbol{\dot {u}}\cdot \delta \boldsymbol{u}\mathrm {d}V_0=\int_{A_0}\boldsymbol{t}_0\cdot \delta \boldsymbol{u}\mathrm {d}A_0+\int_{V_0}\boldsymbol{f}_0\cdot \delta \boldsymbol{u}\mathrm {d}V_0\]

where \(\mathrm {\rho _0}\) denotes mass density and \(\mathrm {\tilde{c}}\) is a viscous damping density function, i.e. damping forces are proportional to velocity. The inertia forces are proportional to the mass and acceleration of the structure. The additional terms in Equation (14) as compared to Equation (12), give rise to the mass and the damping matrices.

The incremental form of the virtual work equation yields

\[\begin{array}{l}\displaystyle \int_{V_0}(\boldsymbol{S}:\delta \Delta \boldsymbol{E}+\Delta \boldsymbol{S}:\delta \boldsymbol{E})\mathrm {d}V_0+\int_{V_0}\rho _0\Delta \boldsymbol{\ddot {u}}\cdot \delta \boldsymbol{u}\mathrm {d}V_0+\int_{V_0}\tilde{c}\Delta \boldsymbol{\dot {u}}\cdot \delta \boldsymbol{u}\mathrm {d}V_0\\\\\displaystyle =\int_{A_0}\Delta \boldsymbol{t}_0\cdot \delta \boldsymbol{u}\mathrm {d}A_0+\int_{V_0}\Delta \boldsymbol{f}_0\cdot \delta \boldsymbol{u}\mathrm {d}V_0\end{array}\]

3. Implementation of the Finite Element Method

3.1. Large Rotation in Space

The finite element nodal points may have up to six degrees of freedom, i.e. three translations and three rotations. In the case of translational degrees of freedom only, the position of the displaced nodal point is uniquely defined by the initial position vector \(\mathrm {\boldsymbol{X}}\) and the nodal displacement vector \(\mathrm {\boldsymbol{v}}\).

\[\boldsymbol{x}=\boldsymbol{X}+\boldsymbol{v}\quad \mathrm {or}\quad x_i\boldsymbol{I}_i=(X_i+v_i)\boldsymbol{I}_i\]

The case of a node that is both translated and rotated must be treated more carefully. This is because large rotations in space are not true vectors that may be expressed by vectorial components in a base coordinate system. The current formulation uses a coordinate system with base vectors \(\mathrm {\boldsymbol{\bar{i}}_i}\) that are frozen to the nodal point and that follow the movement of the node. This coordinate system is parallel to the global base vectors, \(\mathrm {\boldsymbol{\bar{I}}_i}\) in the initial configuration as shown in Figure 3. The orientation of the nodal point in space is uniquely defined by the base vector transformation.

\[\boldsymbol{\bar{i}}_i^{(n)}=\boldsymbol{\bar{T}}_{ij}^{(n)}\boldsymbol{I}_j\]

The rotation of a node is thus not given by angles, but by the components of a rotation matrix which is orthogonal. Consequently, in the current formulation the general motion of a nodal point is described by the three displacement components \(\mathrm {v_i}\) and the nine elements of the rotation matrix \(\mathrm {\boldsymbol{\bar{T}}_{ij}}\).

During the incremental solution procedure \(\mathrm {\boldsymbol{\bar{T}}_{ij}}\) is updated by treating the incremental rotational components about the global axes as rotations for small displacement problems. This will be described in more detail in Updating of the Rotation Matrix.

RIFLEX TheoryManual 42 v0 87
Figure 3. Nodal point with translational and rotational degrees of freedom

3.2. Finite Rigid Bodies at Nodes

Finite rigid bodies may be introduced at the nodes. This feature makes it possible to accommodate for finite, rigid structural parts and to model eccentricities at nodes. This is illustrated in Figure 4, where the element is thought to be attached to the surface of the rigid body. The initial location of the surface point is defined in terms of an eccentricity vector \(\mathrm {\boldsymbol{e}^0}\). After deformation has taken place, the position of the surface point is given by its position vector \(\mathrm {\boldsymbol{x}_p}\).

\[\boldsymbol{x}_p=\boldsymbol{x}+\boldsymbol{e}=(X_i+v_i+e^0_j\boldsymbol{\bar{T}}_{ij})\boldsymbol{I}_i\]

where \(\mathrm {\boldsymbol{e}^0_j}\) are the components of the initial eccentricity vector \(\mathrm {\boldsymbol{e}^0}\) and \(\mathrm {\boldsymbol{\bar{T}}_{ij}}\) is the rotation matrix components of the nodal point defined by Equation (17).

RIFLEX TheoryManual 42 v0 94
Figure 4. Motion of a point P at the surface of a rigid body

3.3. Updating of the Rotation Matrix

The transformation expressed by Equation (18) needs to be updated during the incremental solution procedure. The problem at hand is to establish the new transformation

\[\boldsymbol{\bar{i}}_i^{(n+1)}=\boldsymbol{\bar{T}}_{ij}^{(n+1)}\boldsymbol{I}_j\]

for configuration \(\mathrm {C_{n+1}}\) when the corresponding transformation for configuration \(\mathrm {C_{n}}\) is known. Provided that the incremental rotations from state \(\mathrm {C_{n}}\) to \(\mathrm {C_{n+1}}\) are small, this rotational motion (\(\mathrm {\boldsymbol{\bar{i}}_i^{(n)}arrow\boldsymbol{\bar{i}}_i^{(n+1)}}\)) may be expressed as

\[\Delta \theta =\Delta \theta _i\boldsymbol{I}_i\]

where \(\mathrm {\Delta \theta _i}\) denotes the rotational components about the global axes \(\mathrm {\boldsymbol{I}_i}\). Bergan et al. (1985) and Mollestad (1983) show that the orientation of the nodal coordinate system in \(\mathrm {C_{n+1}}\) is given as

\[\boldsymbol{\bar{i}}_i^{(n+1)}=\boldsymbol{\bar{T}}_{ik}^{(n)}\boldsymbol{\tilde{T}}_{kj}\boldsymbol{I}_j\]

thus

\[\boldsymbol{\bar{T}}_{ij}^{(n+1)}=\boldsymbol{\bar{T}}_{ik}^{(n)}\boldsymbol{\tilde{T}}_{kj}\]

\(\mathrm {\boldsymbol{\tilde{T}}_{kj}}\) denotes the components of the rotation matrix representing the incremental rotations from \(\mathrm {C_n}\) to \(\mathrm {C_{n+1}}\), and given as

\[\boldsymbol{\tilde{T}}_{kj}=\begin{bmatrix}C2\cdot C3&\quad S3\cdot C1+S1\cdot S2\cdot C3&\quad S1\cdot S3-S2\cdot C1\cdot C3\\-S3\cdot C2&\quad C1\cdot C3-S1\cdot S2\cdot S3&\quad S1\cdot C3+S2\cdot S3\cdot C1\\S2&\quad -S1\cdot C2&\quad C1\cdot C2\end{bmatrix}\]

where

\[\begin{array}{l}S1=\sin(\Delta \theta _1)\quad C1=\cos(\Delta \theta _1)\\S2=\sin(\Delta \theta _2)\quad C2=\cos(\Delta \theta _2)\\S3=\sin(\Delta \theta _3)\quad C3=\cos(\Delta \theta _3)\end{array}\]

The above expression for \(\mathrm {\boldsymbol{\tilde{T}}_{kj}}\) is found by taking the rotations \(\mathrm {\Delta \theta _i}\) in the order \(\mathrm {\Delta \theta _1}\), \(\mathrm {\Delta \theta _2}\), \(\mathrm {\Delta \theta _3}\). By taking these rotations in another order, the result in Equation (23) will become somewhat different. However, these differences become negligible when the rotations are small. It should also be kept in mind that Equation (23) is to be used in connection with a linearized incremental formulation which assumes small increments. The most important property of the rotation matrix is that it does not introduce orthogonality errors during the sequence of incremental rotations \(\mathrm {\Delta \theta _1}\), \(\mathrm {\Delta \theta _2}\), \(\mathrm {\Delta \theta _3}\).

The updated rotation matrix in Equation (22) and the translations Equation (16) now defines the new orientation of the node in configuration \(\mathrm {C_{n+1}}\).

Finally, it should be noted from Equation (23) that the following relations hold when the rotations are small

\[\begin{array}{l}\Delta \theta _1\approx \tilde{T}_{23}\\\Delta \theta _2\approx \tilde{T}_{31}\\\Delta \theta _3\approx \tilde{T}_{12}\end{array}\]

3.4. Bar Element

The spatial bar element is described in a total Lagrangian formulation. The description is based on Bergan (1984) and adjusted to a formulation based on integrated cross-section forces and small strain theory.

RIFLEX TheoryManual 42 v0 125
Figure 5. Bar element in initial and deformed configuration.

The element is assumed to be straight with an initial cross-sectional area \(\mathrm {A_0}\) which is constant along the element length. Each of the two nodes has three translational degrees of freedom, which are expressed directly in the global coordinate system (see Figure 5). The element length is denoted \(\mathrm {L_0}\) and \(\mathrm {L}\) in the initial and deformed configuration, respectively. The deformed element length is given by

\[L=\sqrt{\Delta x^2+\Delta y^2+\Delta z^2}\]

where

\(\mathrm {\Delta x=x_2-x_1}\)

\(\mathrm {\Delta y=y_2-y_1}\)

\(\mathrm {\Delta z=z_2-z_1}\)

Based on a total Lagrangian formulation and linear displacement functions, the Green strain is expressed

\[E_f=\frac{1}{2}\frac{L^2-L_0^2}{L_0^2}=\frac{1}{2L_0^2}(\Delta x^2+\Delta y^2+\Delta z^2-L_0^2)\]

The Piola-Kirchhoff stress \(\mathrm {S_f}\) is found through the constitutive law

\[S_f=S_f(E_f,E_0,S_0)\]

where \(\mathrm {E_0}\) is initial strain and \(\mathrm {S_0}\) is initial stress.

The linear displacement function (field function) for the three displacement components is given by

\[\boldsymbol{N}=\boldsymbol{N}_u=\boldsymbol{N}_v=\boldsymbol{N}_w=\begin{bmatrix}1-\xi &\quad \xi \end{bmatrix}\]

where \(\mathrm {\xi =\frac{x}{L}}\) is normalized length coordinate of the element.

In the following small strain theory is used, and it is assumed that \(\mathrm {L_0}\) is the initial stress free element length. Thus, the axial force of the element is given by

\[N=\frac{L-L_0}{L_0}(EA)\]

where \(\mathrm {L_0}\) is initial, stress free element length and \(\mathrm {EA}\) is the axial stiffness. The strain, \(\mathrm {\varepsilon }\), is given by

\[\varepsilon =\frac{L-L_0}{L_0}\]

3.4.1. Internal loads

Based on the virtual work equation, Equation (12), the internal reaction force vector is expressed by

\[\boldsymbol{S}_{int}=\begin{bmatrix}S_{x1}\\S_{y1}\\S_{z1}\\S_{x2}\\S_{y2}\\S_{z2}\\\end{bmatrix}=\frac{N}{L_0}\cdot \begin{bmatrix}-\Delta x\\-\Delta y\\-\Delta z\\\Delta x\\\Delta y\\\Delta z\end{bmatrix}\]

3.4.2. Tangential stiffness

The incremental form of the virtual work principle given by Equation (13) is used to derive the following tangential stiffness relation for the element

\[\Delta \boldsymbol{S}_{int}=\boldsymbol{k}\Delta \boldsymbol{v}=(\boldsymbol{k}_G+\boldsymbol{k}_M)\Delta \boldsymbol{v}\]

where \(\mathrm {\Delta \boldsymbol{S}_{int}}\) denotes increment in the internal load vector and \(\mathrm {\boldsymbol{k}_G}\) and \(\mathrm {\boldsymbol{k}_M}\) denote geometric and material stiffness matrices, respectively. \(\mathrm {\Delta \boldsymbol{v}}\) is the incremental displacement vector. It may be shown that \(\mathrm {\boldsymbol{k}_G}\) and \(\mathrm {\boldsymbol{k}_M}\) are given as

\[\boldsymbol{k}_G=\frac{N}{L_0}\begin{bmatrix}\boldsymbol{I}_3&-\boldsymbol{I}_3\\-\boldsymbol{I}_3&\boldsymbol{I}_3\\\end{bmatrix}\]
\[\boldsymbol{k}_M=\frac{(EA)_T}{L_0^3}\begin{bmatrix}\boldsymbol{D}^T\boldsymbol{D}&-\boldsymbol{D}^T\boldsymbol{D}\\-\boldsymbol{D}^T\boldsymbol{D}&\boldsymbol{D}^T\boldsymbol{D}\end{bmatrix}\]

where

\[\boldsymbol{D}=\begin{bmatrix}\Delta x&\Delta y&\Delta z\end{bmatrix}\]

and \(\mathrm {I_3}\) is the 3x3 unity matrix.

\(\mathrm {(EA)_T}\) is the tangential axial stiffness giving the relation between incremental strain and incremental axial force

\[\Delta N=(EA)_T\Delta \varepsilon\]

3.4.3. Distributed external loads

The consistent external nodal load vector is generally written

\[\begin{array}{l}\boldsymbol{S}_u^E=\displaystyle \int_L\boldsymbol{N}_u^Tp_x(x)\mathrm {d}x\\\\\boldsymbol{S}_v^E=\displaystyle \int_L\boldsymbol{N}_v^Tp_y(x)\mathrm {d}x\\\\\boldsymbol{S}_w^E=\displaystyle \int_L\boldsymbol{N}_w^Tp_z(x)\mathrm {d}x\end{array}\]

where \(\mathrm {p(x)}\) is distributed external loads per unit length in the deformed configuration. \(\mathrm {\boldsymbol{N}_u}\), \(\mathrm {\boldsymbol{N}_v}\) and \(\mathrm {\boldsymbol{N}_w}\) are the displacement functions given by Equation (29).

In most cases it is assumed that the load components are linearly distributed over the element. Using the relation in Equation (29), the external load vector is given by

\[\begin{array}{l}\boldsymbol{S}_u^E=\displaystyle \int_L\boldsymbol{N}^T\boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_u\\\\\boldsymbol{S}_v^E=\displaystyle \int_L\boldsymbol{N}^T\boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_v\\\\\boldsymbol{S}_w^E=\displaystyle \int_L\boldsymbol{N}^T\boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_w\\\end{array}\]

where \(\mathrm {\boldsymbol{p}_u}\), \(\mathrm {\boldsymbol{p}_v}\) and \(\mathrm {\boldsymbol{p}_w}\) are load intensity vectors at the element nodes.

By carrying out the integrals in Equation (39) the external nodal load vector is written

\[\boldsymbol{S}^E=\begin{bmatrix}S_{x1}\\S_{y1}\\S_{z1}\\S_{x2}\\S_{y2}\\S_{z2}\end{bmatrix}=\frac{L}{6}\begin{bmatrix}2p_{x1}+p_{x2}\\2p_{y1}+p_{y2}\\2p_{z1}+p_{z2}\\p_{x1}+2p_{x2}\\p_{y1}+2p_{y2}\\p_{z1}+2p_{z2}\end{bmatrix}\]

The external load vector can alternatively be calculated using a lumped load formulation, which gives

\[\boldsymbol{S}^E=\begin{bmatrix}S_{x1}\\S_{y1}\\S_{z1}\\S_{x2}\\S_{y2}\\S_{z2}\end{bmatrix}=\frac{L}{2}\begin{bmatrix}p_{x1}\\p_{y1}\\p_{z1}\\p_{x2}\\p_{y2}\\p_{z2}\end{bmatrix}\]

It should be noted that volume loads (weight and buoyancy) are calculated directly in the global reference system, while other load contributions are calculated in the local element system before transformation to global system.

For elements penetrating the sea surface, contributions to the external load vector from hydrodynamic loads and buoyancy are calculated considering the wetted element length. In this case, a linear distribution is assumed from the submerged element end to the intersection point between the element and the sea surface.

For floating partly submerged cross-sections (special option) consistent formulation is always used. The contributions from hydrodynamic loading are calculated based on numerical integration of Equation (38).

Note that hydrodynamic load contributions proportional to structural accelerations and velocities, i.e. added mass and hydrodynamic damping forces, are included in the mass and damping matrices respectively.

3.4.4. Mass matrix

The mass matrix contains contributions from structural and hydrodynamic (added) mass according to

\[\boldsymbol{m}=\boldsymbol{m}^s+\boldsymbol{m}^h\]

where \(\mathrm {\boldsymbol{m}^s}\) is the structural mass matrix and \(\mathrm {\boldsymbol{m}^h}\) is the added mass matrix.

The structural mass matrix is calculated according to

\[\boldsymbol{m}^s=m^s\int_{L_0}\boldsymbol{N}^T\boldsymbol{N}\mathrm {d}x\]

where \(\mathrm {\boldsymbol{m}^s}\) is structural mass per unit length and \(\mathrm {\boldsymbol{N}}\) is the displacement function matrix given by

\[\boldsymbol{N}=\begin{bmatrix}N_{u1}&&&&&\\&N_{v1}&&&&\\&&N_{w1}&&&\\&&&N_{u2}&&\\&&&&N_{v2}&\\&&&&&N_{w2}\end{bmatrix}\]

All terms that are not indicated in Equation (44) are equal to zero.

\(\mathrm {\boldsymbol{N}_u}\), \(\mathrm {\boldsymbol{N}_v}\) and \(\mathrm {\boldsymbol{N}_w}\) are defined in Equation (29).

The direction-dependent added mass matrix is calculated in the local system according to

\[\begin{array}{l}\boldsymbol{\boldsymbol{m}}^h_{uu}=\displaystyle \int_{L_0}\boldsymbol{N}^Tm_x^h\boldsymbol{N}\mathrm {d}x\\\\\boldsymbol{m}^h_{vv}=\displaystyle \int_{L_0}\boldsymbol{N}^Tm_y^h\boldsymbol{N}\mathrm {d}x\\\\\boldsymbol{m}^h_{ww}=\displaystyle \int_{L_0}\boldsymbol{N}^Tm_z^h\boldsymbol{N}\mathrm {d}x\end{array}\]

where \(\mathrm {m_x^h}\) , \(\mathrm {m_y^h}\) and \(\mathrm {m_z^h}\) are added mass per unit length referring to local x, y and z-axes, respectively.

For elements penetrating the sea surface, only the submerged part contributes to the added mass.

The mass matrix can alternatively be calculated using a lumped mass formulation, which in the local system is given by

\[\boldsymbol{m}=\frac{1}{2}L\cdot \begin{bmatrix}m^s+m^h_x&&&&&\\&m^s+m_y^h&&&&\\&&m^s+m_z^h&&&\\&&&m^s+m^h_x&&\\&&&&m^s+m_y^h&\\&&&&&m^s+m_z^h\end{bmatrix}\]

All terms that are not indicated in Equation (46) are equal to zero.

Due to directional dependency the mass matrix contains off-diagonal terms after transformation to the global system. To obtain a diagonal mass matrix after transformation, the off-diagonal terms are simply added to the diagonal for the lumped formulation.

For floating partly submerged cross-sections (special option) the added mass may vary over the element length. Consequently, the added mass matrix is calculated based on numerical integration of Equation (45). This option is restricted to the consistent formulation.

3.4.5. Damping matrix

The damping matrix consists of contributions from structural damping, \(\mathrm {\boldsymbol{c}^s}\), and hydrodynamic damping, \(\mathrm {\boldsymbol{c}^h}\).

\[\boldsymbol{c}=\boldsymbol{c}^s+\boldsymbol{c}^h\]

Hydrodynamic damping is only considered for partly submerged cross-sections (special option).

For description of structural damping, confer Structural Damping Models.

The hydrodynamic damping matrix in local system is calculated according to

\[\begin{array}{l}c_{uu}=0\\\\c_{vv}=\displaystyle \int_{L_0}\boldsymbol{N}_v^Tc_y^h(x)\boldsymbol{N}_v\mathrm {d}x\\\\c_{vv}=\displaystyle \int_{L_0}\boldsymbol{N}_w^Tc_x^h(x)\boldsymbol{N}_w\mathrm {d}x\end{array}\]

where \(\mathrm {c^h_y}\) and \(\mathrm {c^h_z}\) are hydrodynamic damping per unit length.

As indicated, the damping may vary along the element. Thus, the damping matrix is calculated based on numerical integration.

3.5. Beam Element

The beam element is formulated using the concept of co-rotated ghost reference as outlined in Section 2. The present formulation is, apart from a few modifications, identical to the beam element described by Bergan (1984). A detailed discussion of this element together with examples demonstrating its capabilities may be found in Mollestad (1983) and Engseth (1984).

As indicated in Figure 6, the beam has 3 translational and 3 rotational degrees of freedom at each node. They are defined in relation to the local x, y, z-system in the \(\mathrm {C_{0n}}\) configuration. The x-axis is oriented along the secant between the two end nodes 1 and 2 of the element and goes through the centroid of the cross section. The y-axis is taken as the average direction of the principal y-axes at the two beam ends, and the z-axis is finally taken to be perpendicular to the x-y plane. The \(\mathrm {C_{0n}}\) configuration is then oriented along the x-axis with cross-sectional principal axis in the y- and z-directions. It is important to note that the rotational degrees of freedom in Figure 6 express deformational rotations in relation to the co-rotated straight element.

RIFLEX TheoryManual 42 v0 201
Figure 6. Nodal degrees of freedom for beam element.

The beam theory is based on the following assumptions:

  • A plane section of the beam initially normal to the x-axis, remains plane and normal to the x-axis during deformations.

  • Lateral contraction caused by axial elongation is neglected.

  • The strains are small.

  • Shear deformations due to lateral loading are accounted for by modifying the bending stiffness as described in Bell (1987).

  • St. Venant torsion model is applied.

  • Torsional warping resistance is neglected.

  • Coupling between torsion and bending is accounted for by a second order approximation of torsion and bending curvature, however, the coupling is not reflected in the stiffness matrix.

As shown in Figure 7, the displacement of an arbitrary point P with coordinates x, y and z may be expressed as

\[\begin{array}{l}\displaystyle {u}(x,y,z)={u}_0(x)-y\frac{\mathrm {d}{v}_0}{\mathrm {d}x}-z\frac{\mathrm {d}{v}_0}{\mathrm {d}x}\\\\{v}(x,y,z)={v}_0(x)-z\theta \\\\{w}(x,y,z)={w}_0(x)+y\theta \end{array}\]

where \(\mathrm {{u}_0}\), \(\mathrm {{v}_0}\) and \(\mathrm {{w}_0}\) are the displacements of the corresponding point on the reference axis. Equation (8) may now be used to find the axial strain. If all quadratic strain terms that are zero on the x-axis are neglected, and if the quadratic axial strain term is assumed to be negligible, then the Green strain may be expressed as

\[{E}_{xx}={u}_{0,xx}-y{v}_{0,xx}-z{w}_{0,xx}+\frac{1}{2}({v}^2_{0,x}+{w}^2_{0,x})\]

The torsional behaviour of the beam is based on the relationship

\[{M}_\theta =GI_t\theta _x\]

where \(\mathrm {\boldsymbol{M}_\theta }\) is the twist moment and \(\mathrm {GI_t}\) is the torsional stiffness.

RIFLEX TheoryManual 42 v0 212
Figure 7. Prismatic beam.

A standard element formulation gives

\[\begin{array}{l}{u}_0=\boldsymbol{N}_u\boldsymbol{v}_u\quad \mathrm {where}\quad \boldsymbol{v}^T_u=\begin{bmatrix}v_{x1}&v_{x2}\end{bmatrix}\\\\{v}_0=\boldsymbol{N}_v\boldsymbol{v}_v\quad \mathrm {where}\quad \boldsymbol{v}^T_v=\begin{bmatrix}v_{y1}&\theta _{z1}&v_{y2}&\theta _{z2}\end{bmatrix}\\\\{w}_0=\boldsymbol{N}_w\boldsymbol{v}_w\quad \mathrm {where}\quad \boldsymbol{v}^T_w=\begin{bmatrix}v_{z1}&\theta _{y1}&v_{z2}&\theta _{y2}\end{bmatrix}\\\\{\theta }_0=\boldsymbol{N}_\theta \boldsymbol{v}_\theta \quad \mathrm {where}\quad \boldsymbol{v}^T_\theta =\begin{bmatrix}\theta _{x1}&\theta _{x2}\end{bmatrix}\end{array}\]

where\(\mathrm {\boldsymbol{N}_u}\) and \(\mathrm {\boldsymbol{N}_\theta }\) are linear interpolation functions, while \(\mathrm {\boldsymbol{N}_v}\) and \(\mathrm {\boldsymbol{N}_w}\) express cubic interpolation functions.

The interpolation functions are given by

\[ \begin{array}{l} \boldsymbol{N}_u=\boldsymbol{N}_\theta = \begin{bmatrix} 1-\displaystyle \frac{x}{L}&\quad \displaystyle \frac{x}{L} \end{bmatrix} \\\\ \boldsymbol{N}^T_v= \begin{bmatrix} 1-3\displaystyle \frac{x^2}{L^2}+2\frac{x^3}{L^3} \\\\ x-2\displaystyle \frac{x^2}{L}+\frac{x^3}{L^2} \\\\ 3\displaystyle\frac{x^2}{L^2} - 2 \frac{x^3}{L^3} \\ \\ -\displaystyle\frac{x^2}{L} + \frac{x^3}{L^2} \end{bmatrix} \\\\\\ \boldsymbol{N}^T_w= \begin{bmatrix} 1-3\displaystyle \frac{x^2}{L^2}+2\frac{x^3}{L^3} \\\\ -x+2\displaystyle \frac{x^2}{L}-\frac{x^3}{L^2} \\\\ 3\displaystyle\frac{x^2}{L^2} - 2 \frac{x^3}{L^3} \\ \\ \displaystyle\frac{x^2}{L} - \frac{x^3}{L^2} \end{bmatrix} \end{array}\]

3.5.1. Internal loads

Based on the internal virtual work equation, Equation (12), using the results expressed by Equation (50) and Equation (52), it may be shown that the expressions for the internal reaction forces are given by

\[\begin{array}{l}\boldsymbol{S}_u=\displaystyle \int_{L_0}\boldsymbol{N}^T_{u,x}N_{xx}\mathrm {d}x\\\\\boldsymbol{S}_v=\displaystyle \int_{L_0}\boldsymbol{N}^T_{v,xx}M_{y}\mathrm {d}x\\\\\boldsymbol{S}_w=\displaystyle \int_{L_0}\boldsymbol{N}^T_{w,xx}M_{z}\mathrm {d}x\\\\\boldsymbol{S}_\theta =\displaystyle \int_{L_0}\boldsymbol{N}^T_{\theta ,xx}M_{\theta }\mathrm {d}x\end{array}\]

where \(\mathrm {N_{xx}}\), \(\mathrm {M_y}\), \(\mathrm {M_z}\) and \(\mathrm {M_\theta }\) denote the cross-sectional force resultants. \(\mathrm {\boldsymbol{N}_u}\), \(\mathrm {\boldsymbol{N}_v}\), \(\mathrm {\boldsymbol{N}_w}\) and \(\mathrm {\boldsymbol{N}_\theta }\) are taken from Equation (53) with the \(\mathrm {_x}\) and \(\mathrm {_{xx}}\) subscripts denoting first and second derivatives of \(\mathrm {x}\). In the current formulation these resultants will be determined from predefined curves giving the moment-curvature relationship and the relation between axial force and elongation for the given cross section. This eliminates the numerical integration over the cross-section. The axial force \(\mathrm {N_x}\) is positive when in tension, and the sign conventions for the moments \(\mathrm {M_y}\) and \(\mathrm {M_z}\) are given in Figure 8. Note that \(\mathrm {M_y}\) and \(\mathrm {M_z}\) denote the moments giving transverse displacements in the y- and z-directions, respectively. Some higher order terms have been neglected in Equation (54).

RIFLEX TheoryManual 42 v0 241
Figure 8. Sign conventions for moments.

The increment in internal virtual work, Equation (13), is used to calculate the element stiffness matrix. The expression for the geometric stiffness matrix becomes

\[\boldsymbol{k}_G=\begin{bmatrix}0&0&0&0\\0&\boldsymbol{k}_{Gvv}&0&0\\0&0&\boldsymbol{k}_{Gww}&0\\0&0&0&0\end{bmatrix}\]

where

\[\begin{array}{l}\boldsymbol{k}_{Gvv}=\displaystyle \int_{L_0}N_{xx}\boldsymbol{N}^T_{v,x}\mathrm {d}x\\\\\boldsymbol{k}_{Gww}=\displaystyle \int_{L_0}N_{xx}\boldsymbol{N}^T_{w,x}\mathrm {d}x\end{array}\]

From this result one may see that the only geometric effect accounted for by the current beam element, is the influence of axial force on the coupling between transverse forces and transverse displacements. Equation (55) shows no coupling between axial forces and twist, which signifies that torsional stability problems are excluded.

The material part of the tangent stiffness matrix takes the following form

\[\boldsymbol{k}_M=\begin{bmatrix}\boldsymbol{k}_{Muu}&0&0&0\\0&\boldsymbol{k}_{Mvv}&0&0\\0&0&\boldsymbol{k}_{Mww}&0\\0&0&0&\boldsymbol{k}_{M\theta \theta }\end{bmatrix}\]

where

\[\begin{array}{l}\boldsymbol{k}_{Muu}=\displaystyle \int_{L_0}{N}_{T}\boldsymbol{N}^T_{u,x}\boldsymbol{N}_{u,x}\mathrm {d}x\\\\\boldsymbol{k}_{Mvv}=\displaystyle \int_{L_0}{B}_{yT}\boldsymbol{N}^T_{v,xx}\boldsymbol{N}_{v,xx}\mathrm {d}x\\\\\boldsymbol{k}_{Mww}=\displaystyle \int_{L_0}{B}_{zT}\boldsymbol{N}^T_{w,xx}\boldsymbol{N}_{w,xx}\mathrm {d}x\\\\\boldsymbol{k}_{M\theta \theta }=\displaystyle \int_{L_0}GI_{\theta }\boldsymbol{N}^T_{\theta ,x}\boldsymbol{N}_{\theta ,x}\mathrm {d}x\\\end{array}\]

In Equation (58) \(\mathrm {N_T}\), \(\mathrm {{B}_{yT}}\), \(\mathrm {{B}_{zT}}\) denote the tangential cross-sectional moduli where \(\mathrm {N_T}\) is the axial stiffness and \(\mathrm {{B}_{yT}}\) and \(\mathrm {{B}_{zT}}\) are the bending stiffnesses. It is seen from Equation (57) that there is no coupling between axial terms and bending terms, and no coupling between bending about the y- and z-axes. This is because the x-axis is chosen to go through the centroid and y and z to follow the principal axes.

The tangential stiffness is expressed as

\[\boldsymbol{k}=\boldsymbol{k}_G+\boldsymbol{k}_M\]

Eccentricties for the area and shear centers are handled as described in Section [element_local_trans].

3.5.2. Distributed external loads

The external nodal load vector due to distributed external loads is generally written

\[\begin{array}{l}\boldsymbol{S}_u=\displaystyle \int_L\boldsymbol{N}^T_up_x(x)\mathrm {d}x\\\\\boldsymbol{S}_v=\displaystyle \int_L\boldsymbol{N}^T_vp_y(x)\mathrm {d}x\\\\\boldsymbol{S}_w=\displaystyle \int_L\boldsymbol{N}^T_wp_z(x)\mathrm {d}x\\\\\boldsymbol{S}_\theta =\displaystyle \int_L\boldsymbol{N}^T_\theta p_\theta (x)\mathrm {d}x\end{array}\]

where \(\mathrm {p(x)}\) is distributed external loads per unit length.

In most cases it is assumed that the load components are linearly distributed over the element

\[p(x)=\boldsymbol{Np}\]

where \(\mathrm {\boldsymbol{p}}\) is the load intensity vector at the element nodes and \(\mathrm {\boldsymbol{N}}\) is linear displacement function given by

\[\boldsymbol{N}=\begin{bmatrix}\displaystyle 1-\frac{x}{L}&\displaystyle \quad \frac{x}{L}\end{bmatrix}\]

By use of the relation in Equation (61) the external load vector is calculated according to

\[\begin{array}{l}\boldsymbol{S}_u=\displaystyle \int_L\boldsymbol{N}^T_u\boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_x\\\\\boldsymbol{S}_v=\displaystyle \int_L\boldsymbol{N}^T_v\boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_y\\\\\boldsymbol{S}_w=\displaystyle \int_L\boldsymbol{N}^T_w\boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_z\\\\\boldsymbol{S}_\theta =\displaystyle \int_L\boldsymbol{N}^T_\theta \boldsymbol{N}\mathrm {d}x\cdot \boldsymbol{p}_\theta \end{array}\]

where \(\mathrm {\boldsymbol{p}_x}\), \(\mathrm {\boldsymbol{p}_y}\), \(\mathrm {\boldsymbol{p}_z}\) and \(\mathrm {\boldsymbol{p}_\theta }\) are load intensities at the element nodes.

It should be noted that contributions to fixed end moments in \(\mathrm {\boldsymbol{S}_v}\) and \(\mathrm {\boldsymbol{S}_w}\) are disregarded. This is not a fully consistent formulation, but it avoids superposition of the particular solution in the nonlinear stress calculation for the element. Further, the principle of superposition is really not valid in nonlinear analysis. On the other hand, this introduces an error strongly dependent on the actual discretization and consequently one should avoid using adjacent elements with large deviations in length.

By carrying out the integral in Equation (63) disregarding fixed end moments the external load vector is written

\[\boldsymbol{S}^E=\begin{bmatrix}S_{x1}\\S_{y1}\\S_{z1}\\S_{\theta x1}\\S_{\theta y1}\\S_{\theta z1}\\S_{x2}\\S_{y2}\\S_{z2}\\S_{\theta x2}\\S_{\theta y2}\\S_{\theta z2}\end{bmatrix}=\begin{bmatrix}\frac{L}{6}(2p_{x1}+p_{x2})\\\frac{L}{60}(21p_{y1}+9p_{y2})\\\frac{L}{60}(21p_{z1}+9p_{z2})\\\frac{L}{6}(2p_{\theta 1}+p_{\theta 2})\\0\\0\\\frac{L}{6}(p_{x1}+2p_{x2})\\\frac{L}{60}(9p_{y1}+21p_{y2})\\\frac{L}{60}(9p_{z1}+21p_{z2})\\\frac{L}{6}(p_{\theta 1}+2p_{\theta 2})\\0\\0\end{bmatrix}\]

The external load vector can alternatively be calculated using a lumped load formulation, which gives

\[\boldsymbol{S}^E=\begin{bmatrix}S_{x1}\\S_{y1}\\S_{z1}\\S_{\theta x1}\\S_{\theta y1}\\S_{\theta z1}\\S_{x2}\\S_{y2}\\S_{z2}\\S_{\theta x2}\\S_{\theta y2}\\S_{\theta z2}\end{bmatrix}=\frac{L}{2}\begin{bmatrix}p_{x1}\\p_{y1}\\p_{z1}\\p_{\theta 1}\\0\\0\\p_{x2}\\p_{y2}\\p_{z2}\\p_{\theta 2}\\0\\0\end{bmatrix}\]

For elements penetrating the sea surface, contributions to the external load vector from hydrodynamic loads and buoyancy are calculated considering the submerged element length. In this case, linear distribution is assumed from submerged element end to the intersection point between the element and the sea surface.

For partly submerged cross-sections, special option, the contributions from hydrodynamic loading are calculated based on numerical integration of Equation (60). In this case linear displacement functions are used for all degrees of freedom.

Note that hydrodynamic load contributions proportional to structural acceleration and velocity, i.e. added mass and hydrodynamic damping, are included in the mass and the damping matrices respectively.

Eccentricities for the gravity and buoyancy centers are handled as described in Section [element_local_trans].

3.5.3. Mass matrix

The mass matrix contains contributions from the structural and hydrodynamical (added) mass according to

\[\boldsymbol{m}=\boldsymbol{m}^s+\boldsymbol{m}^h\]

where \(\mathrm {\boldsymbol{m}^s}\) is the structural mass matrix and \(\mathrm {\boldsymbol{m}^h}\) is the added mass matrix.

The structural mass matrix is calculated according to

\[\begin{array}{l}\boldsymbol{m}^s_{uu}=m^s\displaystyle \int_{L_0}\boldsymbol{N}^T_u\boldsymbol{N}_u\mathrm {d}x\\\\\boldsymbol{m}^s_{vv}=m^s\displaystyle \int_{L_0}\boldsymbol{N}^T_v\boldsymbol{N}_v\mathrm {d}x\\\\\boldsymbol{m}^s_{ww}=m^s\displaystyle \int_{L_0}\boldsymbol{N}^T_w\boldsymbol{N}_w\mathrm {d}x\\\\\boldsymbol{m}^s_{\theta \theta }=m^s_\theta \displaystyle \int_{L_0}\boldsymbol{N}^T_\theta \boldsymbol{N}_\theta \mathrm {d}x\\\end{array}\]

where \(\mathrm {m^s}\) is the structural mass per unit length and \(\mathrm {m^s_\theta }\) is the mass moment of inertia for rotation about the element x axis per unit length.

The added mass matrix is calculated according to

\[\begin{array}{l}\boldsymbol{m}^h_{uu}=\displaystyle \int_{L_0}\boldsymbol{N}^T_um^h_x\boldsymbol{N}_u\mathrm {d}x\\\\\boldsymbol{m}^h_{vv}=\displaystyle \int_{L_0}\boldsymbol{N}^T_vm^h_y\boldsymbol{N}_v\mathrm {d}x\\\\\boldsymbol{m}^h_{ww}=\displaystyle \int_{L_0}\boldsymbol{N}^T_wm^h_z\boldsymbol{N}_w\mathrm {d}x\\\\\boldsymbol{m}^h_{\theta \theta }=\displaystyle \int_{L_0}\boldsymbol{N}^T_\theta m^h_\theta \boldsymbol{N}_\theta \mathrm {d}x\end{array}\]

where \(\mathrm {m_x^h}\), \(\mathrm {m_y^h}\) and \(\mathrm {m_z^h}\), are added mass per unit length and \(\mathrm {m^h_\theta }\) is added mass due to rotation around the local x-axis.

For elements penetrating the sea surface only the submerged part gives contribution to the added mass matrix.

The mass matrix can alternatively be calculated using a lumped mass formulation. Due to direction dependency the mass matrix will contain off-diagonal terms after transformation to global system. To obtain a diagonal global mass matrix, off-diagonal terms are then simply added to the diagonal.

For floating partly submerged cross-section (special option), the added mass may vary over the element length. Consequently, the added mass matrix is calculated by numerical integration of Equation (68). This option is restricted to consistent formulation.

Mass center eccentricities are handled as described in Section [element_local_trans].

3.5.4. Damping matrix

The damping matrix consists of contributions from structural damping, \(\mathrm {\boldsymbol{c}^s}\), and hydrodynamic damping, \(\mathrm {\boldsymbol{c}^h}\).

\[\boldsymbol{c}=\boldsymbol{c}^s+\boldsymbol{c}^h\]

For description of structural damping, confer Structural Damping Models.

Hydrodynamic damping is only considered for floating partly submerged cross-sections (special option). The hydrodynamic damping matrix is calculated according to

\[\begin{array}{l}\boldsymbol{c}_{uu}=0\\\\\boldsymbol{c}_{vv}=\displaystyle \int_{L_0}\boldsymbol{N}^T_vc^h_y(x)\boldsymbol{N}_v\mathrm {d}x\\\\\boldsymbol{c}_{ww}=\displaystyle \int_{L_0}\boldsymbol{N}^T_wc^h_z(x)\boldsymbol{N}_w\mathrm {d}x\\\\\boldsymbol{c}_{\theta \theta }=\displaystyle \int_{L_0}\boldsymbol{N}^T_\theta c^h_\theta (x)\boldsymbol{N}_\theta \mathrm {d}x\\\\\end{array}\]

where \(\mathrm {c_y^h}\) and \(\mathrm {c_z^h}\) are hydrodynamic damping per unit length due to displacements in y- and z-directions and \(\mathrm {c_\theta ^h}\) is hydrodynamic damping due to rotation around local x-axis.

As indicated the damping may vary along the element. Thus, the damping matrix is calculated based on numerical integration.

3.5.5. Coupled bending and torsion model

The governing equations for the coupled bending and torsion model is developed by representing the curvature in terms of a second order approximation of the deformational motion. The same second order approximation of the deformational motion as used by Battini (2002) is applied, which when introduced into the longitudinal Green strain definition gives the following curvature components,

\[\begin{eqnarray} \kappa_x=\theta _{x,x}+\frac{1}{2}[\theta _{y,x}\theta _z-\theta _y\theta _{z,x}]\\ \kappa_y=\theta _{y,x}+\frac{1}{2}[\theta _{z,x}\theta _x-\theta _z\theta _{x,x}]\\ \kappa_z=\theta _{z,x}+\frac{1}{2}[\theta _{x,x}\theta _y-\theta _x\theta _{y,x}] \end{eqnarray}\]

which can be expressed on matrix format as,

\[ \boldsymbol{\kappa}=[\boldsymbol{I}_3-\frac{1}{2}\boldsymbol{W}(\boldsymbol{\theta })] \boldsymbol{\theta }_{,x}\qquad\boldsymbol{\theta }=\begin{bmatrix}\theta _{x}\\\theta _{y}\\\theta _{z}\end{bmatrix}\quad \boldsymbol{\kappa}=\begin{bmatrix}\kappa_{x}\\\kappa_{y}\\\kappa_{z}\end{bmatrix}\quad \boldsymbol{W}(\boldsymbol{\theta })=\begin{bmatrix}0&-\theta _z&\theta _y\\\theta _z&0&-\theta _x\\-\theta _y&\theta _x&0\\\end{bmatrix}\]

where the bending curvature components about the co-rotated \(\mathrm {y}\) and \(\mathrm {z}\) axes are denoted \(\mathrm {\kappa_y}\) and \(\mathrm {\kappa_x}\), respectively, and \(\mathrm {\kappa_x}\) is the torsion. The interpolations used for the cross-section rotation vector \(\mathrm {\boldsymbol{\theta }}\) and its longitudinal gradient \(\mathrm {\boldsymbol{\theta }_{,x}}\) are based on Equation (52) and Equation (53),

\[ \boldsymbol{\theta }=\boldsymbol{H}_{\theta }\boldsymbol{v}=\begin{bmatrix}0 &0 &0 &0 &0 &0 &0 &0 &0 &0 & \boldsymbol{N}_{\theta } &\\0 &0 &0 &0 &0 &0 &-\boldsymbol{N}_{w,x}& & & &0 &0\\0 &0 &\boldsymbol{N}_{v,x}&&& &0 &0 &0 &0 &0 &0\end{bmatrix}\begin{bmatrix}v_{x1}\\v_{x2}\\v_{y1}\\\theta _{z1}\\v_{y2}\\\theta _{z2}\\v_{z1}\\\theta _{y1}\\v_{z2}\\\theta _{y2}\\\theta _{x1}\\\theta _{x2}\end{bmatrix}\]
\[\boldsymbol{\theta }_{,x}=\boldsymbol{H}_{\theta ,x}\boldsymbol{v}\]

where the cross-section rotation interpolation matrix is denoted \(\mathrm {\boldsymbol{H}_{\theta }}\) to avoid confusion with \(\mathrm {\boldsymbol{N}_{\theta }}\). The virtual curvature is selected as,

\[\delta \boldsymbol{\kappa}=[\boldsymbol{I}_3-\frac{1}{2}\boldsymbol{W}(\boldsymbol{\theta })]\delta \boldsymbol{\theta }_{,x}\]

where the term \(\mathrm {-\frac{1}{2}\boldsymbol{W}(\delta \boldsymbol{\theta })\boldsymbol{\theta }_{,x}}\) is discarded because it introduces spurious internal loads for the standard co-rotational beam formulation. The internal virtual work is given by the left-hand side of Equation (12), which for the virtual work contributions due to bending and torsion can be expressed as follows for a single beam element,

\[ \delta W_{bt}=\int_{L_0}\delta \boldsymbol{\kappa}_{,x}^T\boldsymbol{C}\boldsymbol{\kappa}_{,x}\mathrm {d}x\qquad\boldsymbol{C}=\begin{bmatrix}GI_{\theta }&0&0\\0&B_{zT}&0\\0&0&B_{yT}\end{bmatrix}\]

in which the co-rotated \(\mathrm {y}\) and \(\mathrm {z}\) axes are assumed to coincide with the principal axes for the second area moments of the cross-section. Inserting the curvature vector \(\mathrm {\boldsymbol{\kappa}}\) and the virtual curvature vector \(\mathrm {\delta \boldsymbol{\kappa}}\), and expanding the matrix products yield,

\[\delta W_{bt}=\int_{L_0}\delta \boldsymbol{\theta }_{,x}^T[\boldsymbol{C}+\frac{1}{2}\boldsymbol{W}(\boldsymbol{\theta })\boldsymbol{C}-\frac{1}{2}\boldsymbol{C}\boldsymbol{W}(\boldsymbol{\theta })-\frac{1}{4}\boldsymbol{W}(\boldsymbol{\theta })\boldsymbol{C}\boldsymbol{W}(\boldsymbol{\theta })]\boldsymbol{\theta }_{,x}\mathrm {d}x\]

where the first term inside the bracket represents the standard bending and torsion contributions, see Equation (54). The fourth term is neglected because it is one order less than the second and third terms. The internal load vector representing coupling between bending and torsion is therefore expressed as follows,

\[\boldsymbol{S}_{vw\theta }=\frac{1}{2}\int_{L_0}\boldsymbol{H}_{\theta ,x}^T[\boldsymbol{W}(\boldsymbol{\theta })\boldsymbol{C}-\boldsymbol{C}\boldsymbol{W}(\boldsymbol{\theta })]\boldsymbol{H}_{\theta ,x}\boldsymbol{v}\mathrm {d}x\]

The stiffness matrix contribution due to coupling between bending and torsion is neglected because it does not significantly affect the convergence rate of the solution procedure.

3.5.6. Deformational displacements

In the following the procedure that is used to determine the deformational displacements of the element referred to the \(\mathrm {C_{0n}}\) configuration is described. Before doing so, a few additional coordinate systems need to be defined.

Figure 9 below shows the beam element in deformed configuration with the co-rotated reference configuration \(\mathrm {C_{0n}}\) close to \(\mathrm {C_{n}}\). In addition to the two coordinate systems with base vectors \(\mathrm {\boldsymbol{I}_i}\) and \(\mathrm {\boldsymbol{\bar{i}}_i}\), two more coordinate systems are used. These are shown in Figure 9 and are represented by the base vectors \(\mathrm {\boldsymbol{i}_i^0}\) and \(\mathrm {\boldsymbol{i}_i}\) . The first one is a coordinate system that follows the deformation of the beam end. There is one such system for each beam end, however, only the coordinate system at end b is shown in Figure 9. The coordinate system with base vectors \(\mathrm {\boldsymbol{\bar{i}}_i}\) gives the orientation of the ghost reference configuration \(\mathrm {C_{0n}}\). Figure 9 also illustrates that the nodes may be eccentrically connected to the beam ends. Note that the transformation between the beam end coordinate system and the corresponding nodal point system \(\mathrm {\boldsymbol{\bar{i}}_i}\) does not change during deformation.

RIFLEX TheoryManual 42 v0 316
Figure 9. Deformed beam element.

The translations are easily found as

\[\begin{array}{l}v_{x1}=v_{y1}=v_{z1}=v_{y2}=v_{z2}=0\\\\v_{x2}=L-L_0\end{array}\]

The orientation of the beam end coordinate system may be expressed by the following transformation

\[\boldsymbol{i}^0_i=\boldsymbol{\bar{T}}^0_{ik}\boldsymbol{I}_k\]

where \(\mathrm {\boldsymbol{\bar{T}}^0_{ik}}\) is different at the two beam ends. The transformation between the element coordinate system (\(\mathrm {C_{0n}}\)-configuration) and the global system may stated as

\[\boldsymbol{i}_i=\boldsymbol{T}_{ij}\boldsymbol{I}_j\]

The deformational rotation at the beam ends now appears through the rotation matrix between the \(\mathrm {\boldsymbol{i}^0_i}\)-system in configuration \(\mathrm {C_n}\) and the \(\mathrm {\boldsymbol{i}_i}\)-system in the \(\mathrm {C_{0n}}\)-configuration. By inverting Equation (81) and inserting the result into

\[\boldsymbol{i}^0_i=\boldsymbol{\bar{T}}^0_{ij}\boldsymbol{T}_{kj}\boldsymbol{i}_k\]

This result is analogous to the transformation expressed by \(\mathrm {\boldsymbol{\tilde{T}}_{kj}}\) in Equation (23). By using the result from Equation (25) and comparing this to Equation (82), one finds the following expressions for the beam end rotations

\[\begin{array}{l}\theta _x=\bar{T}^0_{21}T_{31}+\bar{T}^0_{22}T_{32}+\bar{T}^0_{23}T_{33}\\\\\theta _y=\bar{T}^0_{31}T_{11}+\bar{T}^0_{32}T_{12}+\bar{T}^0_{33}T_{13}\\\\\theta _z=\bar{T}^0_{11}T_{21}+\bar{T}^0_{12}T_{22}+\bar{T}^0_{13}T_{23}\end{array}\]

There is one such relation for each beam end. The element displacement vector referred to the local base vectors \(\mathrm {\boldsymbol{\bar{i}}_i}\) is now determined by the results in Equation (79) and

\[\boldsymbol{v}=\begin{bmatrix}v_{x1}&v_{x2}&v_{y1}&\theta _{z1}&v_{y2}&\theta _{z2}&v_{z1}&\theta _{y1}&v_{z2}&\theta _{y2}&\theta _{x1}&\theta _{x2}\end{bmatrix}^T\]

where it is emphasized that the components of \(\mathrm {\boldsymbol{v}}\) refer to the nodal points.

3.5.7. Element local transformations

The element displacement vector is re-organized as follows prior to system equation assembly

\[\boldsymbol{\bar{v}}=\begin{bmatrix}v_{x1}&v_{y1}&v_{z1}&\theta _{x1}&\theta _{y1}&\theta _{z1}&v_{x2}&v_{y2}&v_{z2}&\theta _{x2}&\theta _{y2}&\theta _{z2}&\end{bmatrix}^T\]

which formally is expressed by

\[\boldsymbol{{A}}=\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\end{bmatrix}\]

applied to the element load vectors

\[\boldsymbol{\check{S}}=\boldsymbol{{A}}^T\!\boldsymbol{S}\qquad\boldsymbol{S}=\begin{bmatrix}\boldsymbol{S}_u&\boldsymbol{S}_v&\boldsymbol{S}_w&\boldsymbol{S}_{\theta }\end{bmatrix}^T\]

and the element matrices

\[\begin{array}{l}\boldsymbol{\check{k}}=\boldsymbol{{A}}^T\!\boldsymbol{k}\,\boldsymbol{{A}}\\\\\boldsymbol{\check{c}}=\boldsymbol{{A}}^T\!\boldsymbol{c}\,\boldsymbol{{A}}\\\\\boldsymbol{\check{m}}=\boldsymbol{{A}}^T\!\boldsymbol{m}\,\boldsymbol{{A}}\\\\\end{array}\]
CRS7 allecc
Figure 10. General cross-section. Z,Y is the local/beam element system element system. W,V is the principal axis in the cross section area

The constitutive response of general cross-sections is expressed in the principal axes coordinates to allow for direct use of the shear stiffness model in Bell (1987). The principal axes V and W for the second area moment are formally determined by the requirement \(\int_AV\,W\,\,\mathrm {d}A=0\) where \(\mathrm {A}\) is the cross-section area. The angle of the principal \(\mathrm {V}\)-axis measured relative to the beam element \(\mathrm {Y}\)-axis is denoted \(\mathrm {\theta }\), see Figure 10. The transformation between the principal axes and the beam element \(\mathrm {Y}\)- and \(\mathrm {Z}\)-axis is defined by

\[ \boldsymbol{\hat{T}}_{\theta }=\begin{bmatrix}\boldsymbol{T}_{\theta }&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{T}_{\theta }&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{T}_{\theta }&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{T}_{\theta }\\\end{bmatrix}\qquad\boldsymbol{T}_{\theta }=\begin{bmatrix}1& 0 & 0\\0& \cos(\theta ) & \sin(\theta )\\0& -\sin(\theta ) & \cos(\theta )\\\end{bmatrix}\]

where \(\mathrm {\boldsymbol{I}_{3\times 3}}\) and \(\mathrm {\boldsymbol{0}_{3\times 3}}\) denote the identity and zero matrix of dimension 3 times 3. The transformation of the element internal load vector reads as follows

\[ \boldsymbol{\hat{S}}=\boldsymbol{\hat{T}}_{\theta }^T\boldsymbol{\check{S}}\]

while the material stiffness matrix is transformed according to

\[ \boldsymbol{\hat{k}}=\boldsymbol{\hat{T}}_{\theta }^T\boldsymbol{\check{k}}\,\boldsymbol{\hat{T}}_{\theta }\]

The shear forces and torsion moments in \(\mathrm {\boldsymbol{\hat{S}}}\) refer to the shear center, while the axial force and the bending moments in \(\mathrm {\boldsymbol{\hat{S}}}\) refer to the area center. For general cross-sections the shear and area centers may be eccentrically located. This eccentricity is accounted for by the transformation matrix

\[\boldsymbol{T}_{sa}=\begin{bmatrix}\boldsymbol{I}_{3\times 3}&\boldsymbol{e}_{sa}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{I}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{I}_{3\times 3}&\boldsymbol{e}_{sa}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{I}_{3\times 3}\\\end{bmatrix}\qquad\boldsymbol{e}_{sa}=\begin{bmatrix}0 & Z_a & -Y_a\\-Z_s& 0 & 0\\Y_s& 0 & 0\\\end{bmatrix}\]

where \(\mathrm {(Y_a,Z_a)}\) and \(\mathrm {(Y_s,Z_s)}\) are respectively the coordinates of the area and shear centers expressed in the beam element system, see Figure 10. The eccentricity transformation applies to the element internal load vector and material stiffness matrix according to

\[\boldsymbol{\bar{S}}=\boldsymbol{T}_{sa}^T\boldsymbol{\hat{S}}\]
\[\boldsymbol{\bar{k}}=\boldsymbol{T}_{sa}^T\boldsymbol{\hat{k}}\,\boldsymbol{T}_{sa}\]

It is emphasized that \(\mathrm {\boldsymbol{\bar{S}}}\) and \(\mathrm {\boldsymbol{\bar{k}}}\) are energy-conjugate to the element displacement vector \(\mathrm {\boldsymbol{\bar{v}}}\) which refers to the nodal points of the beam element.

Similarly, general cross-sections may have an eccentric mass center that must be accounted for by the following transformation of the element mass matrix

\[\boldsymbol{T}_{m}=\begin{bmatrix}\boldsymbol{I}_{3\times 3}&\boldsymbol{e}_{m}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{I}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{I}_{3\times 3}&\boldsymbol{e}_{m}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{I}_{3\times 3}\\\end{bmatrix}\qquad\boldsymbol{e}_{m}=\begin{bmatrix}0 & Z_m & -Y_m\\-Z_m& 0 & 0\\Y_m& 0 & 0\\\end{bmatrix}\]
\[\boldsymbol{\bar{m}}=\boldsymbol{T}_{m}^T\boldsymbol{\check{m}}\,\boldsymbol{T}_{m}\]

where \(\mathrm {(Y_m,Z_m)}\) is the mass center expressed in the beam element coordinates. The eccentricity transformation applies also to the load vector containing the gravity and buoyancy force

\[\boldsymbol{\bar{S}}=\boldsymbol{T}_{m}^T\boldsymbol{\check{S}}\]

in which the buoyancy and mass centers are assumed to be coincident.

3.5.8. Equilibrium projection of load vector

The internal load vector in Section [internal_loads] is in equilibrium with respect to the undeformed element. In nonlinear analysis the load vector should be in equilibrium with respect to deformed element, which is particularly important if the deformed length \(\mathrm {L}\) is significantly different from the initial length \(\mathrm {L_0}\). The equilibrium projection of the internal load vector is expressed as follows

\[\boldsymbol{\tilde{S}}=\boldsymbol{{P}}^T\boldsymbol{\bar{S}}\]

where the projection matrix can be shown to be given by

\[\boldsymbol{{P}}^T=\begin{bmatrix}\frac{1}{2} &0 &0 &0 &0 &0 &-\frac{1}{2}&0 &0 &0 &0 &0\\0 &0 &0 &0 &0 &\frac{1}{L}&0&0 &0 &0 &0 &\frac{1}{L}\\0 &0 &0 &0 &-\frac{1}{L}&0 &0&0 &0 &0 &-\frac{1}{L}&0\\0 &0 &0 &\frac{1}{2}&0 &0 &0&0 &0 &-\frac{1}{2}&0 &0\\0 &0 &0 &0 &1 &0 &0&0 &0 &0 &0 &0\\0 &0 &0 &0 &0 &1 &0&0 &0 &0 &0 &0\\-\frac{1}{2}&0 &0 &0 &0 &0 &\frac{1}{2}&0 &0 &0 &0 &0\\0 &0 &0 &0 &0 &-\frac{1}{L}&0&0 &0 &0 &0 &-\frac{1}{L}\\0 &0 &0 &0 &\frac{1}{L}&0 &0&0 &0 &0 &\frac{1}{L}&0\\0 &0 &0 &-\frac{1}{2}&0 &0 &0&0 &0 &\frac{1}{2}&0 &0\\0 &0 &0 &0 &0 &0 &0&0 &0 &0 &1 &0\\0 &0 &0 &0 &0 &0 &0&0 &0 &0 &0 &1\\\end{bmatrix}\]

3.5.9. Transformation of load vector to global system

Prior to assembly of the system equations, the element internal load vector \(\mathrm {\boldsymbol{S}^g}\) is transformed to global coordinates according to

\[\boldsymbol{\hat{T}}=\begin{bmatrix}\boldsymbol{T} &\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{T} &\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{T} &\boldsymbol{0}_{3\times 3}\\\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{0}_{3\times 3}&\boldsymbol{T} \\\end{bmatrix}\qquad\boldsymbol{T}=\begin{bmatrix}\boldsymbol{i}_1^T\\\boldsymbol{i}_2^T\\\boldsymbol{i}_3^T\\\end{bmatrix}\]
\[\boldsymbol{S}^g=\boldsymbol{\hat{T}}^T\boldsymbol{{P}}^T\boldsymbol{\bar{S}}\]

where superscript \(\mathrm {g}\) signifies global coordinates and \(\mathrm {\boldsymbol{i}_i}\) is the element system unit vectors in global coordinates, see Section [deformational_displacements].

3.5.10. Tangent stiffness matrix in the co-rotational framework

The element tangent stiffness matrix presented in this section is at present implemented only for the CRS7 cross-section.

A consistent tangent stiffness matrix requires that increments of the internal load vector are related to increments of the total displacements, which enables quadratic convergence rate of Newton’s method in the vicinity of the solution point. The tangent stiffness relation is derived from

\[\mathrm{d}\boldsymbol{S}^g= \mathrm {d}\boldsymbol{\hat{T}}^T\boldsymbol{P}^T\boldsymbol{\bar{S}}+\boldsymbol{\hat{T}}^T\mathrm {d}\boldsymbol{P}^T\boldsymbol{\bar{S}}+\boldsymbol{\hat{T}}^T\boldsymbol{P}^T\mathrm {d}\boldsymbol{\bar{S}}\]

where the semi-tangential moment contribution to stiffness from the internal moments are neglected, see Krenk (2009). The linearization of the projection matrix \(\mathrm {d}\mathrm {\boldsymbol{P}}\) is also neglected as the change of element length within a step of the Newton solution procedure is assumed negligible. The tangent stiffness relation can therefore be expressed as

\[\mathrm{d}\boldsymbol{S}^g= \mathrm {d}\boldsymbol{\hat{T}}^T\boldsymbol{\tilde{S}}+\boldsymbol{\hat{T}}^T\boldsymbol{P}^T\mathrm {d}\boldsymbol{\bar{S}}\]

where the linearization in the first term provides the geometric stiffness matrix and is associated with rigid body rotation of the element. The second term represents change of internal loads due to deformation and can hence be expressed by the material stiffness matrix \(\mathrm {\boldsymbol{\bar{k}}}\) from Section [element_local_trans] as follows

\[\mathrm{d}\boldsymbol{\bar{S}}=\boldsymbol{\bar{k}}\boldsymbol{P}\mathrm {d}\boldsymbol{\bar{v}}\]

which gives the following material tangent stiffness matrix in global coordinates

\[\boldsymbol{k}_M^g=\boldsymbol{\hat{T}}^T\!\boldsymbol{P}^T\boldsymbol{\bar{k}}\boldsymbol{P}\boldsymbol{\hat{T}}\]

The linearization of the transformation matrix \(\mathrm {\boldsymbol{{T}}}\) is expressed as follows

\[ \mathrm {d}\boldsymbol{T}=\frac{\partial {\boldsymbol{T}}}{\partial {\boldsymbol{\omega }}}\mathrm {d}{\boldsymbol{\omega }}=\begin{bmatrix}\boldsymbol{0}_{1\times 3}\\{\boldsymbol{i}_3}^T\\-{\boldsymbol{i}_2}^T\\\end{bmatrix}\mathrm {d}\omega _x+\begin{bmatrix}-{\boldsymbol{i}_3}^T\\\boldsymbol{0}_{1\times 3}\\{\boldsymbol{i}_1}^T\\\end{bmatrix}\mathrm {d}\omega _y+\begin{bmatrix}{\boldsymbol{i}_2}^T\\-{\boldsymbol{i}_1}^T\\\boldsymbol{0}_{1\times 3}\\\end{bmatrix}\mathrm {d}\omega _z\]

where \(\mathrm {d}\boldsymbol{\omega }\) is the incremental element rigid body rotation. The relation can be expressed in terms of the skew-symmetric matrix \(\mathrm {\boldsymbol{W}}\) and the element transformation matrix \(\mathrm {\boldsymbol{T}}\) as follows,

\[\mathrm{d}\boldsymbol{T}= \begin{bmatrix} 0&\mathrm {d}\omega _z&-\mathrm {d}\omega _y\\ -\mathrm {d}\omega _z&0&\mathrm {d}\omega _x\\ \mathrm {d}\omega _y&-\mathrm {d}\omega _x&0\\ \end{bmatrix} \begin{bmatrix} {{\boldsymbol{i}}_1}^T\\{{\boldsymbol{i}}_2}^T\\ {{\boldsymbol{i}}_3}^T\\ \end{bmatrix} =-\boldsymbol{W}(\mathrm {d}{\boldsymbol{\omega }})\boldsymbol{T}\]

and

\[ \mathrm {d}\boldsymbol{T}^T=\boldsymbol{T}^T\boldsymbol{W}(\mathrm {d}{\boldsymbol{\omega }})\]

The geometric stiffness matrix can now be established

\[\mathrm{d}{\boldsymbol{\hat{T}}}^T\boldsymbol{\tilde{S}}= \begin{bmatrix} \mathrm {d}{\boldsymbol{T}}^T\tilde{\boldsymbol{f}}_1\\ \mathrm {d}{\boldsymbol{T}}^T\tilde{\boldsymbol{m}}_1\\ \mathrm {d}{\boldsymbol{T}}^T\tilde{\boldsymbol{f}}_2\\ \mathrm {d}{\boldsymbol{T}}^T\tilde{\boldsymbol{m}}_2\\ \end{bmatrix}\]

where the components of the equilibrium projected internal load vector are decomposed into contributions representing forces \(\mathrm {\tilde{\boldsymbol{f}}_i}\) and moments \(\mathrm {\tilde{\boldsymbol{m}}_i}\) at element node \(\mathrm {i}\)

\[\boldsymbol{\tilde{S}}= \begin{bmatrix} \tilde{\boldsymbol{f}}_1\\ \tilde{\boldsymbol{m}}_1\\ \tilde{\boldsymbol{f}}_2\\ \tilde{\boldsymbol{m}}_2\\ \end{bmatrix}\]

Introducing the expression for \(\mathrm {d}\boldsymbol{T}^T\) gives

\[\mathrm{d}\boldsymbol{\hat{T}}^T\boldsymbol{\tilde{S}}=\begin{bmatrix}\boldsymbol{T}^T\boldsymbol{W}(\mathrm {d}\boldsymbol{\omega })\tilde{\boldsymbol{f}}_1\\\boldsymbol{T}^T\boldsymbol{W}(\mathrm {d}\boldsymbol{\omega })\tilde{\boldsymbol{m}}_1\\\boldsymbol{T}^T\boldsymbol{W}(\mathrm {d}\boldsymbol{\omega })\tilde{\boldsymbol{f}}_2\\\boldsymbol{T}^T\boldsymbol{W}(\mathrm {d}\boldsymbol{\omega })\tilde{\boldsymbol{m}}_2\\\end{bmatrix}=-\boldsymbol{\hat{T}}^T\begin{bmatrix}\boldsymbol{W}(\tilde{\boldsymbol{f}}_1)\\\boldsymbol{W}(\tilde{\boldsymbol{m}}_1)\\\boldsymbol{W}(\tilde{\boldsymbol{f}}_2)\\\boldsymbol{W}(\tilde{\boldsymbol{m}}_2)\\\end{bmatrix}\mathrm {d}\boldsymbol{\omega}\]

The rigid body rotation of the element can be shown to be given by the following matrix relation for small increments

\[ \mathrm {d}\boldsymbol{\omega }= \boldsymbol{G}\,\boldsymbol{\hat{T}}\,\mathrm {d}\boldsymbol{\bar{v}} \qquad \boldsymbol{G}= \begin{bmatrix} 0&0&0&\frac{1}{2}&0&0&0&0&0&\frac{1}{2}&0&0\\ 0&0&\frac{1}{L}&0&0&0&0&0&-\frac{1}{L}&0&0&0\\ 0&-\frac{1}{L}&0&0&0&0&0&\frac{1}{L}&0&0&0&0\\ \end{bmatrix}\]

and thus

\[\mathrm{d}{\boldsymbol{\hat{T}}}^T\boldsymbol{\tilde{S}}=-\boldsymbol{\hat{T}}^T\!\begin{bmatrix}\boldsymbol{W}(\tilde{\boldsymbol{f}}_1)\\\boldsymbol{W}(\tilde{\boldsymbol{m}}_1)\\\boldsymbol{W}(\tilde{\boldsymbol{f}}_2)\\\boldsymbol{W}(\tilde{\boldsymbol{m}}_2)\\\end{bmatrix}\!\boldsymbol{G}\,\boldsymbol{\hat{T}}\mathrm {d}\boldsymbol{\bar{v}}\]

in which the geometric stiffness matrix in global coordinates reads,

\[\boldsymbol{k}_G^g=-\boldsymbol{\hat{T}}^T\!\begin{bmatrix}\boldsymbol{W}(\tilde{\boldsymbol{f}}_1)\\\boldsymbol{W}(\tilde{\boldsymbol{m}}_1)\\\boldsymbol{W}(\tilde{\boldsymbol{f}}_2)\\\boldsymbol{W}(\tilde{\boldsymbol{m}}_2)\\\end{bmatrix}\!\boldsymbol{G}\,\boldsymbol{\hat{T}}\]

\(\mathrm {\boldsymbol{k}_G^g}\) is a rank 3 matrix which adds stiffness only to the 3 rigid body rotations of the element. It is important to note that the derived geometric stiffness matrix does not add stiffness for the deformation modes of the element. The RIFLEX equation solver handles only symmetric matrices while the derived geometric stiffness matrix is non-symmetric. An artificial symmetrization is therefore applied according to

\[\boldsymbol{k}_G^g=\frac{1}{2}(\boldsymbol{k}_G^g+{\boldsymbol{k}_M^g}^T)\]

3.6. Connector Element

The primary purpose of the connector element is to enable modelling of structural joints. The connector element description here refers to ball joints in RIFLEX. Newer flex joint elements use a master-slave formulation to eliminate fixed degrees of freedom, and allow the user to specify stiffnesses in the free degrees of freedom. The connector element is associated with 2 separate nodes which occupy exactly the same initial position in space. Forces in the connector element are introduced by specified spring stiffness related to translation and rotational degrees of freedom, see Figure 11 below.

RIFLEX TheoryManual 42 v0 333
Figure 11. Connector element.

For application the connector element is exclusively connected to a specified adjacent beam element and used to model internal moment-free joints in the structure only. The procedure used to determine local rotations is similar to the method used for the beam element except for choice of reference system.

As for the beam element, a base vector system for each of the link element ends is introduced, denoted \(\mathrm {\boldsymbol{i}^{0a}}\) and \(\mathrm {\boldsymbol{i}^{0b}}\). In Figure 12 only \(\mathrm {\boldsymbol{i}^{0a}}\) is presented for readability. The transformations between these systems and the corresponding nodal point systems \(\mathrm {\boldsymbol{\bar{i}}^{(a)}}\) and \(\mathrm {\boldsymbol{\bar{i}}^{(b)}}\) (only \(\mathrm {\boldsymbol{\bar{i}}^{(b)}}\) is shown) are given by

\[\begin{array}{l}\boldsymbol{i}^{0a}_i=\boldsymbol{T}^{0a}_{ij}\boldsymbol{\bar{i}}^{(a)}_j=\boldsymbol{T}^{0a}_{ij}\boldsymbol{\bar{T}}^{(a)}_{jk}\boldsymbol{I}_k=\boldsymbol{\bar{T}}^{0a}_{ik}\boldsymbol{I}_k\\\\\boldsymbol{i}^{0b}_i=\boldsymbol{T}^{0b}_{ij}\boldsymbol{\bar{i}}^{(b)}_j=\boldsymbol{T}^{0b}_{ij}\boldsymbol{\bar{T}}^{(b)}_{jk}\boldsymbol{I}_k=\boldsymbol{\bar{T}}^{0b}_{ik}\boldsymbol{I}_k\end{array}\]

where

\[\boldsymbol{T}^{0a}_{ij}=\boldsymbol{T}^{0b}_{ij}=\boldsymbol{T}^{0}_{ij}\]

Due to the definitions of the local coordinate systems used, \(\mathrm {\boldsymbol{T}^0_{ij}}\) also describes the transformation between the adjacent beam end base vector system and the nodal system. Keep in mind that \(\mathrm {\boldsymbol{T}^0_{ij}}\) remains constant through the deformations.

RIFLEX TheoryManual 42 v0 345
Figure 12. Base vector systems for the connector element.

Due to the relation with the adjacent beam element, the ghost reference for this beam element is adopted as reference system for the connector element. Employing Equation (116) and the inverse of Equation (81) gives the following transformation between the connector element end coordinate systems and the ghost reference system of the adjacent beam element

\[\begin{array}{l}\boldsymbol{i}^{0a}_i= \boldsymbol{\bar{T}}^{0a}_{ij}\boldsymbol{T}_{kj}\boldsymbol{i}_k\\\\\boldsymbol{i}^{0b}_i=\boldsymbol{\bar{T}}^{0b}_{ij}\boldsymbol{T}_{kj}\boldsymbol{i}_k\\\end{array}\]

The local rotations of the ends of the connector element related to the ghost reference system are now found by inspecting the differences in the coordinate systems \(\mathrm {\boldsymbol{i}^{0a}_c}\), \(\mathrm {\boldsymbol{i}_i}\) and \(\mathrm {\boldsymbol{i}^{0b}_c}\), \(\mathrm {\boldsymbol{i}_i}\).

The rotations in the connector element should be small, or primarily take place in one plane if they are large.

3.6.1. Internal structural reaction forces and stiffness matrix

The internal structural reaction forces are given by

\[ \boldsymbol{S}^E= \begin{bmatrix}S_{x1}\\ S_{y1}\\S_{z1}\\ S_{\theta x1}\\ S_{\theta y1}\\ S_{\theta z1}\\ S_{x2}\\S_{y2}\\ S_{z2}\\ S_{\theta x2}\\ S_{\theta y2}\\ S_{\theta z2} \end{bmatrix}= \begin{bmatrix} -k_x\cdot (v_{x2}-v_{x1})\\ -k_y\cdot (v_{y2}-v_{y1})\\ -k_z\cdot (v_{z2}-v_{z1})\\ -k_{\theta x}\cdot (\theta _{x2}-\theta _{x1})\\ -k_{\theta y}\cdot (\theta _{y2}-\theta _{y1})\\ -k_{\theta z}\cdot (\theta _{z2}-\theta _{z1})\\ k_x\cdot (v_{x2}-v_{x1})\\ k_y\cdot (v_{y2}-v_{y1})\\ k_z\cdot (v_{z2}-v_{z1})\\ k_{\theta x}\cdot (\theta _{x2}-\theta _{x1})\\ k_{\theta y}\cdot (\theta _{y2}-\theta _{y1})\\ k_{\theta z}\cdot (\theta _{z2}-\theta _{z1})\\ \end{bmatrix}\]

where \(\mathrm {k_x}\), \(\mathrm {k_y}\), \(\mathrm {k_z}\), \(\mathrm {k_{\theta x}}\), \(\mathrm {k_{\theta y}}\) and \(\mathrm {k_{\theta z}}\) are spring stiffnesses associated with the different degrees of freedom.

In the current version of the program the spring stiffnesses are defined by the stiffness of the related beam element as

\[\begin{array}{l}\displaystyle k_x=k_y=k_z=\frac{(EA)}{L_0}\\\\\displaystyle k_{\theta x}=k_{\theta y}=k_{\theta z}=\frac{(EA)}{100L_0}\end{array}\]

where \(\mathrm {(EA)}\) is the stiffness in axial direction and \(\mathrm {L_0}\) is the initial length of the related beam element. However, any of the rotational stiffnesses can be specified equal to zero.

The tangential stiffness matrix is computed using the tangential material moduli of the related element

\[ \boldsymbol{k}=\boldsymbol{k}_M=\\ \begin{bmatrix} k_x&&&&&&-k_x&&&&&\\ &k_y&&&&&&-k_y&&&&\\ &&k_z&&&&&&-k_z&&&\\ &&&k_{\theta x}&&&&&&-k_{\theta x}&&\\ &&&&k_{\theta y}&&&&&&-k_{\theta y}&\\ &&&&&k_{\theta z}&&&&&&-k_{\theta z}\\ -k_x&&&&&&k_x&&&&&\\ &-k_y&&&&&&k_y&&&&\\ &&-k_z&&&&&&k_z&&&\\ &&&-k_{\theta x}&&&&&&k_{\theta x}&&\\ &&&&-k_{\theta y}&&&&&&k_{\theta y}&\\ &&&&&-k_{\theta z}&&&&&&k_{\theta z}\\\ \end{bmatrix}\]

All terms that are not indicated in Equation (121) are equal to zero.

4. Cross Sectional Modelling

4.1. Hyper-elastic constitutive models

The detailed cross-sectional structure of many slender marine structures is normally very complex, e.g. flexible pipes, umbilicals and wires. In a global analysis program, it is therefore convenient to model the global cross-sectional behaviour in terms of force displacement relations. Such relations are axial force versus axial elongation, bending moment versus curvature and torsional moment versus torsional twist:

\[ \begin{array}{l} A=A(\Delta x)\\ M=M(R,\Delta p)\\ M_T=M_T(\Delta \theta ) \end{array}\]

where

  • \(\mathrm {A\:\:\:\::\quad }\) axial force

  • \(\mathrm {\Delta x\::\quad }\) axial elongation per unit length

  • \(\mathrm {M\:\:\::\quad }\) bending moment

  • \(\mathrm {R\:\:\:\::\quad }\) radius of curvature

  • \(\mathrm {\Delta p\:\::\quad }\) pressure difference

  • \(\mathrm {M_T:\quad }\) torsional moment

  • \(\mathrm {\Delta \theta \::\quad }\) axial rotation per unit length

The force/displacement relations can be specified as linear by giving the cross-section stiffness, or as nonlinear by specification of force/displacement relations in tabular form, see Figure 13.

It is also possible to specify a hysteretic bending moment/curvature relationship.

As indicated in Equation (122) it can also be proposed to include a pressure dependent bending stiffness to account for cross-sectional ovalization during bending. This effect is, however, not considered in the present program version.

RIFLEX TheoryManual 42 v0 374
Figure 13. Nonlinear cross section parameters.

4.2. Elasto-plastic bi-linear hysteresis bending model

The bending response of the tensile armour in a flexible pipe may be explained by considering an initially straight pipe section that is gradually bent with constant interlayer radial pressures. For small curvatures, the helix tendons are kept in place by friction, and a linear relationship thus exists between the curvature and the change of tendon stress. In this situation, the bending moment is related to the curvature by a factor equal to the tensile armour bending stiffness contribution in an equivalent cross-section not allowing for interlayer slippage. As the curvature increases, the change of tendon stress exceeds the available friction stress and slippage along the helices starts developing at the neutral axis of bending. When slippage is present around the whole circumference, the tendon stress and the resulting bending moment remain unchanged if the pipe is subjected to a further increase of the curvature. This response can be modelled by a tri-linear moment-curvature relation as proposed in work by Savik (2011). However, the transition between start of slippage to full slippage is small in terms of curvature. This motivates use of a bi-linear bending model which allows for a simpler algorithmic treatment that may achieve better numerical performance than more sophisticated models. Further, the bi-linear bending model can be linearized in a straightforward way leading to a consistent tangent stiffness relation.

The implemented bi-linear model assumes an additive split of the bending moment into elastic and hysteresis parts,

\[\boldsymbol{M}_{tot}=\boldsymbol{M}+\boldsymbol{M}_e\qquad\boldsymbol{M}=\begin{bmatrix}M_{tot,y}&M_{tot,z}\end{bmatrix}^{\boldsymbol{\top}}\]

where the constitutive relation for the elastic part reads,

\[\boldsymbol{M}_e=EJ\boldsymbol{\kappa}\]

in which \(\mathrm {\boldsymbol{\kappa}}\) is the bending curvature vector and the bending stiffness \(\mathrm {EJ}\) is illustrated in the figure Figure 14.

bilinearhystmoment
Figure 14. Bi-linear hysteresis bending moment constitutive parameters

The constitutive model for the hysteresis part is formulated in the framework of computational elasto-plasticity. In order to determine whether slippage is present or not, a slip function \(\mathrm {f}\) is defined according to,

\[f(\boldsymbol{M})=|\boldsymbol{M}|-M_f\leq0\qquad\boldsymbol{M}=\begin{bmatrix}M_y&M_z\end{bmatrix}^{\boldsymbol{\top}}\]

where \(\mathrm {f=0}\) during slip, \(\mathrm {f<0}\) indicates no slippage and \(\mathrm {f>0}\) is an inadmissible state. \(\mathrm {M_f}\) is the friction moment limit, corresponding to the bending moment where slip takes place.

The rate of change of bending moment with respect to time is given by,

\[\dot {\boldsymbol{M}}=K_E(\dot {\boldsymbol{\kappa}}-\dot {\boldsymbol{\kappa}_p})\]

where \(\mathrm {K_E}\) is the elastic modulus corresponding to the tensile armour bending stiffness contribution in an equivalent cross-section not allowing for slippage, which in terms of the stiffness data in the figure Figure 14 is expressed by,

\[K_E=SF\cdot EJ-EJ \qquad SF\geq 1.0\]

In Equation (126), \(\mathrm {\dot {\boldsymbol{\kappa}}}\) is the rate of total bending curvature and \(\mathrm {\dot {\boldsymbol{\kappa}}_p}\) is the rate of bending curvature associated with slip. The latter rate is given by an associative slip rule,

\[\dot {\boldsymbol{\kappa}}_p=\dot {\lambda}\frac{\partial }{\partial \boldsymbol{M}}f(\boldsymbol{M})=\dot {\lambda}\frac{\boldsymbol{M}}{|\boldsymbol{M}|}\qquad\dot {\lambda}\geq0\]

in which \(\mathrm {\dot {\lambda}}\) is the plastic rate parameter, determined by the consistency criterion \(\mathrm {\dot {f}=0}\) if slippage occurs and set to zero if \(\mathrm {f<0}\). In case of no slip, the incremental constitutive relation is obtained by straightforward integration of Equation (126)
with \(\mathrm {\dot {\boldsymbol{\kappa}}_p=\boldsymbol{0}}\). When slippage occurs, the incremental relation valid for finite increments is established by linearizing the moment update scheme. A fully implicit backward Euler scheme is employed to integrate Equation (126) and when linearized results in the following consistent tangent stiffness relation,

\[\begin{bmatrix}dM_y\\dM_z\end{bmatrix}=\frac{K_EM_f}{|\boldsymbol{M}^{tr}|^3}\begin{bmatrix}|\boldsymbol{M}^{tr}|^2-M_y^2&-M_yM_z\\-M_yMz&|\boldsymbol{M}^{tr}|^2-M_z^2\\\end{bmatrix}\begin{bmatrix}d\kappa_y\\d\kappa_z\end{bmatrix}\]

where \(\mathrm {\mathbf{M}^{tr}}\) is the trial moment obtained by assuming no slippage during the finite curvature increment, i.e. \(\mathrm {\dot {\boldsymbol{\kappa}}_p}\) is set to zero in Equation (126) according to,

\[\boldsymbol{M}^{tr}=K_E\,\dot {\boldsymbol{\kappa}}\Delta t+\boldsymbol{M}^p\]

where \(\mathrm {\Delta t}\) is the time increment and \(\mathrm {\boldsymbol{M}^p}\) is the bending moment in the previous equilibrium state.

The moment-curvature response obtained when subjecting the cross-section to a sinusoidal moment loading of 4 kNm about either the y- or z-axis is shown in Figure 15. The cross-section was modelled by \(\mathrm {EJ=20kNm^2}\), \(\mathrm {M_f=1kNm}\) and \(\mathrm {SF=50}\).

The moment response for the same cross-section is shown in the figure Figure 16 when subjected to 2 kNm sinusoidal bending moments about the y- and z-axis which are 90 degrees out of phase. Here, the y-moment load is first ramped to its full value before the z-moment is applied. Observe that the bending moment trajectory is circular, whereas the curvature in the figure Figure 17 undergoes a gradual increase before it stabilizes at a circular trajectory.

hysteresis BL 2Dmomentcurvature
Figure 15. Moment-curvature response for uni-axial bending load
hysteresis BL 3Dmoment
Figure 16. Moment response for bi-axial moment loading
hysteresis BL 3Dcurvature
Figure 17. Curvature response for bi-axial moment loading

4.3. Axial-torsion coupled strain model

The axial-torsion coupling effects are implemented for axisymmetric cross-sections with hyperelastic axial tension/elongation and torsion moment/rotation relations.

The effect is calculated based based on the following relations:

1) Torion strain effect on axial strain is calculated by use of the equation

\[{\tilde{E}_{xx}=E_{xx}+\beta\cdot\frac{-\theta_{x1}+\theta_{x2}}{L}}\]

where \(\mathrm {\beta}\) is the tension/torsion coupling parameter, \(\mathrm {\theta x1}\) and \(\mathrm {\theta x2}\) are torsional rotations at element end 1 and 2 respectively and \(\mathrm {L}\) is the stress free element length.

2) The axial load, \(\mathrm {N_{xx}}\), is calculate based on \(\mathrm {\tilde{E}_{xx}}\) which consequently is used in the internal structural load vector for axial lods.

3) The internal structural torsional reaction forces for element end 1 and 2 are given by

\[ \begin{bmatrix} S_{\theta x1}\\ S_{\theta x2}\\ \end{bmatrix}= \begin{bmatrix} -M_T -\beta \cdot N_{xx}\\ +M_T +\beta \cdot N_{xx}\\ \end{bmatrix}\]

4) The axial-torsion coupling terms in the tangential material stiffness matrix are given by

\[ \begin{bmatrix} K_{x1,\theta x1}\\ K_{x1,\theta x2}\\ K_{x2,\theta x2}\\ K_{\theta x1, x2}\\ \end{bmatrix}= \begin{bmatrix} K_{\theta x1,x1}\\ K_{\theta x2,x1}\\ K_{\theta x2,x2}\\ K_{x2,\theta x1}\\ \end{bmatrix}= \begin{bmatrix} \beta \cdot EA/L\\ -\beta \cdot EA/L\\ \beta \cdot EA/L\\ -\beta \cdot EA/L\\ \end{bmatrix}\]

where \(\mathrm {EA}\) is the axial stiffness.

5. Guideline to System Modelling

Advice with respect to element length and bending stiffness for systems that may have zero tension or compression loads.

5.1. Bar elements:

Element are deformed only in the axial directions. Compressive loads may be limited by two mechanisms:

  1. Lateral motion of nodal points, which imply rigid-body rotation of elements. The shorter the elements, the more rapid the lateral displacement may occur. Too short elements may, however, cause numerical instabilities. Hence, use of short elements to limit axial compressive force in a bar model, will require very short time step.

  2. Use of nonlinear tension-strain relationship, specifying a very low stiffness for negative strain. This will allow the element to shrink, and does not involve element rotations.

5.2. Beam elements:

Compressive load on a beam element may lead to stability problems if the element is long compared with its bending stiffness. In regions where zero tension or compressive loads are expected, element lengths should be limited according to

  1. \(\mathrm {1\leq2(EI/q)^{1/3}}\)

  2. \(\mathrm {1\leq\pi (EI/P_c)^{1/3}}\)

where

  • \(\mathrm {EI\quad }\) = cross-section stiffness, weakest axis if not rotation symmetric

  • \(\mathrm {q\quad \:\:\:\,}\) = largest lateral load intensity, e.g. element weight/length

  • \(\mathrm {P_c\quad \:}\) = largest compressive loads

6. References

Battini, J-M. (2002): Co-rotational beam elements in instability problems PhD Thesis, Department of Mechanics, Royal Institute of Technology, Stockholm.

Bell, K. (1987): Matrisestatistikk Tapir Forlag, Trondheim, Norway.

Bergan, P.G. et al. (1985): FENRIS System Manual, Theory-Program Outline-Data Input Report 83-6064, A.S. Veritec, Høvik, Norway.

Bergan, P.G. et al. (1984): FENRIS Satellite 1 Manual, Theory-Program Outline-Data Input Report 83-6062, A.S. Veritec, Høvik, Norway.

Mollestad, E. (1983): Techniques for Static and Dynamic Solution of Non-linear Finite Element Problems Report No. 83-1, Div. of Structural Mechanics, Norwegian Institute of Technology, Trondheim.

Engseth, A. (1984): Finite Element Collapse Analysis of Tubular Steel Offshore Structures Dr.ing. Thesis, Div. of Marine Structures, The Norwegian Institute of Technology, Trondheim.

Krenk, S. (2009): Non-linear Modeling and Analysis of Solids and Structures Cambridge University Press, Cambridge, UK.

Malvern, L.E. (1969): Introduction to the Mechanics of a Continuous Medium Prentice Hall, Englewood Cliffs, New Jersey.

Remseth, S.N. (1978): Nonlinear Static and Dynamic Analysis of Space Structures Dr.ing. Thesis, Div. of Structural Mechanics, The Norwegian Institute of Technology, Trondheim.

Savik, S. (2011): Theoretical and experimental studies of stresses in flexible pipes Computers & Structures, vol. 89, pp. 2273-91.