Static Finite Element Analysis

1. General

The state of the discretized finite element model is completely determined by the nodal displacement vector. The purpose of the static analysis is to determine the nodal displacement vector so that the complete system is in static equilibrium. The static equilibrium configuration is therefore found as the solution of the following system of equations

\[\boldsymbol{R}^S(\boldsymbol{r})=\boldsymbol{R}^E(\boldsymbol{r})\]

where

  • \(\mathrm {\boldsymbol{r}\quad }\) - Nodal displacement vector including all the degrees of freedom for the system, i.e. displacements for a bar model and displacements and rotations for a beam model. Both displacements and rotations are relative to the stress-free reference configuration. See Finite Element Formulation for details regarding treatment of large rotations in 3D space.

  • \(\mathrm {\boldsymbol{R}^S(\boldsymbol{r})\quad }\) - Internal structural reaction force vector found by assembly of element contributions. Contact forces are also treated as internal reaction forces.

  • \(\mathrm {\boldsymbol{R}^E(\boldsymbol{r})\quad }\) - External force vector assembled from all elements. Accounts for specified external forces, rigid body forces for representation of buoys, clump weights etc. and contribution from distributed loading, i.e. weight, buoyancy and current forces.

Both the internal reaction forces and the external loading will in general be nonlinear functions of the nodal displacement vector. Numerically, static equilibrium is found by application of an incremental loading procedure with equilibrium iterations at each load step, i.e. a so-called incremental-iterative procedure with Euler-Cauchy incrementation.

The basic principle in this approach is to accumulate the external loading in a number of small load increments. The static configuration at each load step is then found by iterative solution of Equation (1) for the accumulated external load vector using the displacement vector from the previous load increment as the initial solution. This strategy is the basis for Section 3.

The bending stiffness will in many cases have a limited effect on the overall static configuration. The catenary solution available for standard systems will therefore provide a close approximation to the final static configuration for a wide range of systems. For such systems, it is possible to apply the catenary solution as the initial configuration for the finite element analysis. This approach gives a significant reduction in the computation time. Further details are given in Section 4.

2. Incremental Equilibrium Iterations

The force imbalance vector at incremental load step \(\mathrm {k}\) is introduced as

\[\boldsymbol{R}_k(\boldsymbol{r})=\boldsymbol{R}^S_k(\boldsymbol{r})-\boldsymbol{R}^E_k(\boldsymbol{r})\]

The static equilibrium configuration at load step \(\mathrm {k}\) corresponds to zero imbalance force, and can be found according to the following iteration procedure:

2.1. Computation of initial values

The initial values for the equilibrium iteration are based on the static equilibrium configuration computed at the previous load step:

\[\Delta \boldsymbol{r}^0_k=-\Big[\frac{\partial \boldsymbol{R}_{k-1}}{\partial \boldsymbol{r}}\Big]^{-1}(\boldsymbol{R}^S_{k-1}-\boldsymbol{R}^E_k)\]
\[\boldsymbol{r}^0_k=\boldsymbol{r}_{k-1}-\Delta \boldsymbol{r}^0_k\]

where \(\mathrm {\displaystyle \frac{\partial \boldsymbol{R}}{\partial \boldsymbol{r}}}\) is the Jacobian matrix (i.e. tangential stiffness matrix), \(\mathrm {\Delta \boldsymbol{r}}\) is the incremental displacement vector and subscripts \(\mathrm {_k}\) and \(\mathrm {_{k-1}}\) denote load step \(\mathrm {k}\) and \(\mathrm {k-1}\), respectively. It should be noted that Equation (3) also will allow for correction of a possible force imbalance at the previous load step.

2.2. Newton-Raphson iteration procedure

Several iteration cycles can be performed to improve the initial solution given by Equation (3) and Equation (4). A Newton-Raphson approach is adopted due to the quadratic convergence rate offered by this procedure. The iteration cycle number is in the following denoted \(\mathrm {j}\). \(\mathrm {\:j=0}\) corresponds to the start values described above while \(\mathrm {j=M}\) is the maximum number of iteration cycles allowed.

The expressions for correction of the displacement vector at interaction cycle \(\mathrm {j}\) are given by

