Response modelling of Marine Risers and Pipelines

(1)

Trondheim, December 1990

DIVISION OF MARINE STRUCTURES THE NORWEGIAN INSTITUTE OF TECHNOLOGY

THI_ UNIVERSITY OF TRONDHEIM, iN vRWAY

OF

MARINE RISERS AND PIPELINES

TECHNISCHE UNIVERSITEIT Laboratorium voor Scheepshydromechanica Archief Mekelweg 2, 2628 CD Delft by Tel.: 015 - 786873 - Fax: 015 - 781836 Carl M. Larsen

(2)

ZNIN

do

ACKNOWLEDGEMENT

These lecture notes were partly prepared for the MOT project involving lectures in Panvel, India July 1989, and partly for the WEGEMT School in Trondheim January 1991.

The notes are written on the basis of several sources.

In most cases no references are given as the sources are Unaccessible (project reports. lecture notes in Norwegian etc.). The authors of these sources are hereby acknowledged:

Odd M. Faltinsen (hydrodynamic forces) Alf G. Engseth (3-D beam element) Bernt J. Leira (linearization techniques)

Trondheim, January 1991

itaiesce_A.

Carl M. Larsen

--adu-duadidf "To, dddMod

(3)

INTRODUCTION

In the following common features of marine risers and pipelines are listed and commented with regard to response modelling.

Marine risers and pipelines are long slender structures where cross section properties often are constant in long sections. The structures are typically one-dimensional without

branching, and the axial force - often referred to as the

effective tension - is the

governing stiffness parameter.

Figure 1.1 Riser and pipeline structures

(6)

-Response Modelling of Marine Risers and Pipelines

Introduction 1.2

Key-words to characterize structural behaviour are large deformations

small (elastic) strains

geometric stiffness more important than elastic stiffness (lateral direction)

Fig. 1.2 illustrates the principal difference between elastic and geometric stiffness.

I

4'

re,

4_11

C.

Elastic stiffness: Geometric stiffness:

Equilibrium obtained by change Equilibrium obtained by change

of elastic stresses of geometry.

Figure 1.2 Illustration of elastic vs. geometric stiffness

External loads from waves and current can be described by Morion's equation

1 TED 2

F = _5_.p CD DivIv + pCm a

CD Drag coefficient

Cm - Inertia coefficient, Cm = 2.0 for circular cylinder when X/D large

: Wave length

D : Diameter of member

: fluid velocity

(7)

Figure 1.3 Wave load parameters

Due to the large flexibility of these structures, the velocity parameter found in Eq. (1) must be implemented as the relative fluid-structure velocity, which means that a load non-linearity will be present.

Natural frequencies for risers and pipelines will in many cases be found in the frequency range for wave forces and also vortex shedding. Consequently, dynamic analysis is required for prediction of wave response, and modelling of correct damping is important. Hydroelastic phenomena like lock-in vibrations can also occur.

Risers and pipelines are found in a marine environment and hence subjected to large hydrostatic pressure, buoyancy and loads from waves and current. Methods for load effect analysis must therefore be able to include these effects in a correct and efficient way. In many cases, stochastic analysis is wanted.

The finite element method (I-EM) has shown to be a general and efficient tool for a wide range of structural applications. Other methods like transfer matrices, finite differences and more analytical formulations have proved to be efficient for some applications, but will in

(8)

Introduction 1.4

most cases be less general. From a FEM point of view, these problems seem to be rather simple, often applying beam or bar elements in one single chain.

In the following sections the most important features of the actual beam elements will be outlined, and some other items with respect to finite element modelling will be discussed.

(9)

2 FINITE ELEMENT FORMULATION

2.1

Introduction

The following notations will be applied: Stress-strain relation

a = C

E

G : strain vector

C :

Material law, matrix

E : Strain vector

Internal work, linear elastic system:

f T

U = - E adV

v T

E C EdV

v

Displacement at

arbitrary position in an element is

given as function of node

displacement

(2.1)

(2.2)

(2.3) I

(10)

Response Modelling of marine Risers and Pipelines

Finite Element Formulation 2.2

u = N v

(2.4)

ii : displacements at arbitrary position, vector

N :

interpolation functions

v : node displacement vector

Strains will be given by differentiation of displacement functions

E = Au = ANy By (2.3)

A : differentiation operand

Differentials of interpolation functions

Internal work (one element):

U =

I I

T11 TC B v dV (2.6)

-= T v (2.3)

