Proceedings of OMAE’06

25st International Conference on Offshore Mechanics and Arctic Engineering

June 4-9, 2006 – Hamburg, Germany

OMAE2006-92153

IMPLICIT AND EXPLICIT IMPLEMENTATION OF THE DYNAMIC RELAXATION

METHOD FOR THE DEFINITION OF INITIAL EQUILIBRIUM CONFIGURATIONS
OF FLEXIBLE LINES
Danilo Machado Lawinscky da Silva Breno Pinheiro Jacob Marcos Vinícius Rodrigues
LAMCSO – Laboratory of Computer Methods and Offshore Systems
PEC/COPPE – Post-Graduate Institute of the Federal University of Rio de Janeiro, Civil Eng. Department
Rio de Janeiro, RJ, Brazil

ABSTRACT The research in analytical formulations based on the

Recent activities in the offshore oil exploitation industry catenary equations is very extensive [2] and focus in topics as
require new structural concepts employing flexible lines (both the treatment of lines composed by materials of different
mooring lines and risers). Such systems present increasingly properties; including buoys and distributed floaters;
complex configurations, with dynamic nonlinear behaviour; considering irregular sea bottom; considering the flexural
therefore, the use of efficient numerical solution procedures, stiffness, and arbitrary loads such as current. However, there
based on the Finite Element Method, becomes mandatory for are limitations for the development and efficient use of
their analysis. analytical formulations, when one considers complex systems
The usual analysis procedure for flexible lines by the FEM not only such as those mentioned above, but also systems with
is based in the calculation of an initial, stable static equilibrium several lines connected amongst themselves and to tendons, or
configuration in order to define the finite element mesh. for instance offloading hoses, where the flexural stiffness is
Usually this configuration is obtained by the classic catenary relatively more important and the hydrostatic behavior is
equations. However, in more complex problems these nonlinear and should be treated with a more rigorous
equations cannot be applied. Therefore, the objective of this formulation.
work is to present the use of a more general finite element Therefore, it is seen that procedures based on the catenary
approximation, associated to dynamic relaxation algorithms. equations may not be feasible for more complex systems. The
Such algorithms can be started from arbitrary configurations, determination of static stable configurations for such systems
not necessarily in equilibrium. The resulting procedure is should then use procedures based on more general Finite
accurate, robust, and avoids numerical problems such as the ill- Element approximations. These procedures require robust
conditioning of the tangent stiffness matrix, allowing the static solution schemes that can be started from a FE mesh
equilibrium configuration to be obtained in an efficient way. corresponding to an easily defined configuration, not
necessarily in balance under the action of the loads. Therefore,
INTRODUCTION
in that procedure an initial geometry (free of loads) should be
Systems of flexible lines and cables are frequently used in specified arbitrarily.
civil and oceanic engineering design. Such lines are typically
used in mooring systems to position offshore platforms; in riser However the conventional solution procedures for nonlinear
systems and offloading floating hoses to transport fluids static problems discretized by FE, which is based on the
resulting of the petroleum exploitation process. iterative Newton-Raphson method, can become inefficient for
several reasons. For example, the use of an initial geometry
Flexible lines present strong dynamic nonlinear behavior, free of loads makes the tangent stiffness matrix becomes
with environmental loads dependent on time and nonlinear singular. Supplying an arbitrary state of initial tensions can
boundary conditions, therefore requesting complex numeric eliminate that problem, but this leads again to the initial
solution procedures, based on the Finite Element Method. In problem, how to define a configuration under such state of
the usual procedure, the first stage is the establishment of an tensions. Also, the very low or absent flexural stiffness, as well
initial static equilibrium configuration under dead weight only. as the compression stiffness, leads to a highly ill-conditioned
Traditionally this initial static configuration of flexible lines is tangent stiffness matrix. In static analyses this problem can
obtained by the classic catenary equations [1]. This initial cause the lack of convergence of the iterative Newton-Raphson
configuration is then employed to define the finite element process.
mesh that is subsequently employed in the full static and
dynamic analyses. The objective of this work then consists of establishing
robust and efficient static solution procedures, oriented to the

1 Copyright © 2006 by ASME