\[\Delta \boldsymbol{r}^j_k=\Big[\frac{\partial \boldsymbol{R}_{k-1}}{\partial \boldsymbol{r}}\Big]^{-1}\boldsymbol{R}^{j-1}_k\]
\[\boldsymbol{r}^j_k=\boldsymbol{r}^{j-1}_k-\Delta \boldsymbol{r}^j_k\]

A reliable convergence criterion is important for any iterative solution procedure, as it is the basis for deciding whether the iterations may be terminated or not. In the current formulation, a modified Euclidean displacement norm is adopted. The norm is based on total displacements

\[\|\boldsymbol{r}_k\|=\frac{1}{N_t}\sum^{N_t}_{i=1}r^2_{ki}\]

where the summation is only carried out over translational components. \(\mathrm {N_t}\) is the number of active translational degrees of freedom.

The iterations are terminated when the following convergence criterion is satisfied

\[\frac{\|\Delta \boldsymbol{r}^j_k\|}{\|\boldsymbol{r}^j_k\|}<\varepsilon _d\]

where \(\mathrm {\varepsilon _d}\) is a specified tolerance requirement. The incremental norm between two iterations is calculated as

\[\|\Delta \boldsymbol{r}^j_k\|=\|\boldsymbol{r}^j_k\|-\|\boldsymbol{r}^{j-1}_k\|\]

A convergence criterion based on energy can also be selected. This convergence criterion is defined by, \(\mathrm {\frac{E_k^j}{E_{ref}}<\varepsilon _e}\)

where the incremental residual energy is given by,

\(\mathrm {E_k^j=\frac{1}{N}\sum_{i=1}^{N}|\Delta r_{ki}^jR_{ki}^{j-1}|}\)

in which \(\mathrm {N}\) is the number of active degrees of freedom including rotations. The denominator in the energy convergence criterion is taken as the largest incremental residual energy for \(\mathrm {j=0}\) ever calculated during the solution process, \(\mathrm {E_{ref}=\text{max}(E_0^0,E_1^0,...,E_k^0,...,E_K^0)}\) where \(\mathrm {K}\) is the final step number within a load group for static analysis or the final time step for dynamic analysis. The reference energy \(\mathrm {E_{ref}}\) is reset at the initial step of a new load group in static analysis.

The energy convergence criteria cannot be used alone, i.e. the displacement criterion in Equation (8) must also be applied. This is because the unbalanced force refers to state \(\mathrm {j-1}\), which may give false convergence in simulations with oscillating convergence behavior.

The iteration procedure is terminated when the convergence tolerances \(\mathrm {\varepsilon _d}\) and \(\mathrm {\varepsilon _e}\) are satisfied or when the specified maximum number of iterations is reached.

The numerical stability of the described procedure is to a large extent dependent on Equation (3) and Equation (4) giving a good starting point for the iteration, which can be obtained by selecting sufficiently small load increments.

It should also be observed that the incremental displacement vector, \(\mathrm {\Delta \boldsymbol{r}^j_k}\), according to Equation (3) and Equation (5), is found by solution of a linear system of equations treating incremental rotations as vector quantities. Incremental rotations are, however, treated by an Euler sequence of rotations when updating the displacement vector (Equation (4) and Equation (6)) according to the procedure described in Updating of the Rotation Matrix. These approximations are valid for small incremental rotations only. Selection of sufficiently small load steps that give small incremental rotations is therefore also crucial for the numerical stability of the equilibrium iteration when rotations are involved, i.e. elements with bending stiffness.

The Jacobian matrix commonly denoted the tangential stiffness matrix, \(\mathrm {\boldsymbol{K}}\), can according to Equation (2) be expressed as

\[\boldsymbol{K}=\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{r}}=\frac{\partial \boldsymbol{R}^S}{\partial \boldsymbol{r}}-\frac{\partial \boldsymbol{R}^E}{\partial \boldsymbol{r}}\]