TBT C B (IV (2.8)

Vi

Ki : element stiffness matrix, element i

The global stiffness matrix can then formally be written

K

= E

_aiT _{TiT a}

(2.,)

a : connectivity matrix

Ti:

Transformation matrices Ki : element stiffness matrices

K.

2

(11)

2.2 Two-dimensional beam element

The elastic stiffness matrix for a 2-dimensional beam element is shown below.

Figure 2.1 Linear beam element

EA 0 12E1 0 symmetric Q3 6E1 4E1 c2 -EA

-12E1 6E1 0 12E1

c2

6E1

0 2E1 0 6E1 4E1

A

The second order strains can be expressed as

EA

(2.10)

The geometric stiffness matrix takes into account that the structure will change its geometry when subjected to external load, and that this change will affects equilibrium (see Fig. 1.2). In order to derive this stiffnessmatrix, one have to consider the second order strains in the

beam.

g2

(12)

Finite Element Formulation. 2.4

=

1 _{all() 2}

1()

aV 2 E_xx

2 ax

,

'notations, see Fig. 2.2.,

Fig. 2.2 Beam element

The second term of Eq. 1(12) is the second order contribution to axial strains from lateral displacements, which is relevant for 'development of the geometric stiffness. This matrix is derived from the work form lateral displacements on axial stresses.: The strains of concern is therefore limited to

1

_{= =()}

1 av 2

xx

2a

By assuming the axial force P constant along the element and a linear variation of the bending moment, we have

P(x) P = eonst..

(2.13) M(x) = M 1(1. - + M2

Q'

Normal stresses are then seen to be

Gx(x,y) M(x)

A

The strain energy can theft be found:

IT (2.11) (2,12) (2. 1t)' + =

-=

(13)

- differentials of interpolation, functions or f_.

The higher order strains then seen to be

1 av 2

= --(---)

xx

_{2 ax}

I T = v N_'_x

NT

V

u =

yG dV

(2.11)

where both E and G have x-components only. The integration must however, be performed both over the beam cross section and along the beam. We will here include the second order strains only:

Li = Exx Gx(x,y) dV (2.16)

where Gx(x,Y) is given by Eq. (2.14) and Exxl by Eq. (2.12). Note that Exxl is

constant over the cross section as the lateral displacement v is constant according to Navier's

hypothesis.

The lateral displacements v can be described by

V = v N (2.17)

v - node displacement vector

N -

interpolation function Differentiation gives av _{= v} _N,x

(z./ g)

-1 t V T ax

(14)

U

44._

, 4-Lc_ Ca

2, 8

$

144..e.44

kg cs

: U = 1 lvTN G(x y)N Tv dV v KG = IN,x(3(x,y)N,x dV

When developing this integral, one can separate between the bending and axial terms in the stress function G(x, y). It can be proved that the axial term only is of interest for the present application. See also Felippa /1/. The geometric stiffness matrix can then be developed:

0 SYMM

K=

G

7

KTOT = KE KG

and transformed to the global coordinate system by

0 2e2 10 15 0 0

00

6 0

--0

=-5 10 Q 2g2 0 0 10 30 10 15

This matrix is constant as long as the axial force P is constant. Note that the matrix

represents linearization of the non-linear "rope effect'' as illustrated below, and that effects due to beam curvature is included in addition.

For application on the general 2-dimensional case, the matrices (2.10) and (2.22.) are added.

T.

(2.20)

(2.21)

(2.22)

(2.23)

0 6 4A,L, _r

-L

v 2 +

(15)

-Figure 2.3 "Rope effect"

...

Figure 2.4 Global coordinate system

K_Global,TOT

=T K .

_TOTTT (2.24)

In the general 3-dimensional case one have to include another type of transformation when dealing with large rotations. This will be discussed in the next section.

2.3 3-dimensional beam element

In order to handle the large displacements and rotations found in marine pipelines and risers, a beam element without any displacement limitations is needed. This can only be achieved by applying a so-called co-rotated ghost reference system, which will be described in the

following.

(16)

Fig. 2.5 shows how the 12 nodal degrees of freedom for the beam element are defined in relation to the local x, y, z system. The beam theory is based on the following assumptions:

During deformation a cross section of the beam remains plane and perpendicular to the

x-axis.

Lateral contraction caused by axial elongation is neglected. The strains are small.

Shear deformations due to lateral loading are neglected, but st. Venant torsion is accounted for. e z2 y2 vy2 vz2 Y0x2 X

Figure 2.5 Nodal degrees of freedom for beam element.

on

vx2

The stiffness matrix for this element can be derived in a similar way as for the 2-dimensional case. What is important in the 3-dimensional case, is the description of large rotations and

stress calculation.

Rotations in the 3-dimensional space needs to be handled with care, because large rotations 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 ii that is frozen to the nodal point and that follows the movement of the node. This coordinate system is parallel to the global base vectors Ii in the initial configuration (see Fig. 2.6). The orientation 4. the nodal point in space is uniquely defined by the base vector transformation.

(17)

(n) T(n) (2.25) i

where T is the rotation matrix containing direction cosines for the i vector relative to I.

initial position 7 1 n ) - (3n ) -.-(n ) 12 displayed

and rotated position

Figure 2.6 Nodal point with translational and rotational degrees of freedom.

The rotation of a node is thus not given by angles, but by the components of the rotation

matrix T which is orthogonal. In the current formulation the general motion of a nodal point is described by the three displacement components ui and the nine elements of the rotation

matrix Tij.

During the incremental solution procedure Tij is updated by treating the incremental rotational components about the global axes as rotations for small displacement problems.

Fig. 2.7 shows the beam element in deformed configuration with the co-rotated reference configuration Con close to C. In addition to the two coordinate systems with base vectors Ii and

i i, two more coordinate systems are used. As shown on the figure, these are

.0

represented by the base vectors ti and ii. The first one is a coordinate system that follows the deformation of the beam end. There is one such system for each beam end, however, only the coordinate system at end "b" is shown. The coordinate system with base vectors ii gives the orientation of the ghost reference configuration Con. Fig. 2.7 also illustrates that the nodes

I =

(18)

may be eccentrically connected to the beam ends. Note that the transformation between the beam end coordinate system and the corresponding nodal point system does not change during

deformation ii.

Figure 2.7 Deformed beam element

When calculating internal forces in the element, one makes use of the nodal displacements of the element referred to the local element system. It is then easily seen that

vxi vy1 vz1 = Vy2 Vz2 (2.26)

'?=

L - L0

x

where L denotes the secant length of the deformed beam and Lo is theinitial beam length. The rotations at beam end "b".(0x2, ey2, az2), are found by considering the difference in the ii°b-system and the it-system. A similar procedure give the rotations at beam end "a". When the element displacement vector referred to the local base vectors i1 has been determined, the strains, stresses, forces and stiffness may be calculated according to conventional beam theory. Non-linear material behaviour may be accounted for if wanted.

As a final remark it should be pointed out that the assumption of small deformations used in the stress analysis, does not represent any serious limitation. The formulation is still valid for large displacements and rotations since the reference configuration is rotated along with the

element.

b

2

node b

,ode a = =

(19)

2.4 Consistent vs. concentrated loads

Use of a concentrated force model on long, slender beam subjected to a distributed load is illustrated on Fig. 2.8.

c

(t--.7")

471-

4.,111,

J

rrox- .

Figure 2.8 Use of concentrated load

C:5 4,..07 Lire

CD Li. CI-4.-- c4-C_ cQ.

(."-SLCA-c-trO ct

to ct_ck. _`D'A

The error term 42/12 is seen to be proportional to the element length squared, and is hence strongly dependant on the actual discretization. This term should be compared to the term

M = EI

(2.27)

where El: bending stiffness curvature

For flexible pipe with small El, the moment from distributed element forces will be significant compared to Eq.(Z.Z1) even for short element length. Use of concentrated forces and direct calculation of moments from the stress matrices should therefore beavoided. For steel risers and pipelines, the El term is significantly larger, which means that the concentrated force model can give acceptable result with reasonable element lengths.

Note that the term 42/12 is connected to element properties. This means that this correction will be different for two adjacent beams if the element lengths or load intensities are different. On the other hand, if all elements in a model are of identical length and have almost uniform

ttO

(20)

load, the concentrated force model combined with a post-processing type of correction will give acceptable results. This strategy is implemented in the RISANA program for analysis of tensioned steel risers.

The correct modelling technique is however to apply a consistent load model. "Consistent loads" means that the load formulation is consistent with the stiffness matrix, which means that the nodal loads should be determined on the basis of the same interpolation functions as the stiffness matrix.

Fig.2.9 shows an element subjected to a distributed load q. In the finite element method this load must be represented by a nodal force vectorP.

As seen on the figure, the distributed load is represented by nodal forces and moments. The consistent load vector P can be calculated by

F.

Figure 2.9 Element and nodal loads

P =

JNT

Np dx (2- 1-8)

Notation: P: Consistent nodal loads

N: Interpolation functions for displacements, see Eq. (2.4)

0

(21)

The consistent load vector P is related to one element, and all contributions are added to the global load vector R. If two adjacent elements have identical length and distributed load, the moments in the load vector P will be equal but of opposite sign, which means that they cancel when added to the global load vector. This is illustrated in Fig. 2.10.

The contribution to the global load vector at the node between the two elements will be P = P12 + P21 = 2P12

(2.31)

M = M12 ± M21 = 0

This illustrates the strategy previously mentioned where corrections due to distributed loads were performed as post-processing only.

?

Np: Interpolation function for loads

Q: element length

load intensities at nodes

Note that the distributed load function :-s 1 (x) = N p(x) -

r,

Hence, the consistent load vector can also be expressed as

P

= I

N Tr) dX

i

where both N and ct are functions of the length coordinate x.

-(2.R)

(2.30)

t-t

ti,,,

rl

2,4

(22)

Pcz ?2- ---Ce

Figure 2.11 Consistent loads at node, different elements.

The figure above illustrates the general situation with two different elements. The stress analysis will be performed as follows:

Element 1: M12TOT MG12 4- M12 (2.34i) Element 2: M21 _TOT MG21 ÷ M21 (2.35) 11CS k Acyt,. E.1 c1/4,-t_ csc,,,,,Leca. L-0 t,

-tct L,g,

K r = R

(2.32)

we can find node displacements for all elements, and stress resultants in each element can be found by the relation

Si = k cv

(2.31)

where ki: element stiffness matrix including geometric stiffness

vi: node displacement vector Si: stress resultant vector

The calculated moments must now be corrected due to the effect of load distribution

Pcz

(23)

tvwo (2 01-4".44

t-Ttw-,JA

rtrr tA,C.4-,Cr

cos.4-e,

c41ro dA". ct_se. t:.% me ate

Fig. 2.12 illustrates the error introduced by using the uncorrected results from a concentrated force model as compared to the consistent force model.

For cases where non-linear material properties are included, use of consistent loads as described here is in principle not possible. This follows from the fact that local and global moments are superimposed according to Eq. (2.34) and (2.35). Moments should then be calculated from element curvature (interpolation functions) and not by use of stress matrices.

Concentrated loads will then give a good result.

Another important factor when dealing with stress analysis is the way axial stress is calculated

in a non-linear analysis. In Eq. (2.247) the axial deformation istaken as the difference between

the secant length for the deformed element and the element length in stress-free condition. (See Fig. 2.13a).

Alternatively, axial deformation is taken as the difference between the true curve length of the deformed element and the length in stress-free condition. (See Fig.16b).

The last alternative will introduce a coupling between bending deformations and axial strains, and through the geometric stiffness term (Eq. 2.22) also to the stiffness matrix. This coupling can in many cases give numerical problems, and should be avoided. For applications on marine pipelines and risers the first formulation of axial force calculation is recommended.

Use of hybrid elements has also been proposed, see McNamara /2/. The benefit is claimed to be more stable algorithms for non-linear dynamic analysis of problems with low or even zero tension, i.e. flexible risers on shallow water. This method introduces the axial force as an

(24)

Figure 2.13 Alternative axial elongations

2 o

0---

C&c,--c"

Pe. 4,7 a

Figure 2.12 Bending moment in wave zone, steel riser

-cCi

C.Ismiossoc===isagx,C)

-2-o

5 a)

0

(25)

2.5 Solution strategies

A complete analysis of marine risers and pipelines will normally consist of a static and a dynamic part. In the following, strategies for use of the finite element method and other methods will be outlined for pipelines during installation, anchorlines, tensioned marine risers and flexible risers.

2.5.1 Pipelines during installation

The static problem is geometrically strongly non-linear and can in principle be solved by using a non-linear finite element approach involving an incremental procedure starting from a stress-free configuration. This approach will, however, be unnecessarily time-consuming and should be substituted by a more efficient method, i.e.shooting method (iterative use of transfer matrices) or a stiffened catenary approach (see Sec. 6 and /3/). Other staic solutions are also available /4/. The dynamic analysis and eigenvalue calculation can, however, apply the finite element method. The following strategies are proposed:

A. Static/linear dynamic analysis

Calculate initial (2-Dim.) configuration by stiffened catenary theory.

Determine finite element model on the basis of geometry and effective tension along the pipe. This configuration will serve as a reference for all succeeding analyses, and later results will be superimposed on these deformations and stresses.

Rla.. A./....Z.,

/

Figure 2.14 Pipeline during installation

...I

,

2.,

(26)

Perform if wanted a geometrically nonlinear analysis by some incremental procedure to account for effect excluded in the initial analysis, i.e. current normal to the geometry

plan.

Perform a linear dynamic analysis by using the initial configuration or the configuration after step 3 as basis. This analysis might be performed in frequency domain or in time domain with constant system matrices but including non-linear loads.

B. Static/non-linear dynamic analysis

1. Calculate initial configuration by stiffened catenary theory.

"). Determine a finite element model on the basis of geometry and effective tension along

the pipe. As a part of this model. A "stress-free" configuration must be defined, i.e. the pipe resting on seafloor. This configuration must be known in order to perform a correct calculation of internal forces.

Perform iteration for equilibrium on the finite element model in the initial configuration. Note that the results from step 1 will not necessarily give equilibrium as the finite element model is not identical to the stiffened catenary model. The initial solution will, however be close to the correct FEM solution, which should make the equilibrium iteration fairly fast.

Perform non-linear static and dynamic analyses on the F.E. model. All results will refer to an equilibrium situation and hence give the complete stress picture without any

superposition.

2.5.7 Anchor lines

An efficient method for static analysis of anchor line including axial flexibility is given in ref.

/5/. This method should be used for the initial 3-Dim. calculation. For the FE model, bar

elements is sufficient as no significant bending stiffness exists. The strategy for analysis is further identical to the pipeline analysis.

(27)

Figure 2.15 Anchor line

2.5.3 Tensioned risers

The initial configuration of a tensionedriser is statically determined, which means that no sophisticated analysis is necessary. Axial stresses and strains, node positions and effective tension can be calculated by direct inspection of vertical equilibrium. The finite element model can then be established, and if a stressfree configuration is wanted, (non-linearcase) this can be identified as a weightless, non-tensioned vertical beam.

The dynamic analyses can then proceed along the same lines as for the pipeline case.

Figure 2.16 Tensioned riser

572

In. 0 U-k_cf;

4

(28)

2.5.4 Flexible risers

The initial configuration of a flexible riser can in many cases be rather difficult to calculate. This is specially true for riser configurations involving special anchoring devices, cradles or surface buoys. The most efficient strategy for this analysis seems to be the "shooting method" /6/, which is an iterative solution of the boundary value problem by transforming it to an initial value problem with one or more parameters as unknown. It is, however, difficult to make general computer programs based on this method, which means that a general FE based procedure should also be available. This procedure must start from a stressfree configuration and perform load and displacement increments until the wanted configuration is reached. This analysis might be very time-consuming and can in many cases have problems with regard to stability for intermediate positions. An example of the stress-free and final configuration is shown on Fig. 2.18. This particular analysis involved approximately 500 meters with forced displacements, using more than 100 increments.

For flexible risers material non-linearitiesilarge displacements and tension variation can be important, which means that a non-linear time domain procedure is necessary. Sophisticated

(29)

"fe Mck-f?-.1 C-4c-i la-LA

)\--coLtA LA, co oe,--&

c.

Figure 2.1g 1-EM analysis from stress-free configuration

computer programs based op the second strategy previously outlined has therefore 'been

developed,

An important detail when dealing with design of flexible riser ystern is to have a bending gtiffenet at ends that give a gradually increase of the bending stiffness from the flexible pipe to a -rigid connection. Aspects of design and analysis, procedures in this relation is found in Ref. /7/.

2,6 References,

/1/ Felippa, C.A..:

Refined Finite Element Analysis of Linear and Non-linear Two-dimensional,Structures

Report to national Science Foundation, U:S.A. Berkeley, California Oct. 1966, /2/ McNamara, J.F., O'Brien, P.J. and Gilroy, S.G..:

Non-linear Analysis of Flexible Risers Using Hybrid Finite Elements, OMAE, TOKYO

1986.

13/ Dixon, D.A. and Rutledge, D.R.:

Stiffened Catenary Calculations in Pipeline Laying Problem. Journal of Engineering for Industry, Feb. 1968.

(30)

/ 4/

`?. 1

CLLAC/A ok._

C

ç.c_ca_

151 Korkut, M.D. and Hembert, E.J.:

Some notes on Static Anchor Chain Curve Offshore Technology Conference, Houston

1970.

/6/ Keller, H.B.:

Numerical Methods for Two-point Boundary Value Problems. Blaisdell Publishing Company, 1968.

/7/ $zcZa...L..i.,

c-C

N-9 c.) tt9

,

Response Modelling of marine Risers and Pipelines Finite Element Formulation

2.22

(31)

3 HYDRODYNAMIC FORCES

3.1

Introduction

This section will not cover the subject "hydrodynamic forces" totally. Here, only implications on methods for analysis will be discussed.

The discussion will be limited to slender cylinders in waves and current, and the cylinder is assumed to move. ( Figure 3.1 Notation velocity of cylinder D: cylinder diameter vc: Current velocity vw: wave velocity aw: wave acceleration

ii(x,t): wave elevation

(32)

Hydrodynamic Forces 3.2

Calculation of hydrodynamic forces due to waves and current must always be performed in two steps:

Calculation of velocities and acceleration in the fluid

Calculation of forces on the cylinder from given velocities and accelerations

In the following, both aspects will be discussed.

3.2 Wave kinematics

Linear wave theory:

L /Z (x(4.)

Figure 3.2 Linear wave, longcrested

The linear (Airy) wave theory assumes that the gravity wave can be described by harmonic

functions. This means that sea surface elevation, pressure, velocity and acceleration always can be written as:

p.(t, x, z) = p.a - eKz sin(-wt + 1C.X +E)

NOTE: These sign conventions may vary! 4: parameter (pressure etc.)

p.a: amplitude of parameter at sea surface (mean water level)

(33)

K : wave number =

A,: wave length

co: wave frequency

phase angle

The wave length and the wave frequency is linked by the relation

A. gT2 tgh (27c D) 2n (3.2) T: wave period, T = 2m co water depth

Note that Eq. (3.2) must be solved iteratively as X appears on both sides. For deep water (D > X) one can neglect the last term, and the relation is simplified to

A. gT2 (3.3)

2n

Eq. (3.1) is the general expression for all parameters. Special formulas for vx, vz etc. can easily be derived, but will not be given here.

When applying this equation to the actual problem (see Fig. 3.1), we have to face the

following problems:

1) The structure penetrates the mean water level. This means that the loaded area will vary, due to the presence of wave crests and though as a function of time. The load at one specific point on the cylinder will therefore not be a harmonic function.

(34)

LA_A-Response Modelling of Marine Risers and Pipelines

The key assumption behind the linear wave theory is that the wave

amplitude is

infinitively small. This theory is therefore unable to predict velocities and accelerations in the wave crest.

2

Figure 3.3 Extrapolation

The structure itself will move in the wave potential, which means that the 1x term in Eq. (3.1) will not be constant in time. The velocities and accelerations experienced by the structure will therefore not be harmonic functions.

In practical analyses these effects will in many cases be neglected. Different methods for

analysis can, however, handle these problems differently. This will be discussed in Sect. 4

and 5. Here, some implications with regard to wave crest kinematics will be discussed.

There are many different techniques to overcome the "wave crest problem".

a) Extrapolation of potential

The technique is illustrated on Fig. 3.3. The function

LJA,

c;' 2)

(35)

f(z) = e

(3.4)Kz

is simply extrapolated above mean water level to predict values in the crest, and for situations with water below the mean value, the exposed area is simply reduced without any potential

modification.

Note that this method can lead to unrealistic high value for velocity and accelerations for steep waves, and that the method will predict a considerable mass transport due to waves. The method is impossible to apply for irregular waves, and should in general not be recommended (at least for non-steep waves only).

b) Use of surface value in the wave crest

CA-cA Figure 3.4 Use of mean water level

The difference between a) and b) is for the crest case only. This methodwill also involve some mass transport, but less than case a). The approach can be used for irregular waves, and unrealistic crest values are avoided. However, high frequency components might be given an unrealistic large extension along the riser when using this approach for large stochastic waves. The reason for this is that the expontential decay will not be active in the wave crest, which