definition of static stable configurations for complex systems of support their own weight. When the weight loading is small
flexible lines, under the action of all static load components. compared to the design loads, the dead loaded state is not very
The static configurations obtained can be used in a subsequent different from the “weightless” state. Such structures have an
dynamic analysis, guaranteeing a solution that is stable, intrinsic stiffness.
efficient and free of transient. To overcome the problems Cable Structures (suspended roofs, moorings, large ocean
previously mentioned, associated to the use of analytical structures) are examples of structures which are likely to
catenary formulations or the Newton-Raphson method as present difficult in obtaining an initial state suitable for
solution algorithm for problems discretizated by FEM, dynamic structural calculations. The common denominator in all of these
relaxation algorithms will be considered. is the fact that in an unloaded, weightless condition they are
The Dynamic Relaxation Method (DRM) has been used to essentially mechanisms with negligible intrinsic stiffness. This
obtain the static solution of structural mechanics problems in general type of problem has strong geometric nonlinearities
general. This method is based on the fact that the static solution even if the initial singularity is overcome.
is the steady-state part of the transient response. In this case,
BASIC FORMULATION
the transient part of the solution is not of interest, only the
steady-state response is desired. This technique is classified in The general discrete form of the static equilibrium
the literature as a pseudo-transient or pseudo-dynamic method equations can be written as
[3]. FExt(t) - FInt(t) = 0 (1)
Considering that the static response is the limit case of the where FExt(t) represents the vector of external loads for each
dynamic response damped by a long period of time, DRM d.o.f. of the Finite Element mesh; FInt(t) represents the internal
makes use of "artificial" dynamic effects of inertia and forces. In geometrically nonlinear problems the internal forces
damping to build a conditioning mechanism for the tangent depend on the displacements. When the external loads also
stiffness matrix, in the case of an implicit formulation, or depend on the displacements (as in pressure loading), the
simply in terms of forces in the case of an explicit formulation. loading is nonconservative. Combinations of the two effects are
To accelerate the convergence, artificial damping can be used, common.
leading to the gradual relaxation of the inertia effects. A direct Although equation (1) represents the equilibrium equations,
advantage is obtained in the transformation of a static problem it does not describe a method for finding the values of the
in a dynamic problem; the dynamic terms (inertia and damping) displacements which satisfy equilibrium, nor does it give
act as a conditioning mechanism and the problem of the ill- explicit information on how to calculate the internal and
conditioning of the tangent stiffness matrix disappears. external forces.
The resulting procedure can be started from an arbitrary Newton-Type Methods
configuration, not necessarily in equilibrium under the action The most popular method for solving equation (1) in
of the loads. In short, it can overcome the pointed limitations structural applications is the modified Newton-Raphson
for the catenary equations and for the FEM associated to method. Its general form is
Newton-Raphson Methods, allowing the definition of initial (k -1) (k) (k) (k-1)
stable configuration of flexible lines with low computational K T DU =R - F n+1 (2)
cost. with
(k) (k -1) (k)
CHARACTERISTICS OF THE PROBLEM Un+1 = U n+1 + DU (3)
Structural systems which respond with large deflections are where KT is the tangent stiffness matrix; k is the iteration index;
said to exhibit a geometric nonlinearity. The deformations must (k)
DU are the displacement increments between two successive
be included in the expressions of static equilibrium for such (k) (k -1)
structures. Typically, the dominant terms in the deformations iterations; R and F are respectively the external and
(k -1)
are rotations; however, problems involving large stretching are internal forces; and U n+1 are the displacements measured from
also encountered. The nonlinear equations can be solved using the initial reference state as estimated in the kth iteration.
incremental/iterative methods based on the modified Newton- It is seen that the Newton-Raphson iteration scheme
Raphson method. One of the most popular incremental solution (k)
forms is the so-called incremental self-correcting method. It is involves the approximation of the nonlinear term Fn+1 by a
(k -1) (k)
essentially a series of simple linearized steps with a single truncated Taylor series comprised by F n+1 and KTDU .
Newton-Raphson iteration after each one. Equation (2) is then applied recursively until both the residual
All of these methods rely on an estimate of the initial and the displacement increments are smaller than a preselected
stiffness of the structure. This implies that a stable initial state level. The usual form of the Newton-Raphson procedure is
should be defined. In the most common nonlinear structural obtained by recalculation of the stiffness matrix at the
problems, the initial reference state is quite easy to define. This beginning of each iteration. Modified schemes use
is so because most structures have sufficient stiffness to

2 Copyright © 2006 by ASME

approximations for the stiffness matrix and/or hold the matrix addition to the externally applied discretized forces and the
constant for a number of iterations. internal nodal forces of the continuum or structure, and may be
Two major problems detract from the use of these Newton written as
methods for the class of problems under discussion. First is the FI(t) + FD(t) + FInt(t) = FExt(t) (4)
conditionally stable nature of these methods. There is an **
where FI(t) are the inertia forces, FI(t) = MU; FD(t) are the
interval of convergence around the correct solution, within *