where

  • \(\mathrm {\displaystyle \frac{\partial \boldsymbol{R}^S_k}{\partial \boldsymbol{r}}=\boldsymbol{K}_S\quad }\) Tangential structural stiffness matrix assembled from all elements

  • \(\mathrm {\boldsymbol{K}_S=\boldsymbol{K}_M+\boldsymbol{K}_G\quad }\) where \(\mathrm {\boldsymbol{K}_M}\) and \(\mathrm {\boldsymbol{K}_G}\) are the material and geometric stiffness matrices, respectively. See Finite Element Formulation for details.

  • \(\mathrm {\displaystyle \frac{-\partial \boldsymbol{R}^E_k}{\partial \boldsymbol{r}}=\boldsymbol{K}_E\quad }\) Load correction stiffness matrix accounting for change in external loading due to change in nodal displacement vector. .

\(\mathrm {\boldsymbol{K}_E}\) is nonzero for displacement-dependent loading only (i.e. current-induced loading and variable submergence rate). The only contribution included in the present formulation is the waterplane stiffness for floating, partly submerged elements, see Hydrodynamic Load on Partly Submerged, Floating Elements for further details.

Contributions to \(\mathrm {\boldsymbol{K}_E}\) from current loading is, however, neglected to avoid non-symmetric stiffness matrices. A slower convergence must therefore be expected when current loading is present. Another practical consequence is that small incremental load steps may be required to obtain a stable numerical solution when applying current loading to current-sensitive systems.

2.3. Modified Newton-Raphson iteration

The Newton-Raphson procedure described above is based on recomputation of \(\mathrm {\boldsymbol{K}}\) at each iteration cycle. However, the Newton-Raphson scheme can alternatively be modified by keeping \(\mathrm {\boldsymbol{K}}\) constant during several iteration cycles. The rate of convergence for this method will generally be somewhat slower than for the pure Newton-Raphson scheme, but can be computationally more efficient since the computation time needed to assemble and triangularize \(\mathrm {\boldsymbol{K}}\) is saved.

3. Finite Element Analysis Starting from Stress-Free Configuration

3.1. Incremental Loading Sequence

In static analysis, it is convenient to distinguish between the following basic load types:

  1. Volume forces, i.e. weight and buoyancy

  2. Specified displacements, i.e. displacements from the stress-free configuration to the final position of nodal points with specified boundary conditions

  3. Specified forces, e.g. nodal point loads

  4. Position-dependent forces, e.g. current forces

The incremental loading procedure starting from the stress-free configuration is organized in a sequence of load conditions where each load condition consists of one or more of the basic load types. Each load condition is applied in a user-defined number of load increments.

Specified displacements, specified forces and current are applied in constant increments while variable increments are used for volume forces, see Section 3.2 for details.

The load types are applied in the order specified by the user. They may be applied in one or more load groups.

Possible seafloor contact is modelled in terms of bottom springs and friction forces, see Contact Formulation for details. Vertical springs representing vertical bottom stiffness are always active while friction forces are activated in the load condition specified by the user, see chapter 6 in the RIFLEX User Manual. Friction forces not activated in static analysis will be activated at the beginning of the dynamic analysis.

An example of the configurations after application of load types 1, 2 and 4 is shown in Figure 1 for a steep wave flexible riser configuration.

For special analyses involving elastic contact between lines, e.g. stinger/pipe interaction, it is also possible to include contact forces as an additional load type, see Contact Formulation for further details.

RIFLEX TheoryManual 42 v0 457
Figure 1. Default incremental loading sequence.

3.2. Procedures to Improve Numerical Stability

The system stiffness for slender structures with small bending stiffness is normally very low during the first incremental load steps. This is because the stiffness against transverse displacements of such systems is mainly governed by the geometric stiffness, which initially is zero when starting from a stress-free configuration. Special considerations are therefore necessary during the first incremental load steps to obtain an efficient and stable numerical procedure.

3.2.1. Use of increased initial geometric stiffness

The stability can be improved by use of an artificially high initial geometric stiffness. This can be achieved by introduction of an artificial initial axial force in all elements which is reduced to zero in linear increments during application of the first load condition. An artificial axial force corresponding to a strain of \(\mathrm {\varepsilon _a=0.01}\) in all elements is applied in the present RIFLEX version.

3.2.2. Use of variable load increments