introduces an error that will increase for shorter (highfrequency) waves.

c) Parallel move of potential

This approach will give a symmetric force integral with respect to directions, but a constantly moving load profile. No mass transport due to waves is present. Some modifications of the

(36)

(X-Figure 3.5 Parallel move of potential

profile must be performed at the bottom end for shallow water..

Stretched or modified potential

ot

e) Use of higher order wave theory Figure 3.6 Stretched or modified potential

This technique will use a constant potential at a given depth D, and also the instantaneous water level. The function in between these two points will be described by a simple linear coordinate transformation.

This method is more complicated than the "parallel move" approach, but gives probably a more correct mass transport.

(37)

Higher order wave theories like the 5th order Stoke is available and gives a more coffee./ description Of the gravity wave phenomenon as the boundary condition P is satisfied for the true water surface. It is, however, not possible to use this approach for irregular wave.

Figure 3.7 5th order Stoke wave profile

Note that the 5th order Stoke profile is unsymmetric and will involve some mass transport.

The computer program RISANA can handle smile of these methods and make it possible to compare results.

3.3 Inertia forces

,Inertia forces are linear, and can therefore always be decomposed without any problems. Fig. 18 shows a cross section of a structure subjected to waves and own motions:

The inertia forces Fy and Fz can then be expressed as