which a particular method will converge. When the starting damping forces, FD(t) = CU; FInt(t) are the internal forces,
estimate is outside this interval, convergence cannot be FInt(t) = F(U); and FExt(t) are the external forces, FExt(t) = R, all
obtained without introducing special procedures. These of them being time-dependent. Denoting by n the current time
convergence enhancement schemes are usually limited in scope increment, the semi-discrete finite element equations of motion
and highly specialized so that their effectiveness is quite for time increment n + 1 may be written as
** *
restricted. The interval of convergence is highly problem- MUn+1 + CUn+1 + F(Un+1) = Rn+1 (5)
dependent and sensitive to variations in the solution method. In where M is a mass matrix, C is the damping matrix, F is the
most situations this interval is not explicitly calculated. The internal force vector and R is a vector of external loads. The
second problem is the need for making an initial estimate of the ** *
parameters U, U and U represent the acceleration, velocity and
solution. This boils down to requiring the very thing that is
displacement vectors, respectively. It should be noted that
being sought. A common means of starting a Newton-Raphson
(0) internal force vector F is a function of the displacement as
solution is to assume; Un+1 is zero and use the initial geometry
F(Un+1) = K0Un+1 + Q0(Un+1) (6)
to calculate a linear step. Quite often the stiffness matrix in the
initial configuration is singular or very ill-conditioned. This where K0 represents the linear stiffness matrix and Q is a
means that the first step cannot be solved or the solution will be vector of non-linear terms which may arise from
very far from the starting configuration. It is also quite likely to geometric/material non-linearities, contact/boundary conditions
be far from the correct solution. This ill-conditioning can be and follower forces. The solution of equation (5) may be
compensated for to some extent by introducing an artificial obtained by employing any one of the many direct time
stiffness. Even with this artificial stiffness there is no guarantee integration algorithms.
that a solution starting from a very poor quality initial Conventional Explicit Solution Method
configuration will lead to a first step result that is anywhere The use of an explicit, rather than implicit, time integration
near the correct solution. Without that guarantee, there is a technique to solve the semi-discrete system of equations given
strong likelihood that additional iterations will diverge. by equation (5) is attractive for parallel computations. Central-
The Idea of Dynamic Relaxation difference approximations are typically employed for the
Dynamic relaxation (DR) is a technique by which the static derivatives. That is,
solution is obtained by determining the steady-state response to ** 1
Un = 2 (Un-1 - 2Un + Un+1) (7)
the transient dynamic analysis for an autonomous system. In h
this case, the transient part of the solution is not of interest, * 1
only the steady-state response is desired. Since the transient Un = 2h (Un+1 - Un-1) (8)
solution is not desired, fictitious mass and damping matrices Substituting these approximations into equation (5), the
which no longer represent the physical system are chosen to following equations for displacement at time n + 1 may be
accelerate the determination of the steady-state response. These obtained
matrices are redefined (using existing equations) so as to
produce the most rapid convergence. A clear advantage is ⎛ 12 M + 1 C⎞Un+1 = Rn - F(Un) + 22 MUn
⎝h 2h ⎠ h
gained through this static-to-dynamics transformation; the (9)
1 1
dynamic term(s) (inertia and damping) acts as the conditioning - ⎛ 2 M - 2h C⎞Un-1
mechanism and the problem of an ill-conditioned tangent ⎝h ⎠
stiffness matrix disappear. where h is the time step. These equations are linear (even
The idea of dynamic relaxation has been exploited widely for non-linear problems), and the left-hand side matrix is
in a variety of structural analysis. A history of DR and the constant unless the time step h changes during the solution
development of an adaptive algorithm have been presented by process. Also, if diagonal mass and damping matrix are used,
Underwood [4]. these equations represent an uncoupled system of algebraic
equations in which each solution component may be computed
GOVERNING STRUCTURAL DYNAMICS EQUATIONS independently. For transient dynamic analysis, a time history of
In dynamic analysis the governing semi-discrete equations displacement (system response) is sought. Mass and damping
of motion are obtained by considering the static equilibrium at vector that best model the physical properties of the system are
time t which includes the effect of acceleration-dependent used. Techniques for estimating the maximum allowable time
inertia forces and velocity-dependent damping forces in step are available, such that the time step size may change

3 Copyright © 2006 by ASME

during the transient dynamic analysis. As such, explicit time the initial conditions. The static loads are then applied and held
integration techniques are attractive candidates for constant as the system is allowed to move dynamically until
implementation on parallel computers. These techniques motion dies out. This effectively avoids the singularity problem
generally have low memory and communication requirements if damping and mass terms are appropriately defined.
but are also only conditionally stable numerically. If the physical characteristics are used to model the mass
Conventional Implicit Solution Method and damping, one can expect very large transients which persist
In general, the dynamic response of inertial problems is for long times. This means excessive solution costs and
more efficiently attained with implicit time integration difficulties in controlling spurious oscillations in the solution.
algorithms. The non-linear effects involving the dynamic An obvious remedy is to assume artificial properties (fictitious
equations have been traditionally treated by a variant of the mass and damping) which assure strongly damped responses.
Newton-Raphson iterative technique usually referred to as the In the context that the static response is the limiting case of
modified Newton-Raphson (MNR) method [5]. the damped dynamic response over a long period of time then
A time-discretized, incremental-iterative form of the an important parameter controlling convergence is the damping
equations of motion is written as coefficient. This parameter is determined using basic physical
** (k) * (k) (k -1) (k) (k) (k -1)
principles. Preference is given to values of damping just less
MUn+1 + CUn+1 + (K T DU )n+1 = Rn+1 - F(U n+1 ) (10) than critical so that the solution will oscillate about its
with equilibrium position rather than convergence from one side.
(k) (k -1) (k)
Un+1 = U n+1 + DUn+1 (11) The following sections will summarize the basic
formulations of the explicit and implicit DR methods. It should
In the equation, KT is the tangent stiffness matrix and k is
be stressed that in both cases the basic formulations follows
the iteration index. The static equations are obtained by closely the presented by, respectively, Underwood [4] and
dropping the inertia and damping terms, see equation (2). Wu [6]. Therefore, the contribution of this work relative to
However, occurrence of a singular or highly ill-conditioned KT these earlier publications relies mainly in some (crucial) aspects
will result in collapse of the solution for equation (2). regarding the solution strategy, including for instance the order
By making use of Newmark’s implicit integration formulas of the application of load components; the selection of
* (k) * ** ** (k) parameters, and the choice of the initial mesh, as will be shown
Un+1 = Un + [(1 - γ)Un + γUn+1]h (12)
later and in the numerical examples.
(k) (k)* ** (k)
**
Un+1 = Un + hU n + [(0.5 - β)Un + βUn+1] h2 (13) Explicit Dynamic Relaxation Method
the equation (10) becomes This section describes the basic formulation of the explicit
(k -1) (k) (k) dynamic relaxation method [4]. The fundamental time step
K TE DU |n+1 = DRn+1 (14)
equations will be derived first and then the equations used for
with
evaluating the integration parameters will be presented.
(k -1) (k -1)
K TE |n+1 = (K T + a0M + a1C)|n+1 (15) Derivation of time step equations
(k) (k) (k -1)
DRn+1 = Rn+1 - (KU) n+1 The algorithm of DR is based on the following semi-discrete
(k -1) * ** equations of motion governing structural dynamic response for
+ M[a0(Un - U n+1 ) + a2Un + a3Un] (16) the nth time increment
(k -1) * ** ** *
+ C[a1(Un - U n+1 ) + a 4U n + a 5U n] MUn + cMUn + F(Un) = Rn (19)
(k) (k -1) (k)
Un+1 = U n+1 + DUn+1 (17) where M is a mass matrix, c is a damping coefficient for
In the above equations, as are constants which depend only mass-proportional damping, F is the internal force, and R is a
** *