The overall geometric stiffness will increase as application of the first load condition proceeds. It is therefore convenient to apply the first load condition in gradually increasing incremental load steps allowing for small increments in the first critical phase and larger steps to gain computational efficiency when a realistic geometric stiffness has been established.

The following scaling function is applied:

\[f_p(i)=\frac{1-\mathrm {exp}(\displaystyle \frac{i}{N})^{\displaystyle p}}{1-e}\quad i=1,2,...,N\]

where \(\mathrm {i}\) is the incremental step number, \(\mathrm {N}\) is the total number of load increments in the first load condition and \(\mathrm {p}\) is a shape parameter describing the load increment variation.

The described procedure is in the present version applied for volume forces with \(\mathrm {p=1}\). This approach has been found to give improved efficiency and stability.

4. Finite Element Analysis Starting from Catenary Configuration

4.1. Establishment of Finite Element Model from Catenary Solution

In order to apply the catenary configuration as the initial configuration for the finite element analysis, it is necessary to establish coordinates, nodal rotation matrices, element axial forces and stress-free element lengths from the catenary solution.

Results from the catenary analysis are given in terms of coordinates and tension components at all nodes. The procedures for generation of all relevant finite element data from the catenary solution is summarized in the following.

4.1.1. Coordinates

Coordinates are taken directly from the catenary solution.

4.1.2. Nodal rotation matrices

The nodal rotation matrices defining the rotations relative to the stress-free configuration can be computed directly from the unit vector in the catenary tension direction. The procedure is illustrated in Figure 2 for a catenary configuration in the x-z plane with a horizontal stress-free configuration. Extension to other situations is straightforward.

RIFLEX TheoryManual 42 v0 461
Figure 2. Computation of nodal rotation matrix.

4.1.3. Axial force

The axial element force is found by averaging the tension at the element ends:

\[T_A=\frac{1}{2}\bigg[\sqrt{F^2_{x1}+F^2_{z1}}+\sqrt{F^2_{x2}+F^2_{z2}}\bigg]\]

4.1.4. Stress-free element length

The tensioned element length is found as the chord length, \(\mathrm {l_c}\), of the catenary element, which can be computed directly from the coordinates of the element ends, see Figure 3. The corresponding stress-free element length, \(\mathrm {l^0_c}\), consistent with the axial force, \(\mathrm {T_A}\), is given as

\[l^0_c=\displaystyle \frac{l_c}{(1+\displaystyle \frac{T_A}{EA})}\]

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

It should be noted that this procedure leads to a finite element mesh with stress-free element lengths slightly shorter than specified as input and used in the catenary analysis. The differences are, however, negligible for analysis of most systems when using a reasonable element mesh.

RIFLEX TheoryManual 42 v0 468
Figure 3. Connection between catenary and FEM elements

4.2. Incremental Loading Procedure

In finite element analysis starting from a catenary configuration it is convenient to consider the following basic load types:

  1. Volume forces, i.e. weight and buoyancy

  2. Specified displacements, i.e. from catenary configuration to final position

  3. Specified forces, e.g. nodal point loads

  4. Position dependent forces, e.g. current forces

Application of the remaining basic load types to reach the final static configuration is organized in a sequence of load conditions similar to the approach used for finite element analysis starting from stress-free configuration.

Volume forces are normally applied by the user in the first load condition.

Special considerations and interpretation regarding the basic load types are discussed in the following.

4.2.1. First load condition - volume forces

Volume forces and prescribed displacements from stress-free to final position of terminal points are included in the catenary start solution. Deviations between the catenary and final finite element solution are, however, present due to different mathematical formulations and because bending stiffness is neglected in the catenary solution.

The first load condition will therefore always be a single equilibrium iteration starting from the catenary configuration with volume forces acting. This iteration may fail it there are significant differences between the catenary solution and the final solution due to bending stiffness. It is therefore possible to use a reduced bending stiffness in the first load condition and then apply the bending stiffness in several incremental steps in a later load condition.

4.2.2. Prescribed displacements