Fy dx[(pA Ay) aw,y Ayy 55)

(3.5)

F

z = dx[(pA + Azz) awz Azr

At Cross section area

Ayy, ,Azr Added thas-s, y- and:z-dir. Note that these might be different

-= 0

(38)

Crt, 5s

_C-c-rkso,

In a dynamic analysis the terms

A_yy, A_zz

will be included in the mass of the structure, while the other term will be introducedamong

the external loads. This separation is theoretically correct as long as the correct coefficient are used.

Often the following notation is used for the inertia forces:

Tri)2 ir.D2

H H

F

_{= p(Cm - 1)}

c

m _w

CM: Inertia coefficient = 2.0 for cylinder

DH: Hydrodynamic diameter

This formulation can lead to errors if the hydrodynamic diameter does not give the correct cross section area. By using the (3.5) formulation, one is always safe!

Another complication connected to inertia forces, is the orientation dependant hydrodynamic mass. This is illustrated in the figure below.

(3.6)

(3.7)

Figure 3.8 CO _C".14

(39)

-z

K

The local masses m, and my are defined in a coordinate system that follows

the pipe

orientation a. When representing these in a global coordinate systemX, Z we have

mx = my - cos a + mx sin a

_(3.8)

mz = my - sin a + mx cos a

For constant cross section, mx = 0. If the element rotates (a changes) the local masses will remain unchanged, but this is not the case for the global ones. This effect introduces a non-linearity to a finite element model of the system when large rotations occur.