on the integration parameters, γ and β, and time step, h. vector of external loads. The vectors U, U and U represent the
Equation (15) shows that although the tangent stiffness matrix, acceleration, velocity and displacement vectors, respectively.
KT, may be ill-conditioned, presence of the dynamic terms Internal force F is a function of the displacement and may be
assembled on an element-by-element basis.
prevents the resultant tangent stiffness matrix, KTE, from
A diagonal mass matrix obtained by mass lumping, as well
becoming ill-conditioned. Equilibrium iteration on equation
as diagonal mass-proportional damping matrix are used in this
(14) are initiated at each time step until we achieve
(k) (k -1) research. This approach has significant computational
|| DUn+1 || ≤ eU ; || DR n+1 || ≤ eR (18) advantages. It leads to an uncoupled system of algebraic
DYNAMIC RELAXATION METHOD equations in which each solution component may be computed
independently. The need for assembly and factorization of
A procedure for moving towards a correct solution from a
global matrices is therefore avoided. The lumped mass
nonequilibrium starting guess involves the use of dynamic
technique also improves performance of the algorithm.
equations. The starting guess with zero velocities is taken as

4 Copyright © 2006 by ASME

 The solution of equation (19) is obtained using an explicit ~
time integration method. In this work a central-difference e = || R - Fn || / || R || = || Rn || / || R || ≤ etol (27)
scheme is utilized, where velocities are defined at the mid-point ~
of the time step, and the approximations for the temporal where R is the static load and Rn is the residual force vector
derivatives are given [7]. Namely, for time step n. Note that when convergence is obtained (i.e.
* 1 e ≤ etol), the internal forces balance the external forces and the
Un+1/2 = h (Un+1 - Un) (20)
inertial terms vanish – just as in a static analysis. Alternate
1 * convergence criteria can be formulated based on energy
Un = h (Un+1/2 - Un-1/2)
** *
(21)
considerations or when the norms of the velocity and
where h is a fixed “time” increment. Given that the velocity acceleration vectors approach and remain near zero over
vectors are computed at half-station points, it remains to define several time steps.
*
an expression for Un whish is needed for the damping term. In Evaluation of DR integration parameters
*
[8] it is recommended to average Un over a time step as The DR integration parameters for static analysis consist of the
1 diagonal mass matrix M, damping coefficient c and time step h.
Un = 2 (Un+1/2 + Un-1/2)
* * *
(22)
Performance (rate of convergence) is a function of the spectral
Substituting equations (21) and (22) into equation (19) and radius κ. The most rapid convergence is obtained for the
re-arranging terms in equation (20) yields the fundamental time smallest possible |κ| < 1 which produces uniform convergence
marching equations for advancing the velocity and over the entire range of eigenvalues l0 < l < lm. The following
displacement vectors to the next time step. Thus, expressions for h and c satisfy the optimal convergence
* (2 - ch) * 2h -1 condition [4]
Un+1/2 = (2 + ch) Un-1/2 + (2 + ch) M (Rn - Fn) (23)
*
h≤2/ lm (28)
Un+1 = Un + hUn+1/2 (24)
c @ 2 l0 (29)
where the matrix inverse of M is a trivial computation since
As noted in [4], the first equation corresponds to the
M is diagonal. The above equations involve the use of Un and
* stability limit of the central-difference integrator and the
Un-1/2, therefore calculation of the solution for the first time second equation represents critical damping for the lowest
step requires the use of a special starting procedure. Since U0 eigenvalue. Since l0 and lm are generally not known a priori,
*
and U0 are known initial displacements and velocities, the effective techniques for their estimation must be established.
following relation may be determined from equation (22) An upper bound to the maximum eigenvalue can be obtained
* * *
U1/2 = 2U0 - U-1/2 (25) from Gerschgorin’s theorem which implies that
~|
n |K
Making use of above expression, equation (23) may be
lm < max(i) ∑ M
ij
(30)
simplified to the following special form needed for the first ii
j=1
time step
* 1 * h -1 where K ~ is an element of the global tangent stiffness matrix. It
U1/2 = 2(2 - ch) U0 + 2 M (R0 - F0)
ij
(26)
may be noted from equations (30) and (28) that the time step
The following may be noted in reference to the fundamental size and the mass are not independent. Either of these
time step equations (23) and (24). These equations represent a parameters (h or c) may be fixed and then the appropriate value
linear combination of vectors for computing displacements at of the other determined. In the current algorithm, the time step
the next time step implying that each displacement component size is held constant and arbitrarily assigned a value of one. By
is uncoupled and may be computed independently. It also equating the equations (30) and (28), a relationship is obtained
implies that the computational cost per step is very low and is for defining the diagonal entries of the mass matrix as
2
mostly associated with evaluation of the internal force vector F. h n ~
Since F may be assembled on an element-by-element basis (so Mii ≥ 4 ∑ |Kij| (31)
that the need for a global stiffness matrix is avoided), memory j=1