The iterative approach on boundary conditions used in the catenary analysis will give deviations between specified translational boundary conditions (i.e. x- and z-coordinates) and boundary conditions computed by the catenary analysis. Further, specified boundary conditions for rotations at the supports will not be satisfied by the catenary analysis due to neglecting bending stiffness. The final position is therefore found by application of the prescribed displacements from the catenary solution to specified positions.

4.2.3. Bending stiffness

As mentioned above, a reduced bending stiffness can be applied in the first equilibrium iteration to improve numerical stability. The bending stiffness is therefore considered as a basic load type in order to enable application of bending stiffness in several incremental steps to reach the correct solution.

Variable increments are crucial for numerical stability. Increments according to Equation (11) with \(\mathrm {p=2.5}\) are used in the present version.

5. Static Parameter Variation Analysis

The intention of the static parameter variation analysis is to study the static system behaviour by considering variation of the supernode position, support vessel position, current profile and specified forces. The variations are relative to the final condition obtained from the static analysis described in Section 3 and Section 4. All parameters can be varied simultaneously in constant increments or alternatively one by one by application of the restart option. The basic analysis principle is similar to the finite element analysis described in Section 3.

For analysis of `troublesome' systems it can also be very useful to consider the static parameter variation combined with the restart option as an extension of the static incremental loading procedure to obtain the final static solution.

6. Guidance to Static Analysis

6.1. General Advice Regarding the Incremental Loading Procedure

The most convenient way to perform static analysis is strongly system-dependant. The following should always be taken into account when considering the strategy for static analysis:

  • system topology

  • seafloor contact

  • definition of stress-free configuration

  • order of load application in the incremental loading sequence

  • static analysis parameters, e.g. number of increments and accuracy

  • sensitivity to the basic static load types

Thus, it is in general very difficult to give quantitative advice for specification of an incremental loading procedure that ensures a stable, efficient numerical solution. Some general quantitative guidelines are, however, given below.

6.1.1. Definition of stress-free configuration for arbitrary systems

  1. Efficiency consideration

Application of specified displacements will normally represent the major computational effort in static analysis starting from stress-free configuration. For the sake of efficiency, the user-defined stress-free configuration for arbitrary (AR) systems should therefore be selected to minimize specified displacements.

  1. Stability considerations when seafloor contact is included

When seafloor contact is involved it is important to define the stress-free configuration so that seafloor contact is avoided during the first critical load steps. This is because seafloor contact in the beginning of the load incrementation will contribute to numerical instability problems.

6.1.2. Order of load application

The order of application of static load conditions should always be selected so that physical instability problems caused by compression or `snap-through' behaviour is avoided. Such problems are normally most likely to occur during application of volume forces and current.

Two typical examples that frequently lead to physical instability are:

  1. Application of volume forces as the first load condition for vertical tensioned risers. This situation will in most cases lead to physical instability due to compression. This problem can easily be solved by simultaneous application of the upper specified tension end and volume forces.

  2. Application of current forces to current sensitive systems can lead to snap-through problems as indicated in Figure 4. This problem can be avoided by considering a modified load application sequence, see Figure 4.

RIFLEX TheoryManual 42 v0 469
Figure 4. Application of current loading.

It is important to observe that the described examples can also occur locally in more complex systems.

6.1.3. Specification of incremental load steps

Each load condition is applied in a user-defined number of incremental load steps, according to Section 3 and Section 4. The required number of load steps to obtain a stable numerical solution depends on the sensitivity to the actual load condition. Load conditions that introduce large relative displacements, rotations and/or large changes in geometric stiffness will require a relatively large number of increments. Typical examples are prescribed displacements and the application of current forces to current sensitive systems (see also comments in Section 3.2). Fewer load steps are normally required for application of volume forces.

From this discussion it can be concluded that a basic intuitive understanding of the numerical procedures described in Section 3 and Section 4 as well as a physical understanding of the system behaviour during application of the static load conditions are crucial for a successful static analysis. Special considerations regarding static analysis involving contact formulations are given in Contact Formulation.

The purpose of this section is to give more specific guidance for analysis of commonly-used flexible riser systems and anchor lines in normal operating conditions. The advice is based on experience gained during static analyses of a wide range of system configurations considering representative cross-sectional properties and moderate current sensitivity.