3.4 Drag forces

Morison's equation gives a quadratic relation between force and velocity: FDRAG _1p CD Div riv ,

2

(3.9)

where vr is the relative structure-fluid velocity.

The quadratic term makes life rather complicated! Key-words are as follows:

Decomposing of problems with non-circular cross sections Combination of current and wave loads

Load dependant on response

Introduction of other load frequencies than the wave frequency

AAA. _-.,

0

'b"tr 2 /%'- -:-- .f. ,--t

-I)

(40)

Comments on these key-words will be given in the following

3-D problem, circular cross section and cross sections of arbitrary shape

1 R 2

FR = dr

p CD DRuy) +

Figure 3.10 Non-circular cross section

Such cross sections are often found for flexible riser systems with many flowlines integrated in one bundle. For this case the CD is not independent on direction, and use of Eq. (3.10) is therefore possible only if CD is known as a function of incoming flow direction. Instead of Figure 3.9 \b-A..0 c

jaR

R ' resultant of relative velocity normal to pipe

axis

FR : Force resultant, in same direction as uR

This procedure is straightforward as the drag coefficients for z and y directions are identical.

Problem occurs if a cross section has the shape as shown on Fig. 3.10.

re 1?

(41)

aiming at establishing this kind of function, one should have CD given for the two main directions z and y (Fig. 3.10) and calculate Fy and Fz independently.

R R

Fy =

p q Dy 1

uy 1 uy

Fz = 1 P cL Dz 1

u: 1 u:

(3.11)

Note that using these formulas for a circular cross section gives another result as Eq. (3.10).

Combination of current and wave load

Normally, current is regarded as a static load, while wave effects are dynamic. When using a static analysis followed by a linear dynamic analysis, one shall normally include dynamic

effects only in the dynamic analysis and forget about the

static condition.

This is

''remembered" by the static deformation and stresses, and dynamic results are simply

superimposed to the static. Due to the quadratic drag term, however, we must always include the current velocity when calculating dynamic forces. Strictly spoken, forces in the dynamic analysis shall always be the difference between total forces and forces present in the static

(42)

Hydrodynamic Forces

3J2

Figure 3.11 Submerged element

The element on Fig. 3.11 has a static force given by

n 1

F

= _ pC D D Ivc- sin a

Ivc - sina

s ₂

1 n n

Fn

_TOT =

7 _{pCD DvR}

Iv

1

R

v---`'`

In the dynamic case, the relative velocity normal to the pipe axis is given by

n

VR = Vc sina + Vy,w

sina + v

z,w - coscc - i-K.

^

(3.12)

(3.13)

The total force for the combined current, wave and structural response case will than be found

as

(3.14)

The force that should be introduced in the linear dynamic analysis will now be found as

F Drz YN = FTnOT - F'2

(3.15)

Note that the current effect can not be left out in this calculation.

(43)

(12)

c.--Figure 3.12 Wave zone loads

In the wave through (Fig. 3.12 (1)), there is no load in the static case. This means that when performing a linear dynamic analysis the static load must be subtracted at time increments for

which elements are found in a wave through.

1

F DYN=pCDDIV

TOT TOT

2

This can easily be performed in a time domain approach but is impossible to include in frequently domain.

In the wave crest (Fig. 3.12 (2)), there was no load present in the static case. Due to the wave, water will arrive, and the velocity in the crest must include the combined current and wave potential. Consequently, the dynamic load must be found by

R

where the currently velocity is included in V Top This is nothing else but application of Eq. (3.14) and (3.15) for the case that there are no static loads.

,

k.rto ./Ca

(3.17)

(44)

Load dependant on response

Drag forces must be calculated from the square of the relative velocity. Consequently, the drag force becomes a function of the response itself. Due to the square term, this can not be separated and introduced partly as a damping effect and partly on in the left hand side:

mY + c + kx = F(t) (3.18)

The dynamic equilibrium equation must however be written

mY + c + kx = F(t,X,x) (3.19)

The .t term in F comes from the quadratic drag effect, while the x term is due to the change in position in the wave potential during dynamic response. Geometric effects due to rotation are also present.

In the dynamic analysis the solution of Eq. (3.19) must be found by some kind of iteration. An alternative approach is to apply some kind of linearization techniques. These items will be outlined in the following sections.

Frequencies in quadratic drag forces

The function that describes the quadratic drag force reads

F(t) = k

sincot sinwt

2 = Fosin2C0t sign(sinwt), F = _1 p CD D uo

2

The graph of this figure is shown in Fig. 3.13.

(3.2C)

(45)

F(4)

Figure 3.13, Graph of quadratic drag force

Eq. (3.20) can be represented by a Fourier type series:

00

F(ot) E bsin not

1=1 It b =

IC

f(ot)sinnot

it it

21

= f(c)z)sfrincOt o 1

cr.D u2 sin2 cot - sinnom - dot_a

= p

- Ua

'CDD 2

n cut sinncot ,d(cot) it 0 ( (3.21) - -= p

(46)

etc.

Hence, the force can be written: F(wt) = Fo -2 [-4 sincot - 4

TE 3

4

sin 30)

-15 105

The non-linear drag term in Morison's equation including relative velocity reads:

q(t)

=

CD D j

-

1 (uw - i-)

(3.22)

(3.24)

This means that the drag force will have other frequencies represented than the frequencies directly found in the wave.

Note that for the current/wave case, load components will also appear for the even frequencies 20), 40) .. etc..

This effect can be important incases where eigenfrequencies are found close to 30), 5 co..etc.

3.5

Linearization, harmonic waves

Integrals:

TC

fsin2 cot

sin cot dot = 4

0

IC

sin2 cot

sin n cot dot = 0 for n even

0

sin2 cot sin Rot dcot = 4

0 15

rc

sin2cot

sin 5 cot dot =

-4

0 105

sin 5 0.)t (3.23)

(47)

where is the wave induced velocity

is the structure's velocity, both referred to the same direction.

This representation makes it impossible to split up the drag force into one damping term as a function of i- and one term for external forces given by uw only. Such a representation is, however, needed in frequency domain analyses. If such splitting should be possible, the force must be represented by the following type of function

q L(t) =

1p

CD D - KL (uw - t)

(3.25)

2

KL: linearization coefficient, to be derived in the following

For the harmonic wave case, the wave induced velocity will be representedby

LL ua OM (3.26)

The structure's response i- will then also be harmonic, but in general out of phase with the wave velocity:

r = r1 cos cot + r, sincot

(3.27)

and the velocity:

= cori sin Wt + cor2 cos cot (3.28)

The relative velocity can then be found:

uw -

(ua - or2)cos

+ cori sin om = A(cos cot + (3.29)

where A and 4). are easily found from vector analysis

