Finite Element Formulation

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 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}\)

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.

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 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

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 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

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

As for the strains, it may be shown that

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

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

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 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}}\).

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

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.

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).

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

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.

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

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

The Piola-Kirchhoff 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

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.

The beam theory is based on the following assumptions:

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 x-axis 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

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.

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

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.

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:

where

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.

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,

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 elasto-plasticity. 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 cross-section 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 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.

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

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 axial-torsion coupling terms in the tangential material stiffness matrix are given by

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

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

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

where