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 crosssectional 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 crosssectional 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
where \(\mathrm {\boldsymbol{x}}\) denotes the position of the particle at the time \(\mathrm {t}\). The displacement vector \(\mathrm {\boldsymbol{u}}\) is defined through
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 corotated 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}\)
In RIFLEX the bar element is formulated using a total Lagrangian description, while the beam element formulation uses a corotated 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
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
The strain tensors may also be expressed as
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 corotated 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
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
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 PiolaKirchhoff 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 PiolaKirchhoff stress tensor referred to\(\mathrm {C_0}\) may be expressed as
Similar to what was done for the strains, this stress tensor may also be referred to a corotated ghost reference \(\mathrm {C_{0n}}\). In this case the stress tensor may be written as
As for the strains, it may be shown that
which means that no transformation of the stress components is necessary when the corotated 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 PiolaKirchhoff stresses, this equation may be written
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.
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 corotated 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)
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
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}}\).
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.
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.
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}\).
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).
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
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
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
thus
\(\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
where
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
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 crosssection forces and small strain theory.
The element is assumed to be straight with an initial crosssectional 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
where
\(\mathrm {\Delta x=x_2x_1}\)
\(\mathrm {\Delta y=y_2y_1}\)
\(\mathrm {\Delta z=z_2z_1}\)
Based on a total Lagrangian formulation and linear displacement functions, the Green strain is expressed
The PiolaKirchhoff stress \(\mathrm {S_f}\) is found through the constitutive law
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
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
where \(\mathrm {L_0}\) is initial, stress free element length and \(\mathrm {EA}\) is the axial stiffness. The strain, \(\mathrm {\varepsilon }\), is given by
3.4.1. Internal loads
Based on the virtual work equation, Equation (12), the internal reaction force vector is expressed by
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
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
where
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
3.4.3. Distributed external loads
The consistent external nodal load vector is generally written
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
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
The external load vector can alternatively be calculated using a lumped load formulation, which gives
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 crosssections (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
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
where \(\mathrm {\boldsymbol{m}^s}\) is structural mass per unit length and \(\mathrm {\boldsymbol{N}}\) is the displacement function matrix given by
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 directiondependent added mass matrix is calculated in the local system according to
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 zaxes, 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
All terms that are not indicated in Equation (46) are equal to zero.
Due to directional dependency the mass matrix contains offdiagonal terms after transformation to the global system. To obtain a diagonal mass matrix after transformation, the offdiagonal terms are simply added to the diagonal for the lumped formulation.
For floating partly submerged crosssections (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}\).
Hydrodynamic damping is only considered for partly submerged crosssections (special option).
For description of structural damping, confer Structural Damping Models.
The hydrodynamic damping matrix in local system is calculated according to
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 corotated 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, zsystem in the \(\mathrm {C_{0n}}\) configuration. The xaxis 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 yaxis is taken as the average direction of the principal yaxes at the two beam ends, and the zaxis is finally taken to be perpendicular to the xy plane. The \(\mathrm {C_{0n}}\) configuration is then oriented along the xaxis with crosssectional principal axis in the y and zdirections. It is important to note that the rotational degrees of freedom in Figure 6 express deformational rotations in relation to the corotated straight element.
The beam theory is based on the following assumptions:

A plane section of the beam initially normal to the xaxis, remains plane and normal to the xaxis 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
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 xaxis are neglected, and if the quadratic axial strain term is assumed to be negligible, then the Green strain may be expressed as
The torsional behaviour of the beam is based on the relationship
where \(\mathrm {\boldsymbol{M}_\theta }\) is the twist moment and \(\mathrm {GI_t}\) is the torsional stiffness.
A standard element formulation gives
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
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
where \(\mathrm {N_{xx}}\), \(\mathrm {M_y}\), \(\mathrm {M_z}\) and \(\mathrm {M_\theta }\) denote the crosssectional 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 momentcurvature relationship and the relation between axial force and elongation for the given cross section. This eliminates the numerical integration over the crosssection. 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 zdirections, respectively. Some higher order terms have been neglected in Equation (54).
The increment in internal virtual work, Equation (13), is used to calculate the element stiffness matrix. The expression for the geometric stiffness matrix becomes
where
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
where
In Equation (58) \(\mathrm {N_T}\), \(\mathrm {{B}_{yT}}\), \(\mathrm {{B}_{zT}}\) denote the tangential crosssectional 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 zaxes. This is because the xaxis is chosen to go through the centroid and y and z to follow the principal axes.
The tangential stiffness is expressed as
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
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
where \(\mathrm {\boldsymbol{p}}\) is the load intensity vector at the element nodes and \(\mathrm {\boldsymbol{N}}\) is linear displacement function given by
By use of the relation in Equation (61) the external load vector is calculated according to
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
The external load vector can alternatively be calculated using a lumped load formulation, which gives
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 crosssections, 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
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
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
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 xaxis.
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 offdiagonal terms after transformation to global system. To obtain a diagonal global mass matrix, offdiagonal terms are then simply added to the diagonal.
For floating partly submerged crosssection (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}\).
For description of structural damping, confer Structural Damping Models.
Hydrodynamic damping is only considered for floating partly submerged crosssections (special option). The hydrodynamic damping matrix is calculated according to
where \(\mathrm {c_y^h}\) and \(\mathrm {c_z^h}\) are hydrodynamic damping per unit length due to displacements in y and zdirections and \(\mathrm {c_\theta ^h}\) is hydrodynamic damping due to rotation around local xaxis.
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,
which can be expressed on matrix format as,
where the bending curvature components about the corotated \(\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 crosssection rotation vector \(\mathrm {\boldsymbol{\theta }}\) and its longitudinal gradient \(\mathrm {\boldsymbol{\theta }_{,x}}\) are based on Equation (52) and Equation (53),
where the crosssection rotation interpolation matrix is denoted \(\mathrm {\boldsymbol{H}_{\theta }}\) to avoid confusion with \(\mathrm {\boldsymbol{N}_{\theta }}\). The virtual curvature is selected as,
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 corotational beam formulation. The internal virtual work is given by the lefthand 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,
in which the corotated \(\mathrm {y}\) and \(\mathrm {z}\) axes are assumed to coincide with the principal axes for the second area moments of the crosssection. Inserting the curvature vector \(\mathrm {\boldsymbol{\kappa}}\) and the virtual curvature vector \(\mathrm {\delta \boldsymbol{\kappa}}\), and expanding the matrix products yield,
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,
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 corotated 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.
The translations are easily found as
The orientation of the beam end coordinate system may be expressed by the following transformation
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
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
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
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
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 reorganized as follows prior to system equation assembly
which formally is expressed by
applied to the element load vectors
and the element matrices
The constitutive response of general crosssections 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 crosssection 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
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
while the material stiffness matrix is transformed according to
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 crosssections the shear and area centers may be eccentrically located. This eccentricity is accounted for by the transformation matrix
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
It is emphasized that \(\mathrm {\boldsymbol{\bar{S}}}\) and \(\mathrm {\boldsymbol{\bar{k}}}\) are energyconjugate to the element displacement vector \(\mathrm {\boldsymbol{\bar{v}}}\) which refers to the nodal points of the beam element.
Similarly, general crosssections may have an eccentric mass center that must be accounted for by the following transformation of the element mass matrix
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
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
where the projection matrix can be shown to be given by
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
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 corotational framework
The element tangent stiffness matrix presented in this section is at present implemented only for the CRS7 crosssection.
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
where the semitangential 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
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
which gives the following material tangent stiffness matrix in global coordinates
The linearization of the transformation matrix \(\mathrm {\boldsymbol{{T}}}\) is expressed as follows
where \(\mathrm {d}\boldsymbol{\omega }\) is the incremental element rigid body rotation. The relation can be expressed in terms of the skewsymmetric matrix \(\mathrm {\boldsymbol{W}}\) and the element transformation matrix \(\mathrm {\boldsymbol{T}}\) as follows,
and
The geometric stiffness matrix can now be established
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}\)
Introducing the expression for \(\mathrm {d}\boldsymbol{T}^T\) gives
The rigid body rotation of the element can be shown to be given by the following matrix relation for small increments
and thus
in which the geometric stiffness matrix in global coordinates reads,
\(\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 nonsymmetric. An artificial symmetrization is therefore applied according to
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 masterslave 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.
For application the connector element is exclusively connected to a specified adjacent beam element and used to model internal momentfree 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
where
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.
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
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
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
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
All terms that are not indicated in Equation (121) are equal to zero.
4. Cross Sectional Modelling
4.1. Hyperelastic constitutive models
The detailed crosssectional 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 crosssectional 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:
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 crosssection 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 crosssectional ovalization during bending. This effect is, however, not considered in the present program version.
4.2. Elastoplastic bilinear 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 crosssection 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 trilinear momentcurvature 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 bilinear bending model which allows for a simpler algorithmic treatment that may achieve better numerical performance than more sophisticated models. Further, the bilinear bending model can be linearized in a straightforward way leading to a consistent tangent stiffness relation.
The implemented bilinear model assumes an additive split of the bending moment into elastic and hysteresis parts,
where the constitutive relation for the elastic part reads,
in which \(\mathrm {\boldsymbol{\kappa}}\) is the bending curvature vector and the bending stiffness \(\mathrm {EJ}\) is illustrated in the figure Figure 14.
The constitutive model for the hysteresis part is formulated in the framework of computational elastoplasticity. In order to determine whether slippage is present or not, a slip function \(\mathrm {f}\) is defined according to,
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,
where \(\mathrm {K_E}\) is the elastic modulus corresponding to the tensile armour bending stiffness contribution in an equivalent crosssection not allowing for slippage, which in terms of the stiffness data in the figure Figure 14 is expressed by,
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,
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,
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,
where \(\mathrm {\Delta t}\) is the time increment and \(\mathrm {\boldsymbol{M}^p}\) is the bending moment in the previous equilibrium state.
The momentcurvature response obtained when subjecting the crosssection to a sinusoidal moment loading of 4 kNm about either the y or zaxis is shown in Figure 15. The crosssection was modelled by \(\mathrm {EJ=20kNm^2}\), \(\mathrm {M_f=1kNm}\) and \(\mathrm {SF=50}\).
The moment response for the same crosssection is shown in the figure Figure 16 when subjected to 2 kNm sinusoidal bending moments about the y and zaxis which are 90 degrees out of phase. Here, the ymoment load is first ramped to its full value before the zmoment 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.
4.3. Axialtorsion coupled strain model
The axialtorsion coupling effects are implemented for axisymmetric crosssections 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
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
4) The axialtorsion coupling terms in the tangential material stiffness matrix are given by
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:

Lateral motion of nodal points, which imply rigidbody 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.

Use of nonlinear tensionstrain 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

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

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

\(\mathrm {EI\quad }\) = crosssection 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, JM. (2002): Corotational 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, TheoryProgram OutlineData Input Report 836064, A.S. Veritec, Høvik, Norway.
Bergan, P.G. et al. (1984): FENRIS Satellite 1 Manual, TheoryProgram OutlineData Input Report 836062, A.S. Veritec, Høvik, Norway.
Mollestad, E. (1983): Techniques for Static and Dynamic Solution of Nonlinear Finite Element Problems Report No. 831, 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): Nonlinear 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. 227391.