(48)

A = \1a

-

(0r2)2 + 0)2r22

(3.30)

As (3.29) represents the relative velocity, the phase angle is no longer of interest and Eq. (3.29) can be written

i- = A coscot (3.31)

w

Note that Eq. (3.31) has the consequence that the load will appear at the same frequency as the wave. The 300, 56) etc. force components are hence neglected.

The non-linear and linearized drag force are

qNL = ___1pD CD A21 cos Wt I cosom

2

q L =

_pDCDKcos(0t

1

2

Note that K L* here is not identical to KL in Eq. (3.24). The error from the linearization is

then:

e

2

= (qL -

qNL)2 (3.34)

This error can now be integrated over one wave period, and KL in Eq. (3.33) can then be chosen so that the integrated error has a minimum value, formally

T

-2

f 2

e = e 0 - (3.35) T =

f1

2 2 21 * 2

D CD kL - A

1 cos (iv 4 0 1 12

and the condition for best possible KL is then

aE2

= o

afq

which gives an equation to determineKL.

COS2avdt

(3.32)

(3.33)

(3.36)

(49)

-By use of integration/differentiation rules we have

8A2

3m

KL*

The linearized drag force can then be written

1 8A 2

qL = 2PD CD

coscot

This force is the interaction force, partly damping and partly exiting, depending on time and phase angle. The force can be split into two harmonic forces, one damping and one exiting by introducing Eq. (3.31):

1

qL = _p D CD

LA (uw - i-)

2 3m

1 8A

= ___.p D CD U (external, exiting force) (3.39)

2 37t w

1

8A.

- _p D CD

r

(damping force)

2 3m

Note that both forces are function of the response parameter A, defined by Eq. (3.30). For

a S.D.O.F.-system, the dynamic equilibrium equation can now be given as

(3.37) (3.38)

11

/ill' + (G +

_ p

D CD _8A)

f+ kr = F(t) +

2p

D CD LA 2 37c 37c uw

Eq. (3.39) has to be solved by iteration as A is not explicitly known at the start of the

calculation.

For a system with small response amplitudes compared to water particle amplitudes, the hydrodynamic damping term can be omitted. Eq. (3.40) will then change:

m? + ci- + kr = F(t) +

D CD 8 U2

a cos cot

2 37c

(3.40)

(50)

This is easily seen as A becomes identical to the amplitude of the water particle velocity according to Eq. (3.30) and uw is defined by Eq. (3.26).

3.6 Linearization, stochastic waves

When dealing with offshore structures one must in many cases apply stochastic techniques for response analysis. This is due to the stochastic nature of ocean waves. This section does not include a complete description of how irregular sea is described as it is assumed that this is known from other sources.

Briefly summarized: Sea surface elevation is assumed to be a Gaussian process, and as a consequence of that, the velocity and acceleration at one specific point in the fluid will also be Gaussian processes. The question is then: How can non-linear drag forces be linearized in this case. This problem must be solved in order to apply frequency domain methods on offshore structures with significant drag force loading.

Basic assumption:

1) The water particle velocity at a specific point in the fluid (A) is a Gaussian process with mean value 0 and variance at , denoted N(0, ).

uw

2) For a linear system, the response process (displacement r) will also be Gaussian, and so will the response velocityi-. Both will have zero

2 2

(51)

3) The difference between two Gaussian processes will also be a Gaussian process. Hence,

the relative velocity vR will be a Gaussian process, also this process will have zero

mean.

The variance of the relative velocity will be givenby

2

G2 =

62 +-2

- 2cri, ai.p

vR u.

where p is the correlation coefficient between the two processes

GU CY?

Here, COV (uw, t) is the covariance for the two processes.

The difference between the non-linear and linearized force processes can now be written as

e = (qL - (NI)

1

= p D CD [KL -1 VR OVR

2 (3.42) (3.44)

In Eq. (3.44) both qL, qN-L and vR must be looked upon as stochastic processes, all Gaussian with zero mean.

The criterion for determining KL is no that the expected value of the error should be as small as possible, formally

E[e2] 0 (3.45)

afq

From these assumptions it can be shown that the linearization coefficient will appear as

1) =

COV(uw, fr) _(3.43)

r

(52)

KL =

where a" is the standard deviation of the relative velocity.

As for the harmonic wave case, the linearization coefficient depends on the response, which means that dynamic equilibrium must be found by an iterative procedure. During this

iterationCY,_R - the standard deviation of the relative velocity - has to be found.

As for the harmonic case, the drag force will be split up into one damping force following i-and an external force controlled by uw.

Stochastic linearization of multi-degrees-of freedom systems involves a more complicated formulation as correlation between loads on different degrees of freedom has to be considered. For more details, see refs. /1/ and /2/.

3.7 References

/1/ Leira, B.J. and Olufsen, A.:

"Biplanar linearization of drag forces with application to riser analysis". Offshore Technology Conferences, Houston 1986.

')/

Krolikowski, P. and Gay, T.A.:

"An improved linearization technique for frequency domain riser analysis". Offshore Technology Conference, Houston 1980.

(3.46)

(53)

4 METHODS FOR DYNAMIC

ANALYSIS

4.1

Introduction

In this section methods for solving the dynamic equilibrium equation when using the finite element method will be outlined. Comments on applications relevant to slender marine structures will be given.

The equation for dynamic equilibrium reads:

Mr + Cr + Kr = Q(t)

(4.1)

Definitions:

M: mass matrix, might be concentrated (then a diagonal matrix) or consistent, (then of

same character as K)

C: Damping matrix, might represent drag damping, or other damping sources. Can be

proportional to K and M due to practical

(mathematical) reasons, or specially

established from physical considerations.

K: The global stiffness matrix containing the elastic and geometric stiffness terms

(54)

Methods for Dynamic Analysis 4.2

r, t, I':

Vectors with displacements, velocities and accelerations.

The equation represents dynamic equilibrium for all degrees of freedom represented in the

system. One could say that inertia, damping and restoring forces are balancing external

forces. The equation is a generalization of the

static equilibrium case where restoring

(stiffness) forces are balancing external, time-invariant forces. Details on how these matrices are established will not be discussed, but some specific points relevant to slender structures will be made.

The mass matrix M contains terms both for the structural mass and hydrodynamic mass. Structural mass (inertia) is direction independent, but as the added mass term is not (see sect. 3), this means that M becomes a function of r, displacements. This is a kind of non-linearity that can be included in non-linear time domain solution algorithms only.

Often the mass matrix for offshore structures are frequency dependant as added mass for large-volume bodies is connected to wave generation. This frequency dependence is not relevant for slender bodies, which means that this type of complexity does not apply to the

systems of concern.

Note that the load vector also has contributions from forced displacements on the structure. Such displacements are typically found at the

end of an anchor chain or a flexible riser

connected to a floating vessel. When using tensioned risers on deep floaters, position control might be necessary at several levels, which means that forced (prescribed) displacements are found at several degrees of freedom. How the load vector terms are calculated from the displacements will not be detailed here. Some important facts should however be noted.

Forced displacements means that velocities and accelerations also will be prescribed. This means that the load vector will have contributions not only form displacements and stiffness as for the static case, but also from velocities/damping and acceleration/mass. These terms will appear for non-diagonal damping and mass matrices only.

Forced displacements can be correctly handled in all methods for analysis discussed in the following sections.

(55)

4.2 Eigenvalue analysis

For an undamped system without external loads, restoring and inertia forces are balancing

MF+Kr= 0

(4.2)

This is referred to as "free oscillation" and to solve this problem one has to solve the

eigenvalue equation

[K - 0)2i" MI 4 = 0

(4.3)

This equation is obtained by assuming that the displacements can be written as

r =

i sin cu,t (4.4)

Algorithms for solving the eigenvalue equation (4.3) will not be discussed. The solution will,

2 2 2

however, appear as a set of eigenvalues c02, o.)3 6)2] where n is number of degrees