requirements are very low. An estimate of the minimum eigenvalue may be obtained
The objective of a static analysis using DR is to obtain the from the following mass-stiffness Rayleigh quotient
steady-state solution of the pseudo-transient response. The l0 @ (Wn)T K*n Wn / (Wn)T MWn (32)
mass and damping parameters need not represent the physical where Wn is a selected weighting matrix. For linear
system. Instead, they are defined so as to produce the most
rapid convergence to a steady-state solution, where problems K*n is equivalent to the linear stiffness matrix K0. For
convergence herein is based on a relative error of the force non-linear problems, (K*n)ii represent the diagonal estimators of
imbalance or the directional stiffness after step n and are given by

5 Copyright © 2006 by ASME

 efficiency. They are affected by many factors, such as the form
*
(K*n)ii = {(Fn)i - (Fn-1)i} / h(Un-1/2)i (33)
* and amount of artificial damping and its control procedure, the
Weighting vector options include Un , Un-1/2 and Rn. In this size of load and time increments. With stability being the
*
research Un-1/2 is utilized. This choice appears to produce the priority in mind, some form of adaptive control of these factors
most reliable results, and it also allows for cancellation of according to the level of dynamic response is necessary to
*
Un-1/2 in equation (33) which eliminates the possibility of zero improve efficiency. These are discussed in detail below.
in the denominator. Numerical stability
For non-linear problems which exhibit structural Numerical experiments indicated that both computational
instabilities (snap-through buckling applications, for example), damping and artificial damping were needed to achieve
the stiffness matrix may lose positive definiteness. When this convergence [9]. In addition, it was necessary to adjust the
occurs, the lowest eigenvalue, and therefore the argument of artificial damping and the time step size adaptively as time was
the square root in equation (29), may become negative. Under incremented.
these conditions, the damping coefficient c is set equal to zero Computational damping is deliberately introduced by
as recommended in [4]. Also, some strategies to update selecting γ > 0.5 in Newmark’s integration formulas. The
adaptively the dynamic relaxation parameters may be used [4]. choice β = 1 and γ = 3/2 provides the strongest possible high-
Implicit Dynamic Relaxation Method frequency dissipation. This achieves very poor accuracy in the
This section begins by presenting the basic formulation of low modes and thus it is not useful for transient analysis.
the implicit dynamic relaxation method (without theoretical However, the strong damping characteristics may be exploited
analysis), following the lines already discussed in Wu [6]. in quickly bringing a transient solution to equilibrium, as
Again, the idea of dynamic relaxation is discussed here in the desired.
context of utilizing the kinetic effect for the construction of an Choice of the artificial damping term, though arbitrary, is
efficient conditioning mechanism. This paper presents some important. In general, the assembled global matrices are very
new contributions relative to the earlier work of Wu [6], related narrowly banded. For simplicity it is only necessary to contrive
mainly to some aspects of the solution strategy, that were a diagonal damping matrix. For example, a damping term
presented in [9]. cii = 2m siimii could be used, with s and m being the diagonal
While in the explicit formulation the use of diagonal mass terms of the stiffness matrix and mass matrix, respectively.
and damping matrices is mandatory to maintain the explicit For a given set of γ > 0.5 and β > 0.25(0.5 + γ)2 in
(uncoupled) nature of the algorithm, now in the implicit Newmark’s integration, the present numerical procedure
formulation there is no need to employ contrived diagonal, becomes an implicit single step scheme with two principal
lumped matrices, and the true structural mass matrix can be parameters, namely, m and h. Control of the artificial damping
used instead. In both explicit and implicit relaxation term was left to the adjustment of m. In general this value
formulations the damping term is artificial; while in the explicit should be chosen such that with the adjustment of load level
formulation the choice of a mass-proportional damping C = cM and time step, the number of equilibrium iteration on DU(k) is
is natural, for the implicit formulation, where the mass matrix kept below 10. In practice, only two or three iterations are
does not need to be diagonal, an artificial diagonal damping necessary. Adjusting the time step size affects both the inertia
term is added that does not involve the product of a coefficient effect and damping level.
by the mass matrix, as will be described below. The equations
of motion are solved by using the Newmark’s integration Acceleration of convergence
formulas in which computational (numerical) damping is also For improved efficiency of the numerical algorithm, overall
deliberately introduced. For stability and acceleration of convergence can be accelerated by adaptively adjusting the
convergence, detailed control procedures for the damping level, following parameters.
load increment, time steps, and the integration parameters are Load increment. For non-linear analysis, better stability is
presented. obtained by applying the external loads and given displacement
Relaxation of kinetic effect is achieved gradually by incrementally. Also, it is preferred that at each load level
adaptively decreasing the artificial damping coefficient and convergence on equilibrium iterations can be achieved very
increasing the time increment. The final solution is obtained quickly. The two parameters associated with load increment are
when at full load level the motion, force residual, and increase the size of load increment and the moment at which load level
in displacements are all zero within the specified tolerance is incremented. A constant load increment can be used, but
level. convergence can be greatly accelerated by adaptively
Theoretically, the transient response of the damped motion controlling the increment. This adaptive control, as proposed
will die down and the final static solution will be achieved if a by Wu [6], is done in such a way that the load increment is: (a)
sufficiently long time is allowed. Practically, however, this is unchanged if the number of equilibrium iterations (Niter) to
(k)
not always possible. There are the questions of stability and achieve convergence on DUn+1 is between 5 and 10; and (b)