6.2.1. Recommended incremental loading procedures

Recommended loading sequences and number of load increments are summarized in Table 1. The default loading sequences are recommended for FEM analysis using beam and/or bar elements as well as for the CATFEM analysis using bar elements. It is further recommended that incremental application of bending stiffness is considered for CATFEM analysis when beam elements are used.

It can therefore be concluded that a stable numerical solution is obtained using a straightforward load application for most systems. The recommended number of load increments should, however, be considered as rough guidelines only.

Table 1. Recommended incremental loading procedures for standard analyses
Load Conditions Number of load steps: FEM, beam and/or bar element Number of load steps: CATFEM bar element Number of load steps: CATFEM beam element

1) Volume Forces

5-10

1 (fixed)

1 (fixed)

2) Scaling of bending stiffness

Not relevant

Not relevant

5-10

3) Prescribed displacements

50-200

1

1-5

4) Current

1-10

1-10

1-10

6.2.2. Recommended equilibrium iteration procedure

The true Newton-Raphson iteration scheme should in general be preferred.

It is recommended that the required accuracy measured by the convergence norms described in Section 2 are selected as

  • \(\mathrm {\varepsilon _d=10^{-4}-10^{-6}\quad }\) for a beam element model.

  • \(\mathrm {\varepsilon _d=10^{-7}-10^{-10}\quad }\) for a bar element model.

  • \(\mathrm {\varepsilon _e=10^{-4}-10^{-7}\quad }\) for a beam element model.

The number of equilibrium iterations is normally selected in the range 5-15. If a realistic accuracy criterion is not satisfied in 15 iterations (approximately) this is an indication that the system is unstable, in most cases due to the use of too large load steps.

6.2.3. Selecting method of analysis and system modelling

The CATFEM analysis, as shown in Table 1, is far more efficient than FEM analysis for standard systems. The CATFEM approach should therefore be utilized whenever possible in order to reduce the computational efforts.

The CATFEM approach is in particular efficient when the catenary solution provides a good approximation to the final solution (i.e. systems moderately influenced by bending stiffness). Application of the CATFEM approach will, however, require that the system topology can be described by one of the standard systems, see chapter 2 in the RIFLEX User Manual. Results from FEM and CATFEM analyses will for all practical purposes be identical.

6.3. Error Diagnostics

Failure during incremental application of static load will result in one of the following error messages:

  1. Zero pivot element in system of equations

  2. System of equations failed in stability test

  3. Too large incremental rotations

  4. The algorithm is diverging

The load condition and incremental load step number which gives failure are also specified as a part of the error messages.

Error messages 1 and 2 state that the tangential system stiffness matrix is singular (1) or close to singular (2).

A singular system stiffness matrix is caused by erroneous system modelling, which typically will occur if zero stiffness is present for one or more rigid body motion modes. Such modelling mistakes are in most cases related to boundary condition specifications, connector modelling or inconvenient combinations of beam/bar elements.

Error message 2 signifies that the stiffness matrix is numerically ill conditioned, which also in most cases is due to modelling mistakes, typically by use of beam elements with very high axial stiffness and low bending stiffness.

Error messages 3 and 4 are in most cases (tentatively 95%) caused by using too large incremental load steps, but can also be caused by physical instability problems, see Section 6.1. For identification of physical instability problems it is very useful to study plots of the system configuration during application of the static loads, which is available as output from OUTMOD also when the static analysis fails (see section 8.2 data group B2 in the RIFLEX User Manual).

7. References

Bergan, P.G. Horrigmoe, G., Kråkeland, B. and Søreide, T.H. (1978): "Solution Techniques for Nonlinear Finite Element Problems", International Journal for Numerical Methods in Engineering, Vol. 12.

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

Engseth, A., Bech, A. and Larsen, C.M. (1988): ``Efficient Method for Analysis of Flexible Risers'', BOSS’88, Trondheim, Norway.

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

Striclin, J.A., Haisler, W.E. and von Riesemann, W.A. (1973): ``Evaluation of Solution Procedures for Material and/or Geometrically Nonlinear Structural Analysis'', AIAA Journal, Vol 11, pp. 292-299.