of freedom in the system. To each eigenvalue coi there is associated an eigenvector Hence, n eigenvectors are found. The physical interpretation of these results are as follows: Each eigenvalue represents an eigenfrequency co or natural period T

27t

T = (4.5)

which is a possible free oscillation frequency. The eigenvectors represents the shape of the oscillation connected to the actual eigenfrequency. Hence, for a riser problem, the solution will typically be as illustrated on Fig. 4.1.

Note that the eigenvector Oi represents a shape only, not the numerical values of amplitudes. The amplitude for free oscillation is undetermined as all amplitudes of a modeshape will represent a possible equilibrium between restoring and inertia forces.

In the RISANA computer program one can calculate the eigenfrequencies and modeshape for a tensioned riser and have a plot of the modeshapes.

An important feature of the eigenvectors is that they represent a complete set of functions for the n-dimensional space represented by the degrees of freedom. This means that any possible

(56)

45k-t,

114-Ei T.

-100

lo-cLA-C_I-cie-ca`VA

Figure 4.1 Solution of eigenvalue problem, tensioned riser

state of deformation can be described as a weightened sum of the modeshape:

r - Ea

(4.6)

i=1

Note that this is not an approximation but a mathematically correct wayof describing any deformation shape. Approximations will, however, often be made by using a reduced number of modeshapes in the representation:

171

r

= E

ai 171 < (4.7)

1=1

This kind of approximation is commonly used, and the accuracy obtained depends on number of modeshapes that are included. This must be seen in relation to problem parameters

such as load frequency and load pattern. In a dynamic problem, the displacements are time dependant. Eq. (4.7) will then be written

as

r(t)

= E $Y

(t) -y

(4.8)

i=1

The matrix now contains the selected eigenvectors while the y vector is the time-dependent weights. The matrix multiplication is illustrated on Fig. 4.2.

a

\i

(57)

Figure 4.2 Illustration of matrix multiplication

Another feature of the eigenvectors is that they are orthogonally with respect to each other, and also with respect to the stiffness and mass matrices. This means that the following

relations are valid

iT

= 0 for i

j

7 '

ci) M tit

= 0 for i

j

(4.9)

K

= 0 for i

j

For all the above relations one have that the products are different from zero for i =j.

4.3 Mode superposition

The basic idea of the mode superposition method is to describe the displacements according

to Eq. (4.8).

By introducing this equation to Eq.

(4.1) (dynamic equilibrium) and premultiplying the result by the transposed eigenvector matrix (131. we have

1My+TC457 §T

K §

=

TO

(4.10)

The matrix multiplication is illustrated on Fig. 4.3 for a system using m out of n eigenvectors.

7a7-'

z

4)4 (1),i

(ti

(Pit

(58)

Figure 4.3 Illustration of matrix multiplication

Remembering the orthogonal properties of

the eigenvector we have that the M matrix

resulting from the multiplication

4.7-Mti.=

(4.11)

must be a diagonal matrix M 17711, mm]

(4.12)

The same will be the case for the K-matrix as the eigenvectors are orthogonal also to

K. If

the damping matrix C is of the proportional or Rayleigh type:

C =

M + a2 K

(4.13)

the eigenveetors will be orthogonal also toC. Hence, the originally coupled equation (4.1) is no transformed to a set of m uncoupled equations, each reading:

mi + yi

y1 = Q(t)

(4.14)

CD

:

where Q1(t) is the result of the vector multiplication

T T (4.15)

Q1(t) = Q(t)

often referred to as a generalized load.

The equations (4.14) can now be solved individually by a method depending on the nature of the load vector, and the final response will then be found by linear superposition:

0

r11-1_ ,c

®

ycr

r---c_,

. m2, +

(59)

r(t) = E 4. y i(t)

(4.16) ,i=1

In the equation above, the modeshapes are of course constant for in time while y(t) are found from Eq. (4.14).

For cases where 41 is not orthogonal to we have to solve the coupled equation of the type

11119 CS'

ie;

= ( t)

(4.17)

This equation have m free parameters instead, of n present in Eq. 'Vt. 0

The, method is efficient for linear systems Where, number of eigenveetors required to obtain an acceptable Solution is significantly smaller than number of degrees of freedom. Note however that stresses can only be found from the complete displacement vector (4,16), which means that this operation must be performed prior -to stress, calculation. For drag loaded structures the equations- (4,14) will be 'coupled through the load terms as the loads will be dependent on the total relative velocity square. The solution must therefore follow the coupled Scheme, and one must always refer back to the ,complete representation of displacements and velocities by use of Eq. (4.16) for load caleulatibn in a time domain procedure. This makes the method less attractive for riser application.

Note that the y parameters in Eq. (4.17) does na represent displacements, but so-called "generalized 'displacements", which in faa are scaling factors for modeshapes, defining each

modeshape's representation in the final 'displacement shape (Eq. (4.16)).

Number of modeshapes necessary to give good displacements can be significantly lower than what is, needed for stress analysis. This is due to the fact that stresses (at least bending stresses) are functions of the displacement second derivatives, which are significantly more influenced from higher order modestiapes as is the Case for displacements.

IS

(60)

4.4 The Frequency Response Method

The frequency response method solves the equation of dynamic equilibrium in the frequency domain, and the method can easily be used on problems with many degrees of freedom subjected to stochastic loads. The method is therefore quite popular among offshore structure analysists, and have been successfully applied both on fixed structures and risers.

The load is one arbitrary degree of freedom is described as

Qp) = Q

ei(6)t c9

(4.18)

= Qcv (cosai +

isin apei't

The real part of this complex representation is assumed to be the actual load.

Here:

Qoj load amplitude, d.o.f. j frequency of load

c:

phase angle of load

Note that the phase angle ai is not the same for all degrees of freedom, which means that loads at the same frequency can have different phases from one place (d.o.f) to another. The load vector will now be complex and a function of the load

frequency:

Total load vector: Xeia)

Element j in load vector: Xi = Q0i (cosaj + i sin ad

(4.19)

XRi + X11

(61)

MF+Ce+Kr=Xeic1)1

The particular solution of this equation can be written

r = x

ehot

where x is a complex vector

= roj (cosOi + i sin Of)

The response amplitude roi is given by

2 2

r01 = _\IxRj + x .

-and the phase angle Elj is given by 0 arctg(xulxRj)

The complex response vector x is found from Eq. (4.20).

[K- co2M + icoCix = X

or

Bx=X

(4.20) (4.21) (4.22) (4.23) (4.24) (4.25) (4.26)

where the B matrix is complex and a function of the frequency co (Eq. 4.25) Eq. (4.26) can also be written as

x(co) = H(co) - X(co) (4.27)

where

H=B-1

(4.28)

or

BH= 1

(4.29)

From, Eq. (4.26) and (4.29) we can see that the H matrix can be found by solving Eq. (4.26) for a system "loaded" by the unity matrix 1. The physical interpretation of H is therefore that H1j(03) is the response in degree of freedom i caused by a unit load at degree of freedom j. The load is harmonic and have the frequency co. See Fig. 4.4.

Note that the equation system (4.29) has to be solved for each frequency of interest. For

systems that are loaded by irregular waves, 50 to 80 frequencies might be needed. Number of load vectors will, however, not necessarily be equal to the total number of degrees of

(62)

LAA ZEST' c. E Lc ki>

t

`t"s