6 Copyright © 2006 by ASME

increased/decreased by a factor of 10/Niter (maximum 2) The compound effect of adjusting both the time step size
otherwise. Quick convergence was achieved with a load and damping level is reflected by the steady decrease of the
increment varying from 1 to 10%. force residual. Convergence to the static solution is achieved
when the norms of displacement increment between steps n and
Significant advantages are obtained by applying the full
n + 1, the velocity and acceleration vectors and force residual
dead weight before starting the external load increments, as
were all zero within the specified tolerance level.
discussed in [9]. Also, the dead weight is not applied
incrementally; the full dead weight is applied in one step only. The damping parameter is also adjusted by the step-back
After that, the first load increment is started and the load level strategy described above [9].
is increased with time successively until full load level is IMPLEMENTATION AND NUMERICAL EXAMPLES
reached. That is, the load level is increased at the end of The proposed techniques have been incorporated into a
equilibrium iteration at the current time step. Once full load computer program for the nonlinear static and dynamic analysis
level is achieved, it is held constant over all the subsequent of lines [10]. The program is integrated to a pre and pos-
time steps until the final static solution is obtained. processing interface that generates the models (including the
Time step. The size of time step has a bearing on the level FE meshes) and visualizes the results. The integration of the
of damping and more importantly on inertia effect. The speed interface with the DRM method allows a completely automated
at which equilibrium iteration converges at each time step is generation of the FE meshes that correspond to an initial stable
directly affected. Presently, it was adjusted in two different static configuration of the system.
ways according to the load level. Before reaching the full load The user needs only to define length and properties of the
level, it was adjusted the same way as the load increment. Once line segments, and the position of the connections. All other
full load level was reached, it was adjusted as follows: it was parameters (including analysis parameters such as time step h
increased by 35% if Niter ≤ 5; it was unchanged if 5 < Niter ≤ 10; and damping coefficient c) are automatically selected and
it was decreased by 15% µ (Niter/10)1/2 if Niter > 10; and adaptively varied. Some details of this automated generation
notwithstanding, it was reduced by 25% if || Un+1 || > 2 || Un ||. procedure will be illustrated in the following examples.
The objective is to ensure convergence on the equilibrium Several small preliminary problems have been run to test
iteration and a reasonable rate at which this convergence is the validity of the theoretical predictions [9]. A variety of
achieved at each time step, minimizing the number of time examples involving complex configurations and nonlinear
steps and the total number of equilibrium iterations needed to boundary conditions were also analyzed. In each case, several
increase the time step size as the motions die down and the initial starting positions were used, including initial unstretched
solution moves towards the final static equilibrium position. positions. These unstretched positions can be conveniently
The step by step movement can be monitored by checking the specified arbitrarily; this advantage cannot be taken in
norm of the force residual || Rn+1 - (KU)n+1 || at the end of each conventional numerical techniques because the initial tangent
time step, say its maximum absolute value. As the motions die stiffness matrix is singular. For more complex mooring
down this should be reduced gradually. systems, the facility of generating the initial mesh becomes
The time step is a very important parameter in spite of the more prominent. Convergence was achieved in all cases,
analysis not being a dynamic analysis. The choice of the initial demonstrating satisfactory numerical robustness of the method.
time step value and its control procedure not always guarantee The efficiency varied with the quality of the initial starting
the convergence. The time step h = 0.1 was found to give position, as expected. For convenience, a straight line (free of
satisfactory convergence speed [9]. In addition, another loads) was used as initial mesh and then its connections are
strategy is incorporated [9], involving a step-back procedure moved to the design positions. This approach has given a
consisting in storing in memory the information of the previous satisfactory numerical efficiency for all cases.
step (such as nodal displacement and element curvatures); if the Offloading Floating Hose
current step fails to converge, that information is recovered and Among the studies presented in [9], a special case is a
the procedure restarts with a smaller time step. floating offloading hose, composed by materials of different
Damping. The damping parameter, m, is adjusted properties, with flexural stiffness relatively more important and
simultaneously with the size of time step to control the level of nonlinear hydrostatic behavior. Offloading hoses are pipes used
artificial damping introduced. Before full load was reached, m to transport fluids resultant of the explotation process, from a
was kept the same from one time step to the next. Once full floating offshore system to another, in general from a storage
load level is achieved, m is adjusted from one time step to the unit to a transport unit [11], as shown in Figure 1.
** **
next as follows. It was reduced by 25% if || Un+1 || > || Un ||. The floating hose treated here has a total length of 278.5m
This reduction is disallowed if || Un+1 || > || Un ||. If, however, connecting two vessels distant 250m. The hose comprises 53
Niter > 10, it is increased to 5n (n = Niter/10) times of its current segments of different diameters and materials, presenting
value for the next time step instead. segments with flotation. Their connections with the FPSO are
articulated. The FE discretization employs nonlinear frame

7 Copyright © 2006 by ASME

elements based in a co-rotational formulation that allows the the line leave the XZ plane. The final configuration assumed by
consideration of the flexional stiffness of the hose. The the line is shown in Figure 5.
environmental load consists of a current profile with 1 m/s at
the water surface and null at the sea bottom.

Figure 4 – Static stable configuration (without current

and initial imperfection)