I 0 0

Figure 4.4 Illustration of the H matrix

freedom in the system. This follows from the simple fact that all dof's are not subjected to external loads.

The main advantage of the frequency response method is the elegant way the phase angle problem is treated. Once the H(co) is established, calculation of response due to other load patterns or intensities is straightforward and fast on the computer. Most of the computation effort is connected to the H(co) calculation (solving the equation system), specially if many frequencies are needed. The H matrices will therefore normally be stored in the computer (on disk) when using program systems based on this method. The methods works without calculation of eigenfrequencies, but such information is often needed in order to understand

the results. This method is well suited for stochastic analysis, often combined with fatigue life calculation.

The method is linear, and special considerations are needed when dealing with non-linear drag forces (stochastic linearization). Procedures for this type of analysis willnot be outlined here. Some comments are found in Sect. 5.

Transient response due to impulsive loading cannot be performed, as the particular solution due to harmonic load only is considered. (Steady state solution.)

1

(63)

Mass, damping and stiffness that are dependant on frequency can easily be handled in this type of analysis.

4.5 Time domain procedures

Time domain procedures means here to solve the dynamic equilibrium equation (4.1) by a time integration technique often referred to as "step-by-step time integration". Other time domain methods are also available like convolution or Duhamel integration, often referred to

as the impulse response method". Methods based on Hamilton's principle is also possible

but of marginal practical interest. Time integration can be performed on the coupled set of equations (like (4.1)) or on the individual equations derived by the mode superposition method

(Eq.(4.14)). Two method different approaches exist use of finite differences or numerical differentiation use of numerical integration

Of type a), Houbolt's method is the most well-known.

This method has, however, a

considerable (or at least significant) damping, which will appear as a decay of the oscillation when calculating response time history of a free vibrating system without damping.

The methods based on time integration are the most popular ones in structural analysis and will be further described here. Initially, the numerical procedure itself will be outlined, and after that algorithms for integration of the coupled equation (4.1) will be given.

S CA-0.0A

The basic principle is to in time and calculate dynamic equilibrium at given points in time. The time increment is problem dependent, but for most slender marine structure applications increments between 0.1 - 0.5 seconds are seen. Equilibrium between the discrete points is not ensured. Once the equilibrium condition in one time increment (k) is found in terms of displacements, velocities and accelerations, the condition in the next time increment (k+l)can be found.

01)

(64)

All integration methods are based on the same basic equations. These are two equations arriving from integration of an assumed acceleration time history in the actual time increment, and a third equation taken from dynamic equilibrium

k+1 1:11c+1 +

f ii(t)dt

(4.30) tk tk+ 1 uk+i = u k + f it(t)dt (4.31) tk u(t)

=1[Q(t) - c 140 - k

Lt(t)] (4.32)

Figure 4.5 Illustration of time integration

The different integration methods differs in the assumption of acceleration time history in the time interval [tk, tk+1].

The principle of time integration will be illustrated by going into some

details on the

procedure based on constant average acceleration. The basic assumption is shown on Fig. 4.6

c s.

4

(65)

-Methods for Dynamic Analysis 4.13

GISA

C4l-TrO

₎?

Wit,

(-rt.-UpL

(+)

cot-t, cA Lon, C. LA.C=eCrOr.r...k

, a(t)

_(uk

2 ,

it(t) = a

+ ku _.k+i)t

4, cp

CO

By integration as shown in Eq. (4.30) and (4.31) we have::

Atf

I a) 1:1 JO =

-2kuk Uk+111

b) Uk+ :=Lik + Atilk + _(1:4 _'+1))

1 01

lik+.1 = --k k +1 = C it,k+1 IkUk+i)

The equations (435) a, b and c contain three unknown parameters, uk+i and iik+i, and

can therefore be solved. The most convenient way of doing this, is to rewrite theequations:

4

,

Uk+1 --ikUk+1 - At-, 4

uk) -

At k Uk ) C ) (2.,. (4), t....14-t2 (4.35) (4.36) (4.33) (4.34) = + At2 + k + ..

(66)

2

fik+1tkuk+i

uk) Uk

4

uk+i = ±[Qk+i + m + c) uk + (_4m + c) ak + mak]

A

_At'

At At

4

where A = m +

2c

+ k

At2 At

These explicit equations gives an unconditional stable method that are frequently used in

structural dynamics.

When discussing methods for time integration, three keywords are important

stability damping

period error

Stability is directly connected to the procedure for numerical integration and not to system parameters like damping or stiffness. Loss of stability makes the solution diverge from the correct path as indicated on Fig. 4.7.

P.A.brot4t,

Figure 4.7 Loss of stability

Integration methods can be divided into two categories with respect to stability:

conditionally stable methods unconditionally stable methods

(4.37) (4.38) (4.39)

Response modelling of Marine Risers and Pipelines

OF

MARINE RISERS AND PIPELINES

ACKNOWLEDGEMENT

itaiesce_A.

TABLE OF CONTENTS:

INTRODUCTION

branching, and the axial force - often referred to as the

I

C.

2

FINITE ELEMENT FORMULATION

Introduction

a = C

C :

f T

U = - E adV

E C EdV

arbitrary position in an element is

given as function of node

u = N v

N :

I I

= E

(2.,)

Ti:

2.2 Two-dimensional beam element

A

1()

2 ax

,

= =()

2a

= --(---)

2 ax

NT

u =

yG dV

N -

(z./ g)

U

44._

2, 8

$

kg cs

K=

7

00

--0

T.

-L

...

=T K .

2.3

3-dimensional beam element

i i, two more coordinate systems are used. As shown on the figure, these are

'?=

L - L0

b

node b

2.4 Consistent vs. concentrated loads

c

(t--.7")

4.,111,

J

M = EI

P =

JNT

(2.31)

?

r,

P

N Tr) dX

-(2.R)

t-t

ti,,,

rl

-tct L,g,

K r = R

Si = k cv

_{= =()}

_{2 ax}

_C-c-rkso,

_{= p(Cm - 1)}

_{pCD DvR}