Figure 1 – Offloading floating hose scheme Figure 5 – Static stable configuration (without current,
The pre-processing interface automatically generates an with initial imperfection)
arbitrated initial configuration that consists of a straight line In summary, three static stable configurations were
resting in the water surface. The initial mesh that corresponds obtained: the first one under the action of current, and the other
to this initial configuration is completely free of loads, without two without current.
initial tensions or curvatures. Prescribed displacements are then When current load is not applied, the line remains in the XZ
automatically applied to move their connections to the user- plane. That configuration, besides generating compression in
specified positions, as shown in Figure 2. At this time, some the line, generates curvatures and excessive moments in the
options are available, such as, to apply or not the current load; proximities of their connections. The convergence of DRM is
to use adaptive strategies in the solution procedure, or to hindered and a small loss of efficiency is seen in the solution
manually specify the analysis parameters. procedure.
The introduction of the imperfection in the initial straight
mesh caused the line to leave the XZ plane, forming an "S"
shape in the XY plane. In that configuration, compression is
practically eliminated and the curvatures are reduced. However,
this "S" configuration is not adequate for the subsequent
Figure 2 – Initial mesh dynamic analysis that will consider all environmental loadings
The first result to be presented is the final configuration including current. The significant flexural stiffness hinders the
under action of self weight, buoyancy and current, shown in line to assume its actual position under current action, similar
Figure 3. to the one shown in Figure 3, determined by the DRM with the
application of current.
When current is applied, the generation of the static
configuration is simpler, the solution procedure is more
efficient and the final configuration is adequate for the dynamic
analysis to be done subsequently.
Lazy-S Riser
Figure 3 – Static stable configuration (with current)
Next, two alternative strategies are used to obtain the stable Figure 6 presents a typical model of a Lazy-S riser
configuration of the line without application of current. In this configuration. This configuration presents an intermediate
case, the effect of introducing initial arbitrary imperfections is section passing through an arch with floaters (whose buoyancy
studied. alleviates the weight supported by the floating unit, and
contributes with restoring moments under lateral loads). It also
Firstly, Figure 4 presents the static configuration assumed includes a tendon that supports the buoyancy of the floaters.
by the line without the application of any imperfection to the
initial mesh. It is seen that the line remains in the vertical XZ The Lazy-S model can be seen as three interconnected
plane. lines: a catenary line connecting the platform to the floaters
(with a length of 100.0 m), a tendon (with a length of 71.1 m),
Then, another analysis is performed imposing a small and catenary line connecting the floaters to the seabottom,
imperfection to the straight mesh (a sine wave, with amplitude (120.0 m length). The discretization of the lines employs a FE
1x10-8m in the horizontal XY plane). That imperfection makes mesh with 176 nodes and 177 elements. The environmental

8 Copyright © 2006 by ASME

load consists of a current profile with speed 1.12m/s (to north) The mooring line considered here is comprised by 14
at the deep water surface and 0.57m/s (to northeast) at the sea segments of different diameters and materials, including chain
bottom (85m). and polyester cable, with a total length of 1662.0m. The FE
discretization of the line employs a mesh with 674 nodes and
675 truss elements. The following Table presents the
environmental conditions applied on the mooring line.
Table 1 – Environmental conditions
Current Profile
Depth (m) Speed (m/s) Going to Azimuth (degrees)
0 1.53 S 180
100 1.46 S 180
350 0.89 SW 225
Figure 6 – Typical Lazy-S Riser Configuration and Initial mesh 500 0.78 N 0
In this application, the strategy for the generation of the 1000 0.49 N 0
static equilibrium configuration does not consider initially 1260.9 0.00 N 0
straight meshes for the upper and lower catenary lines; rather, While the application of the NRM to this model was not
the initial meshes for the DRM method are defined as follows: successful, due to lack of convergence, both DRM formulations
a) The mesh of the tendon is defined as a vertical line; b) The (implicit and explicit) were capable to reach the static stable
mesh of the upper and lower catenary lines are defined configuration with equivalent accuracy. The static
employing the catenary equations (under self weight only), configuration is shown in Figure 9, along with the initial mesh
assuming all ends fixed. Therefore, the three lines in the initial (derived by the classic catenary equations).
mesh are equilibrated individually, but the whole system is not.
The static stable configuration determined by the DRM, is
shown in Figure 7. This configuration can then be employed as
the initial mesh of the subsequent dynamic analysis.

Figure 9 – Static Stable Configuration and Initial Mesh

Figure 7 – Static stable configuration An additional test is then performed, where the initial mesh
is determined not by the analytical catenary formulation, but by
Installation of a Mooring Line the use of the DRM in a FE analysis, without application of
The Lazy-S application of the previous section could be current, from a straight mesh as shown in Figure 10. This
amenable to both the DRM and the Newton-Raphson method straight mesh is generated by placing the buoy in the direction
(NRM). The accuracy is similar for both methods, and the of the azimuth of the line, with the distance between the buoy
computational efficiency is also similar. In this section, we now and the anchor equals to the line length, thus resulting in a
present an application where the NRM fails to converge due to straight line linking the anchor and the buoy. The analysis then
stiffness matrix singularity – the simulation of a step of the consists in releasing the buoy until the system finds its
installation procedure of a mooring line, where it remains in equilibrium position under dead weight only. In that case, as
position connected only to a buoy (Figure 8). For illustration expected, the equilibrium position determined by the DRM
purposes, this Figure also shows the installed configuration. coincides with position calculated by the catenary equations.

Figure 10 –Straight Initial Mesh

Figure 8 – Mooring Line

9 Copyright © 2006 by ASME

FINAL REMARKS In summary, the presented procedures based in the DRM to
This work had as objective to demonstrate the use of generate initial static stable configurations for flexible lines
dynamic relaxation algorithms for generation of initial static were shown to be quite efficient and robust, and comprises an
stable configurations for flexible lines. The use of dynamic important contribution to the analysis and design of offshore
relaxation algorithms not only overcame the limitations systems with flexible lines such as mooring lines, risers and
imposed by the classic catenary equations to the generation of hoses.
initial configurations of lines, but also avoided the numeric REFERENCES
problems that arise in the traditional FE solution methods for [1] CARDOSO, D. C. T., LIMA, B. S. L. P., JACOB, B. P.,
nonlinear problems, by supplying a conditioning mechanism ALBRECHT, C. H., MASETTI, I. Q., “Analytical-
for the tangent stiffness matrix. numerical methods for calculation of equilibrated
Two formulations of the Dynamic Relaxation Method were configurations of mooring lines with relief buoys” (in
presented, one implicit and another explicit. Both formulations Portuguese), Annals of XXXI South American journey of
have presented good results in the numeric experiments. In Structural Engineering, vol.3, pp. 1081-1090, may 17-21 ,
general, the following observations can be made: Mendoza - Argentina, 2004.
(a) It is extremely important to start the dynamic analysis [2] FERREIRA, F. M. G., SILVEIRA, E. S. S., MASETTI, I.
from a static stable configuration under the action of all static Q., MENEZES, I. F. M., “A strategy for the generation of
components of the loads, including current. The performance 3D initial configurations for dynamic analysis of morring
gain in the dynamic analysis is considerable, due to the lines” (in Portuguese), XXV CILAMCE – Iberian-Latin
reduction of the unnecessary transient part of the solution. Americam Congress on Computational Methods in
Engineering, Recife, Brazil, November/2004.
(b) The main advantage of DRM is its robustness, being
[3] OAKLEY, D. R., KNIGHT Jr., N. F., “Adaptive Dynamic
able to reach the static equilibrium configuration in situations
Relaxation Algorithm for Non-Linear Hyperelastic
where the Newton-Raphson Method fails (as is usual in
Structures – Part I: Formulation”, Comput. Methods Appl.
complex configuration of flexible lines). Although in some
Mech. Engrg., vol.126, pp. 67-89, 1995.
cases DRM may request a little more CPU time than the NRM,
[4] UNDERWOOD, P. G., “Dynamic Relaxation”,
the cost of the static solution is still low, insignificant in
Computational Methods for Transient Dynamic Analysis,
comparison with the total time of the dynamic analysis.
eds. T. Belytshko e T.J.R. Hughes, North Holland,
(c) The use of adaptive strategies for definition of time Amsterdam, pp. 246-265, 1983.
steps increases significantly the efficiency of the DRM [5] JACOB, B. P., EBECKEN, N. F. F., “An Optimized
algorithm. Implementation of the Newmark/Newton-Raphson
(d) The DRM presents an additional advantage in relation to Algorithm for the Time Integration of Nonlinear
the classic methods of solution for nonlinear problems: it is Problems”, Comm. Numer. Methods Engrg.,vol.10, pp.
very easy to implement. 983-992, John Wiley & Sons, UK/USA, 1994.
The implementation of the DRM in the graphic pre- [6] SHUKAI WU, “Adaptive Dynamic Relaxation Technique
processing interface (as a tool for generation of FE meshes of for Static Analysis of Catenary Mooring”, Marine
initial static configurations) leads to an automated model Structures,vol.8, pp. 585-599, 1995.
generation procedure that is completely independent of all [7] PARK, K.C., “Practical Aspects of Numerical Time
subsequent static and/or dynamic analyses that can be Integration”, Comput. Struct.,vol.7, pp. 343-353, 1977.
performed with the generated model (since it incorporates the [8] PARK, K.C., UNDERWOOD, P.G., “A Variable-Step
applicable initial state of tensions and curvatures). Therefore, Central Difference Method o Structural Dynamics
once the initial configuration/mesh is defined, several analyses Analysis – Part 1: Theoretical Aspects”, Comput.
can be performed with the same model. Methods Appl. Mech. Engrg., vol.22, pp. 241-258, 1980.
It should be stressed that, while the current implementation [9] SILVA, D. M. L., “Generation of initial stable
of the DRM is performed in an in-house software code, it could configurations of flexible lines by Dynamic Relaxation
also be used as a pre-processor module to generate FE meshes Methods” (in Portuguese). M.Sc thesis, COPPE/UFRJ,
of initial stable configurations for commercial FE codes. Rio de Janeiro, RJ, Brazil, 2005.
[10] JACOB, B. P., “PROSIM: Numerical Simulation of
Also, the time spent to generate this configuration is very
Moored Floating Units, Data input manual – version 3.0”
small: the generation procedure is automatically performed by
(in Portuguese), COPPE/UFRJ, Civil Engineering
the interface in a matter of seconds, at the click of a button. On
Department, Rio de Janeiro, Brazil, 2006.
the other hand, the performance gain in the dynamic analysis is
[11] COSTA, A. P. S., ROLO, L.F.A., GOULART, M.P.,
considerable when it starts from a static stable configuration,
SILVA, S.H.S.C., Offshore Loading Trends In Brazil,
due to the reduction of the unnecessary transient part of the
World Maritime Technology Conference, 2002.
solution.

10 Copyright © 2006 by ASME