Finite-element analysis of the marine riser

Finite-el.ement marine riser analysis of the M. H. Patel and S. Sarohia London Centre for Marine Technology, Department of Mechanical Engineering, University College London, London, UK K. F. Ng CanOcean Engineering Ltd, Calgary, Canada (Received February 1983) A two-dimensional finite-element computational method is presented for determining marine riser displacements and stresses due to self-weight, buoyancy, internal and external pressures, surface vessel motions and environmental forces arising from currents and waves. A dynamic analysis is performed in the frequency domain for regular waves by linearizing the hydrodynamic damping term. A fully nonlinear time step integration of the equations of motion using the Newmark constant acceleration method is also presented. Extensive use is made of a substructuring technique in the dynamic analysis to reduce the number of degrees of freedom and, therefore, achieve a substantial reduction in computer time and storage without a discernible performance penalty. A number of test cases are presented to compare the linearized and nonlinear dynamic analyses and to highlight the performance of this method in relation to other marine riser analysis techniques. The application of this finite-element method to mu Ititube production risers with complex cross-sections is discussed. Key words: offshore engineering, marine risers, finite elements, wave dynamics Introduction A marine riser is a compliant vertical single or multiple pipe connection between the surface vessel and subsea wellhead on offshore drilling or production facilities. The surface vessel may be fully floating or of the tethered buoyant nature, where vertical motions are restrained by vertical tethers with substantial compliancy remaining in the horizontal surge, sway and yaw motions. The structural, dynamic and hydrodynamic behaviour of marine risers has been extensively researched. The emphasis has fallen on large fmite-element based com- putations of riser response to currents and waves. These computations are mainly formulated and solved in two dimensions, t-s although some three-dimensional methods also exist. 6'7 There are two problem areas that remain. The first is associated with the almost total absence of reliable proto- type or experimental data on riser motions and stresses in waves. This feature is generally accepted as a serious drawback to the credibility of the large number of fmite- element computational methods that are available. Reference 8 describes some experiments which will partially fulfd this need for data. The second problem area arises in the long computer times and large core storage required by the current generation of finite-element riser computations. The use of linearization techniques in frequency-domain dynamic analysis does lead to a decrease in computation time, but at the cost of a significant approximation in the calculated motions and stresses. On the other hand, a time-series computation models the non- linear behaviour more accurately, but is much more expensive to implement. The requirement to perform a large number of optimization computations in riser design pushes up the costs and poses a significant limitation. This is particularly true for production risers where the com- putation of fatigue damage in irregular waves requires long time-domain computer runs. A finite-element analysis of the marine riser in current and waves is presented in this paper. It is implemented with a linearized frequency-domain solution and a time-domain integration scheme. Emphasis is placed on reducing corn- 0141-0296/84/03175-10/$03.00 © 1984 Butterworth & Co. (Publishers) Ltd Engng Struct., 1984, Vol. 6, July 175 FE analysis of the marine riser: 114. 14. Patel et aL puter time and storage by the use of substructuring techniques in the dynamic analysis which remove just under two-thirds of the degrees of freedom from the calculation without a discernible loss of accuracy. A num- ber of test cases are presented to compare the linearized and nonlinear dynamic analyses and to illustrate the performance of the technique in relation to other finite- element analyses. Static analysis The governing equation for marine riser analysis has been derived by Young and Fowler. 9 Their detailed formulation includes the effects of both external hydrostatic pressure and internal fluid pressure in the riser. The governing equation is presented as: E1 --(..r + PoAo--PiAi) fiX ~ ~ dx 2 A -- [WR + WE +PigAi --PogAo]O~Z : q (1) with the reference axes as shown in Figure I (nomenclature is defined at the end of the paper). The coefficients A and B include the effects of external and internal fluid pressure around the riser pipe. The modified physical tension in coefficient A is presented as an 'effective' tension which must be taken to apply on the riser cross-section. There are a number of special features associated with the internal forces and structural loading of marine risers which have been extensively disctassed recently in the literature. A brief review is presented here. It is important to distinguish between the tension and non-tension contributing riser internal contents. Since the marine riser is a long, slender structure with relatively small bending stiffness, it needs to be kept in tension to prevent buckling collapse. Thus a tension is applied to the riser at its top; and it is the weight in air of the riser pipe, associ- ated choke and kill lines and a hydrostatic modifying term which cause a variation in tension along the riser's length. The weight of the separately tensioned drill pipe and the riser fluid contents do not directly affect the tension variation except through the fluid densities in coefficient A of equation (1). However, the non-tension contributing elements in a riser cross-section must be accounted for when computing an effective lateral force component (coefficient B in equation (1)). The derivation presented in reference 9 includes the effects of riser self-weight and buoyancy on the marine riser governing equation. These effects have also been described in reference 10 from a more physical point of view. For deep-water risers, the top tension requirement to prevent buckling collapse can become excessive. In order to reduce top tension, buoyancy modules can be attached along the length of the riser. The distribution of buoyancy modules influences the tension variation in the riser; thus altering its structural response and internal stresses. However, the increase in diameter of a riser cross- section due to buoyancy modules also increases current and wave forces. This introduces considerable scope for opti- mizing the intensity and distribution of buoyancy modules in deep-water applications. The finite-element analysis presented here is based on a governing equation of the form given by equation (1). The h n / E lement number t~-I i i i / I t I d~~ i' cl~ef x ' c 3 - -y ///////I//////// - '~6 - - - -~4 6 degrees of freedom 2 b ~, =1 Figure 1 Element and global node description: (a) riser element nodes; (b) single beam element riser pipe is idealized as an assembly of beam elements as shown in Figure 1. Each element possesses six degrees of freedom, two translations and one rotation at each end. Consequently, the numerical computation is two-dimen- sional with all external forces on the riser, including forces due to current and waves, acting in one plane. The current loading q per unit length along the riser due to a lateral drag force is: q = ~PoCDdlUclU c (2) where Po is the density of sea water, Co is a drag coeffi- cient, d is an effective riser external diameter, and U e is the vector of current velocity. The variation of current velocity with depth needs to be specified so that Ue is known for any point along the riser. The static analysis has also been used to relate riser deflections and stresses to current and wave loadings in a quasi-static manner. For a known regular wave height and period, the current velocity U c is superimposed on to the horizontal component of the wave particle velocity U w at any instant. The quasi-static hydrodynamic loading can then be written as: q = ½ Po CDdlUe + Uwl (u~ + Uw) + ¼ poTrd2Cm ~]w (3) Ow is the horizontal local component of wave particle acceleration and Cm is an applicable inertia coefficient. The total stiffness matrix Sm for each beam element is derived as the sum of the standard elastic stiffness matrix S e and a geometric stiffness matrix Sg, which is a function of deflected element geometry and the axial force on the element. Thus: Sm= So + Sg (4) 176 Engng Struct., 1984, Vol. 6, July For an element of length L and an 'effective' axial tension T', Sm in member axes is given by the sum of: E/ AL 2]I symmetric 0 12 0 6L 4/, 2 --AL 2/I 0 0 AL 2[1 0 --12 --6L 0 0 6L 2L 2 0 12 --6L 4L 2 (5) "0 0 6/5 T' 0 L/IO s'=T o o 0 --6/5 .0 L/IO 2L2/15 0 0 --L] 10 0 6/5 --L2/30 0 --L/IO symmetric 2L2/15 (6) where A i~ the area of steel cross-section, E is Young's modulus and ! is an appropriate second moment of area. The local effective axial tension T' is calculated by accounting for the modification due to hydrostatic pressures in the surrounding fluid, as described earlier. The fixed end action vectors Amr are obtained by using an assumed shape function N(x) in conjunction with a total lateral load distribution w(x). This load is due to both the hydrodynamic loading q and an effective lateral load derived from term B in governing equation (1) with dy]dx obtained from an initially assumed undeflected riser configuration. Thus: l AmL = - - /w(x) N(x) dx (7) o where x is the vertical distance from the bottom ball joint and l is the total riser length. The final static member and actions A m are then obtained from: A m=Am L +SmD m (8) where D m is the nodal displacement matrix. These com- bined end actions are applied incrementally in order to account for the changes in term B and the nonlinear behaviour caused by large deflections of the riser pipe. ThusA m is divided into a specified number of equal increments AA m which are applied progressively to obtain the incremental displacements AD through the equation: AD = S-IAAm (9) where D and S are the overall displacement and overall stiffness matrices in global coordinates. The overall stiff- ness matrix is re-evaluated after each load increment to account for the change in geometry due to large deflections. The static analysis is executed in different ways depending on the analysis which is to follow. For the frequency-domain dynamic analysis, the effects of current are considered in the preceding static analysis. However, the nonlinear time-domain method requires that the steady current, unsteady wave velocities and riser velocity be summed before applying the nonlinear square-law drag FE analysis of the marine riser: M. H. Pate/et al. force. For this reason, the static analysis preceding the time-domain dynamic analysis only accounts for self- weight, buoyancy and pressure forces on the riser and excludes the current velocity. The loading due to current is then accounted for in the time-domain dynamic analysis. Dynamic analysis Basic equations The differential equation of motion for a system with many degrees of freedom and having a mass matrix M T can be written as: MTff) + CD + SD= F (10) where D is the matrix of nodal displacements, and C and S are the structural damping matrix and the overall stiffness matrix respectively; all are defined in global riser axes. The external force matrix F due to wave action on the system is obtained from a modified form of Morison's equation: F = poV~]+ poV(Cm --1) (U-- /)) +BIU--DI (U--D) (11) where V is the vector of elemental volumes, B is the matrix of hydrodynamic drag coefficients, and U and 0 are the horizontal components of wave particle velocities and accelerations. It is assumed here that the fluid induced forces on a structure are given by the linear superposition of a drag force and an inertia force. The first two terms of equation (11) signify the Froude-Krylov and the added- mass force effects respectively, while the last term describes the drag force. By substituting equation (11) into (10) and replacing [34 + po(Cm --1) V] by M T, poCmV by M H and rearranging, we get: MTD+CD+SD=MHO+BIU--151(U--L) ) (12) The above matrix equation cannot be used directly for incorporating the boundary condition at the surface vessel, which requires that the riser top end must follow the surge motion of the surface platform. This known horizontal riser nodal translation at the surface (denoted by suffix B) can be separated from all the other unknown degrees of freedom (denoted by suffix A), through the following matrix partitioning: LSBA SBBJ MH BA MH BBJ UB [BBA BBBJ × [ iu,, + 0 [j Here, FB is a force required to cause the specified surge motion at the surface. The dynamic response of the riser structure in terms of the remaining degrees of freedom can be obtained solely from the upper set of equations from (13), which do not contain F a. Engng Struct., 1984, Vol. 6, July 177 FE analysis of the marine riser: M. t4. Patel et al. Element property formulation In the formulation of the beam element mass matrix, the lumped mass or the consistent mass approach may be used. In the former, the entire mass is assumed to be concentrated at nodes where the translational degrees of freedom are defined. For such a system, the mass matrix has a diagonal form. Off-diagonal terms disappear since the acceleration of any nodal point mass would only produce an inertia force at that point. The consistent mass formulation, however, makes use of the finite-element concept and requires that the mass matrix be computed from the same shape functions that are used in deriving the stiffness matrix. Coupling due to off-diagonal terms exists, and rotational as well as translational degrees of freedom need to be considered. In theory, this consistent mass approach can lead to greater accuracy, although this improvement is believed to be small. On the other hand, the lumped mass formulation is easier to apply because fewer degrees of freedom are involved, leading to a simpler definition of element properties. The lumped mass formulation is chosen for this analysis because the advantages of a small improvement in accuracy for the consistent mass approach are outweighed by the additional computational effort entailed in its implementation. It having been noted that off-diagonal terms OfMH, M x and B are zero for the lumped mass formulation, the following equations are obtained from (13): MT AADA + CAAD A + SAAD A = MH AA UA + BAA I(UA --DA)I (UA --DA) -- CABOB ISABD B (14) At the end of the static analysis the stiffness matrix of the structure in its deformed position is available. In modelling the dynamic response about this mean statically deflected shape, the stiffness matrix is assumed to remain constant throughout the dynamic analysis. In the lumped mass approach, all the rotational degrees of freedom need to be substructured out. Since vertical wave forces are not significant for the riser system, the vertical translation degrees of freedom can also be elimi- nated. This is the feature which leads to a substantial reduction in computer time and storage in the dynamic analysis. The horizontal degrees of freedom having been segregated, the force-deflection equations can be written in partitioned form as: SNH SNNJ where subscripts H and N denote the horizontal and the other group of vertical and rotational degrees of freedom respectively. From equation (15), DN = --SN~ SNH DH (16) The condensed stiffness matrix suitable for use in the equations of motion is then: SI~IH = SHH -- SHN S~l~q SNH (17) The matrix is further partitioned to separate out the top surge degree of freedom: SAA SAB (18) S~H= SBA SBB where subscript B denotes the vessel motion as before. The mass matrix for each element is buih up by co)l- centrating half of the total mass of mud, pipes and buoy- ancy material at each end of the elenlen(. For a fully submerged vertical element, lhe adde¢~ )nass associated with unit horizontal body acceleration is Po(Cm -- 1) V, where C m is an inertia coefficient which includes the Froude-Krylov effect. Taking half the added mass to be lumped at each node, tile added mass sub- matrix for each element is: ~po(Cm ~l) v (19) This added mass matrix and the real mass matrix are summed together to give the total mass matrix M T AA- The manner in which the partially submerged element at the water surface is idealized depends on the amount by which the element is wetted at the mean sea level. If the wetted length L s is less than half the element length, all the added mass is lumped at the lower node and the element submatrix becomes: where A x is the total cross-sectional area of the riser element, including buoyancy elements when present. Should L s be greater than half the element length L, the added mass associated with the lower half of the element is concentrated at the lower node, while the rest of the hydrodynamic effects are taken to act on the top node. The element submatrix for such a situation is: ½Po(Cm--1)AxL 0 0 P°(Cm--1)Ax(Ls--L/2) 1 (21) For the riser structure, this appears to be a simple and logical way to treat the element at the water surface in the lumped mass formulation. The hydrodynamic mass matrix MH AA, which includes the Froude-Krylov forces, is built up from element submatrices in a similar manner. The submatrices corresponding to equations (20) and (21) respectively are: 00] 0 poCmA×Ls 0 0 PoCmAx(LS --L/2)] (22) (23) Due to the unit relative horizontal velocity (U--/)), the horizontal drag force on a fully submerged element is ½ Po CDLd. The hydrodynamic damping submatrix for such an element is: 0 1_[¼p°CDLdo ¼poCoLd] (24) The corresponding submatrices for a partially immersed 0] ' forLs L/2 ~poCD(Ls --L/2) d ' (26) element are: [ ½ PO CoDLS d 178 Engng Struct., 1984, Vol. 6, July The structural damping matrix may be explicitly defined as C = otoM T + atS (27) To obtain the coefficients ao and al, the damping ratios ~1 and ~'2, in any two modes need to be specified. An eigen- value analysis is carried out to find the natural frequencies corre.sponding to the two modes chosen. For Rayleigh damping: 11 ~1 - - O31 0t0 o, From equation (28): 2(h 092 - ~'2 09 0 So = (29) 092 091 091 092 2(h/092 - ~2/091) al = (30) 092 O91 091 092 A damping ratio of 5% in the first two modes has been chosen for all the analyses carried out in this work. The actual level of structural damping that should be specified is rather unclear in the current literature. Solution in the frequency domain A linearized form of the equation of motion may be obtained by replacing the drag term in equation (14) with a suitable equivalent linear damping term which is pro- portional to the relative velocity (U w --DA). For such a linear system: MT AA/)A + (CAA +Beq AA)/)A + SAADA =M H AA(J+BeqA.AU--CABI)B--SABDB (31) Since the current velocity imposed is not sinusoidal, only the wave particle velocity U w and the structure velocity/)A can be included in the fluid interaction term. The stiffness matrix in the frequency analysis will therefore be obtained from the final statically deformed shape caused by current and riser internal forces. From linear wave theory, the elevation ~ of a single wave train may be represented by: = r cos (ky -- 09t) (32) The corresponding horizontal wave particle velocities and accelerations are given by: cosh (x- - l + h) k Uw = 09r cos(ky -- 090 (33) sinhkh 092r cosh(x -- l + h) k Uw - sin(ky -- cot) (34) sinhkh Rewriting equation (33) in complex form gives: = Re [a~, coshk(x - - I + h).ei~e_itot] Vw (35) t - sinhkh J or Re(Uw e -tt°t') (36) FE analysis of the marine riser:. M. H. Pate/et al. where 0 w is a complex amplitude. Similarly: /)w = Re(--i09Uw e -it°t) (37) In the steady state, the response of the system repre- sented by equation (3 1) to a sinusoidal wave will also be proportional to e -iwt. Thus: D A = Re(/) A e -iwt) (38) where/5 A is also complex. Differentiating equation (38) and substituting equations (36) and (37) into (3 1) gives: [SAA --MT AA 092 --i09(CAA +geq AA)]/)A = M H AA(--i09Uw) +Beq AAUw =/Y (39) where ff is a complex forcing function and Beq is an equi- valent damping matrix described in the following section. Since the matrix Beq contains a term in D, available only from the final solution, an iterative calculation scheme needs to be derived. Starting from a trial solution for the velocities/), Beq is estimated and the simultaneous complex equations (39) are solved for a new set of displacements and velocities DA and/)A: These velocities are compared with the original values (D) and the whole calculation is repeated with a better estimate of Beq until the real and imaginary parts of/) and/)A differ by a small specified tolerance. Linearization of the damping term. Since damping forces are responsible for the dissipation of energy in a vibratory system, the obvious, and the most common, way of obtaining Beq is to equate the work done by the linearized and the nonlinear forces such that: Beq(U-- D ) -=- BI(U--/9)1 (U-- / ) ) (40) For the purpose of illustration, a convenient node where y is assumed to be zero is chosen. From equation (33), the wave particle velocity is: U = R cos09t (41) where 09r coshk(x -- l +h) R - sinhkh Let the corresponding riser nodal velocity be defined by: /) = Q cos(09t - 4) (42) where Q is the amplitude of vibration and 4 is an arbitrary phase difference. The relative velocity is: (U--D) =R cos09t-Q cos(09t-4) = R T cos(09t -- g,) where R T = (R 2 -- 2RQ cos4 + Q2)1/2 (43) -Q sin4 tanxP - R -- Q cos4 The work done by the damping force BI (U-D) I (U - f ) ) over an elemental displacement dD may be written as: dW= BIRT cos(09t- @)IRT cos(09t--~) x Q cos(09t-4) d(09t) (44) On substituting/3 for (09t-- ~0), we can express the work Engng Struct., 1984, Vol. 6, July 179 FE analysis of the marine riser: M. H. Pate1 et al. done over a complete wave period by this nonlinear term as: 277-q w= s BIR, cosflIRT co@ -* xQcos(fl+?EI-$)d/3 Zn-* = s BIRT co@1 RT co@ Q cos(p + y) dp -* = f QBR; cosy (45) by splitting up the limits of integration to account for the modulus sign, and assuming that y = (\k-$) is tune independent. The work done by an equivalent linearized damping force B,,(U-D) over a wave cycle is readily obtained from: 2n W = I B,, RT cos(wt - \k) Q cos(wt - $) d(wt) 0 = aQBeq RT cosr (46) Finally, equating the work done by the two damping terms gives: B,, = $BR~ Hence, the damping coefficient can be used in equation (39): B,, =$I(~-&,I (47) Initial trial solution. To ensure that equations (39) con- verge rapidly to the final solution, a reasonably accurate initial estimate of the displacements D (and thus the velocities D) is required for evaluating the equivalent damping matrix from the total mass matrix and the diagonal terms of the stiffness matrix. Then, assuming a damping ratio of lo%, the initial estimate of B,, is taken to be: Be, = 0 I I 0 0.2 (M&&)r’2 I 0 --__ r-,--- nn ) (MT AASAA) I ,,* This matrix is substituted into equations (39), which are subsequently solved for the initial trial solution. This method leads to rapid convergence with only two or three iterations required for forcing frequencies away from the structure resonant frequencies. Up to 10 iterations may be necessary at and around these resonant frequencies. Time series analysis Wave and wave force specification. The basic method of analysis here involves integrating equation (10) through discrete steps in time and accounting for the nonlinear drag loading without a linearization approximation. In the equation of motion (equation (12)) the general- ized fluid velocity V can be decomposed in to the static current velocity UC and a wave particle velocity ZJ, . Thus equation (12) becomes: MT~+C~+SD=MHU,+BI(U,+DT,-~)~ x (U, + u, -6, (49) where oC is taken to be zero for the current velocity. The requirement to sum the current and wave velocities before applying the resultant loading through the square-law relationship requires that the current velocity be ignored in the static analysis that precedes this time-domain calculation. The time step integration of the equation of motion also allows irregular wave sequences (and the corresponding surface vessel surge responses) to generate the dynamic excitation forces on the riser. This wave sequence can be specified in two ways. A wave elevation spectrum of the incident irregular wave can be used to compute the corre- sponding spectra of the subsurface wave velocities and acceleration as well as the spectrum of surface vessel surge motions. These spectra can be Fourier transformed to generate corresponding time series of these quantities for use in the dynamic analysis. However, this procedure is cumbersome and computationally time-consuming. There- fore, a simple alternative method is used in this analysis. The incident wave elevation is specified as a ‘frequency comb’ sum of individual sinusoidal components with randomly distributed phase angles. The subsurface wave kinematics and surface vessel surge response can then be readily computed by summing the effects of all the sinu- soidal components in the wave spectrum. Newmark time integration scheme. The numerical time integration technique used here was initially proposed by Newmarkr2 and can be considered as an extension of the linearly varying acceleration approximation over a given time interval. The following relations can be used: d . . . . . t+At = 4 + [(l--6) D, + 6D,+A,l At D t+Ar = D, +&At + [(f -0’) ii, + fl’tijt+&] At2 (50) where /3’ and S are parameters which can be varied to achieve an acceptable integration accuracy and stability. Subscript t denotes the variable at the beginning of the time interval At. It is generally accepted that the direct integration analysis relies on an appropriate time step selection. The time step must be small enough to obtain sufficient accuracy, although a time step smaller than necessary would reflect on the cost of the solution. Bathe and Wilson13 analysed the stability and accuracy of various numerical integration schemes and suggested that, for reasonable accuracy, the time step-to-period ratio be not more than l/6 for the highest significant mode. In its standard form, the Newmark fl technique is unconditionally stable. The two parameters 6 and 0’ introduced in equation (SO) indicate how the acceleration is modelled over the time interval. 6 = l/2 and 0’ = l/6 correspond to a linearly varying acceleration. Newmark’s original scheme which is pursued here uses S = l/2 and 0’ = l/4 and gives a con- stant-average-acceleration based integration scheme. Using 180 Engng Struct., 1984, Vol. 6, July these latter values in equation (50) and rearranging gives: 4 fit+At = ~ [o t+~,t - ot - (A t ) btl - f i t (51) 2 z)t+~t = ~ [LJt+~t - f i t ] - ~t Then expressing equation (49) explicitly at instant (t + At) and using the lumped-mass approach with the top vessel surge motion duly separated as in equation (14), we get: MT AAI~A t+at + CAADA t+At + SAADA t+At = MH AAUwA t+at +BAA I(UwA + UCA --L)A)It+At x (UwA + UCA-- DA)t+,,t -- CABDB t+At -- SABDB t+At (52) Substituting equations (51) into (52) and rearranging gives: 2 = MH ~ grWA t+At + BAA I (UwA + UCA --/)A)It+at x (UwA + UCA --DA)t+~t -- CABDB t+~t 4 --SABDBt+At+[-~t MTAA+CAA]DAt 4 = Ft+at (53) This is the basic equation used in the time step integration scheme. The solution scheme assumes that the displacement, velocity, and acceleration vectors at time zero, denoted by subscript 0, are known and the solution is required from time zero to time r. The given time span r is subdivided into equal time intervals At (where At = r divided by the number of time intervals). The algorithm calculates the solution at the next required time from known information at the previous time steps. The process is repeated until the solution at all discrete time points is known. To initialize the numerical solution, the acceleration corresponding to zero time is derived from the reduced form of equation (53), giving: 6A = M;kA [M. ~OWA0 + B~ I(UwA + VCA)ol X (UwA -I- UCA)o -- CAB/)B. --SABDBo ] (54) In arriving at equation (54), the unknown value of velocity/)A t+At of the forcing vector of eouation (53) has been approximated to D A t. The approximation gives an acceptable degree of accuracy provided the time step chosen is sufficiently small. An alternative approach to this would require an elaborate iterative scheme with a significantly greater computation effort. From the set of simultaneous equations (53), the dis- placements are simply obtained from: DA t+At = J-lEt+At where 4 2 J = MT CAA + SAA At 2 AA +~'~ (55) (56) FE analysis of the marine riser: M. H. Pate/et aL The inversion of matrix J in the above equation can be made more efficient by the use of banded equation solvers as suggested by Bathe and Wilson. 13 However, J is inde- pendent of time and needs to be inverted once only. When D a t+At is known, the accelerations and velocities at (t + A t) are derived from equations (51). Results and discussion Computation validation procedures The finite-element calculations presented here have been validated by a number of methods. For the static analysis, the finite-element formulation has been checked by comparison with the analytic result for an idealized weightless tensioned beam./4 A typical such comparison is shown in Table 1. These comparisons confirmed the validity of the computational procedure as well as indicating the number of f'mite elements required for an acceptable level of accuracy. The American Petroleum Institute Committee on the Standardization of Offshore Structures defined a set of test risers as a basis for comparing the performance of riser analysis methods for both static and dynamic loadings. Nine anonymous participants to this study submitted solutions for the various test cases and API Bulletin 2J (reference 15) gives the overall comparisons. Table 2 gives the input data for one of the API test cases, and Table 3 displays the corresponding static analyses results. These are displayed in terms of maximum bending stress value and position, maximum total stress (axial plus peak bending), as well as upper and lower riser angles from the vertical. Results of the analysis presented here are given in Table 3 with the mean values from the nine API test cases displayed in parentheses. It is clear from these comparisons that the API program runs and i.he present method agree reasonably weU. The frequency-domain and time-domain dynamic analyses presented here have also been compared with the dynamic analyses in the API bulletin. Figures 2 and 3 show typical results for one of the AH test risers; the plotted API values are the maximum and minimum of the com- bined results from the nine calculations compiled in the bulletin. The frequency-domain analysis is computed conventionally using a regular wave period of 9 s and wave Table 1 Ten-element idealization of weightless tensioned 500f t beam Parameters: Total length 152.4 m Applied tension 54.422 tf Uniform load intensity* 6.036 x 10 -3 t/m Results: Analytical FE solution idealization Slope at ends (red) 0.007 1 0.007 1 Maximum moment (if-m) 0.863 1 0.868 9 Lateral displacements (m) Node 1, 11 0.0 0.0 2 ,10 0 .1044 0 .1044 3 ,9 0.191 4 0.191 3 4, 8 0.254 9 0,254 8 5, 7 0.293 3 0,293 1 6 0.306 1 0.306 0 * Equivalent to load caused by 0.Sm/s current Engng Struct., 1984, Vol. 6, July 181 FE analysis of the marine riser: M. H. Patel et al. Table 2 Input parameters for API cases Distance from mean sea level to riser support ring Distance from sea f loor to bottom ball joint Water depth Riser pipe outer diameter Riser pipe wall thickness Choke line outer diameter Choke line wall thickness Kill line outer diameter Kill line wall thickness Buoyancy material outer (tiameter Modulus of elasticity of riser pipe Density of sea water Density of mud Drag coefficient Added mass coefficient Effective diameter for wave/current load Density of buoyancy material Current at surface Surface vessel static offset Weight per unit length of riser joint in air Wave height Wave period Vessel surge amplitude Vessel surge phase angle 15.24 m 9.144 m 152.4 m 0.406 4 m 0.015 87 m 0.101 6 m 0.016 51 m 0.101 6 m 0.016 51 m 0.609 6 m 2.1 x 107t/m 2 1.025 t/m 3 1.438 t/m 3 0.7 1.5 0.660 4 m 0.160 2 t/m 3 0.257 4 m/s 4.572 m 0.256 5 tf/m 6.096 m 9s 0.609 6 m 15 ° Table 3 Results of static analyses of API cases. Mean figures from the API Bulletin 2J are quoted below each solution in parentheses Case Max. bending stress Max. total Angles from Value Location* stress vertical (KSI) (ft) (KSI) (deg) Lower Top BJ 500-0-1 3.99 104 5.46 3.64 0.44 (2.53) (111 ) (4.34) (2.94) (0.82) 500-0-2 1.62 104 6.98 2.58 0.96 (0.94) (116) (6.80) (2.20) (1.21) 500-20-1S 5.92 442 9.43 4.35 -1 .17 (5.86) (461) (9.51) (3.66) (--0.79) 500-20-2S 3.90 442 10.381 2.79 0.04 (4.27) (463) (10.54) (2.51) (0.24) * Above lower ball joint 200 150 E _ J ~ loo o c5 50 Figure 2 mum displacement; o, API results; zx, static displacement; quency domain; . . . . . , time domain J 1 I I I 1 2 3 4 5 6 X displacement (m) API case 500-20-1D. +, minimum displacement; x, maxi- , fre- 200 150 m _1 ~ loo . o cJ c ~5 5O 0 ~ -40 -30 -20 -10 0 10 20 Bending stresses,(N/mm 2) Figure 3 API case 500-20-1D. zx, static value; o, API results; ~, static + dynamic (min); +, static + dynamic (max); ~ , f requency domain; ... . . , time domain height of 6.096 m. The time-domain analysis uses a single- frequency 'comb' to produce equivalent data but with the nonlinear drag force due to current and wave velocities included in the calculations. It should be emphasized that none of the results published in the API bulletin has to our knowledge been directly validated by measurements on full-scale risers. Nevertheless this comparison gives an indication of agreement between the other methods and the analysis presented here. A further comparison is presented here between this analysis and an alternative full finite-difference calculation method 16 which maintains all six degrees of freedom per element throughout the static and dynamic analysis. Table 4 gives details of the riser test case with two features that differ from the test cases described in Table 2: these are that the riser has no bottom ball joint and is built in at the lower end, and further that there is no current profile loading the riser for the static case. Figures 4 and 5 show comparisons between the finite-element frequency- and time-domain calculations and the f'mite-difference methods displacements and combined (axial plus bending) stresses respectively. There is reasonable agreement between the three techniques for a 20 element idealization, confn-ming the validity of using a substructured approach with hori- zontal degrees of freedom only in the analysis presented here. Influence of nonlinearity on structural response A comparison of the time-domain and frequency-domain analyses presented in Figures 2-5 gives an indication of the effects of nonlinear fluid loading on the riser response. In the case of Figures 2 and 3, a static current prof'tle is included, and so time-domain and frequency-domain results differ markedly owing to the effect of the square- law drag force with and without linearization. However, the frequency-domain results are at lower values for the induced stresses. For the test case in Figures 4 and 5, the absence of a static current reduces the differences due to the drag force linearization, and the results for frequency- and time.domain calculations are closer together. However, a discrepancy still exists in the wave excitation zone where the loading due to wave velocities relative to the riser is treated by linearized and nonlinear techniques in the 182 Engng Struct. , 1984, Vo l . 6, Ju ly Table 4 Comparison case for DYNRISER Water depth Mean water level to riser support ring Wave height Wave period Tanker surge response Phase Top end Flex joint with 650 Ib ft/deg Bottom end: Encastre Riser pipe outside diameter Riser pipe inside diameter Connections: 200 Ib (0.091 t) every Hydrodynamic coefficients Top tension Offset Densities: seawater mud riser pipe material current velocity Modulus of elasticity E of riser pipe 302.6 ft (92.26 m) 2.624 ft (0.8 m) 27.5 ft (8.382 m) 8.3 s ±14.3 ft (±4.359 m) --103 ° 89.85 × 10 -S tm/deg 5.563 in (0.141 3 m) 4.063 in (0.103 2 m) 50 ft (15.24 m) C D = 1.2; C M = 2.0 27.7 tons (28.14 t) 0.0 m 1.025 t/m S 0.721 t/m S 7.850 tim ~ 0.0 m/s 2.1 x 107 tf/m 2 1OO 8O v ~" 60 .J o 0 # I 1 I I - ' 0 1 2 3 4 5 X. d i sp locement , (m) Figure 4 Comparison of UCLRISER and DYNRISER computer programs. 4, UCLRISER (frequency domain); v, UCLRISER (time domain); +, DYNRISER (time domain) IOC -x 6o g ~5 2O O I I I I ~ 0 50 100 160 200 250 300 Stress, ( N /mm 2 ) Figure 5 Comparison of UCLRISER and DYNRISER computer programs, zx, UCLRISER (frequency domain); v, UCLRISER (time domain); +, DYNRISER (time domain) FE analysis of the marine riser: M. H. Pate/et al. frequency- and time-domain calculations. Reference 19 describes additional results from these analyses, illustrating the kind of discrepancies due to nonlinear effects which are observed here. Discussion and conclusions The finite-element analysis and the frequency-domain and time-domain solutions outlined in this paper attempt to balance the small computing cost advantages of lineariza- tion against the additional accuracy available from the nonlinear time-domain calculation. The frequency-domain analysis uses the linearization approximation of equal energy dissipation between the nonlinear damping and the equivalent linear damping in the solution. An alternative linearization technique for frequency- domain analysis has been tested by Krolikowski and Gay (references 17 and 18). It is based on a statistical minimiza- tion of the mean squared error between the nonlinear damping force and its linear representation used in the analysis. The statistical approach seeks to use linearization at the discrete frequency components of a wave spectrum to arrive at a global linearized damping force with a least- squares minimized error. This technique allows a frequency- domain method to be applied over a wider frequency range, in contrast to the linearization method used in the analysis presented here which is used for regular waves only. The technique of linearization by least-squares minimiza- tion is not followed up in the frequency-domain analysis presented here. This is because both riser methods devel- oped here have been aimed at computing riser motions and stresses, the latter for feeding into fatigue calculations based on linear elastic theory or fracture mechanics. The fracture mechanics approach demands that representative stress time histories for a marine riser in waves be known in detail, particularly in terms of the sequences of stress cycles that are likely to occur. A computationally efficient time-domain analysis is capable of producing this informa- tion, whereas frequency-domain analyses, whatever their level of sophistication in linearization, operate in the frequency domain where the phase information which governs wave sequencing is lost. A further feature which has promoted the use of an efficient time-domain analysis for riser calculations is based on the comparative performance of the frequency-domain and time-domain analyses. A compendium of these results is presented in reference 19 and demonstrates substantial differences in peak stresses between the two analyses. These discrepancies may be reduced by a more sophisticated linearization technique in the frequency-domain analyses, but the discrepancies do highlight the importance of modelling the nonlinear fluid loading on the riser cross- section in a physically representative manner. An additional problem associated with marine risers occurs in the analysis of multi-tube production risers of complex cross-sectional geometries. These may be made up of a central structural riser with a number of large diameter satellite flow lines or as a bundle or array of flow lines. The beam finite-element analysis techniques described in this paper need to be extended to these production risers. Reference 20 suggests one solution by equivalencing a production riser of complex cross-section to a simpler single-tube marine riser, which is then used for the finite. element analysis. This approach is sufficient for a global riser analysis, but it needs to be used with care when localized riser fluid forces or member stresses are required. Engng Struct., 1984, Vo l . 6 , Ju ly 183 FE analysis of the marine riser: M. 14. Pate/et aL References 10 11 12 13 14 15 16 17 18 19 20 Harris, L. M. 'An introduction to deepwater floating drilling operations', Petroleum Publishing Co., 1972 Tucker, T. C. and Murtha, J. P. 'Non-deterministic analysis of a marine riser', Proc. Offshore Technology Conf. (OTC 1770), 1973 Burke, B. G. 'An analysis of marine risers for deep water', Proc. Offshore Technology Conf. (OTC 1771), 1973 Heuze, L. A. 'A 4000-foot riser', Proc. Offshore Technology Conf. (OTC2325), 1975 Larsen, C. M. 'Marine riser analysis', Norwegian Maritime Research no. 4, 1976 Gnome, E., Signorelli, P. and Giuliano, V. 'Three-dimensional static and dynamic analysis of deep-water sealines and risers', Proc. Offshore Technology Conf. (OTC 2326), 1976 Natvig, B. J. 'A large-angle-large-displacement analysis method for marine risers', presented at Intermaritec, Hamburg, Germany, September 1980 Hartnup, G. C., Airey, R., Wilkinson, P., Patel, M. H. and Sarohla, S. 'Hydrodynamic tests on marine risers', Proc. Offshore Technology Conf. (OTC 4319), 1982 Young, R. D. and Fowler, J. R. 'Mathematics of the marine riser', presented at the Energy Technology Conference and Exhibition, Houston, Texas, 5-9 November 1978 Morgan, G. W. and Peret, J. W. 'Applied mechanics of marine riser systems: Pts 1 to 13', Petroleum Engr, Oct. 1974 to May 1976, monthly parts Clough, R. W. and Penzien, J. ' Dynamics of structures', McGraw-Hill, 1975 Newmark, N. M. 'A method of computation for structural dynamics', J. Eng. Mech. Div. ASCE, 1959, 95 (EM3), 67 Bathe, K. J. and Wilson, E. L. 'Numerical methods in finite element analysis', Prentice Hall, 1976 Timoshenko, S. 'Strength of materials: Pt II - Advanced theory and problems', Van Nostrand, 1956 American Petroleum Institute 'Comparison of marine drilling riser analyses', API Bulletin 2J, 1st edn, January 1977 DYNRISER: Dynamic offshore marine riser analysis, SIA Computer Services, Ebury Gate, 23 Lower Belgrave Street, London SWIW 0MJ Krolikowski, L. P. 'Modern production risers: Pt 9 - The frequency domain method for production riser analysis', Petroleum Engr Int., 1981, July, 88 Krolikowski, L. P. and Gray, T. A. 'An improved lineaxisation technique for frequency domain riser analysis', Proc. Offshore Technology Conf. (OTC 3777), 1980 Patel, M. H. and Sarohia, S. 'The influence of nonlinear marine riser behaviour on methods of analysis and design', presented at the ASME Offshore Deep Seas Symp., New Orleans, March 1982 Patel, M. H. and Sarohla, S. 'On the dynamics of production risers', Proc. Third Int. Conf. on Behaviour of Offshore Struc- tures (BOSS 82), 1982, 1,599 Notat ion A Ai Ao Am AmL Ax B Beq C Co Cm d D b 15 Dm E cross-sectional area of steel in riser pipe inner cross-sectional area of riser pipe outer cross-sectional area of riser pipe vector of equivalent nodal actions vector of fixed end actions total overall cross-sectional area of riser element (including buoyancy element) matrix of hydrodynamic damping coefficients matrix of equivalent damping coefficients structural damping matrix hydrodynamic drag coefficient hydrodynamic coefficient (including Froude- Krylov force) outer diameter of riser element vector of nodal displacements in global axes vector of nodal velocities vector of nodal accelerations vector of nodal displacements in member axes modulus of elasticity F dynamic force loading vector g acceleration due to gravity h water depth I second moment of area of riser cross-section k wave number L riser element length Ls immersed length of riser element at free surface l total riser vertical length M H hydrodynamic mass matrix M T total mass matrix N matrix of shape functions Pi local hydrostatic pressure inside pipe Po external local hydrostatic pressure q lateral current load per unit length of riser QbR bRT velocity amplitudes in equations (41) and (42) r wave amplitude S overall stiffness matrix in global axes Se elastic beam stiffness matrix Sg geometric beam stiffness matrix Sm stiffness matrix in member axes t time T tension T' 'effective' local axial tension on riser element (= T + PoA o - -P iAi) U vector of nodal fluid velocities /J vector of nodal fluid accelerations Uc vector of nodal current velocities along riser length Uw vector of nodal horizontal wave particle velocities Uw vector of nodal horizontal wave particle accelerations V vector of elemental volumes W work done by damping force w total lateral load intensity on riser WE weight per unit length of independently strung drill pipe and contents inside riser w R weight per unit length of riser pipe x vertical distance from bottom end riser mounting y horizontal distance from vertical centre-line above sea bed connection ao, a~ structural-damping coefficients AA incremental load vector AD incremental displacement vector At time step fl a variable 3' a phase difference Pi density of riser contents Po density of sea water ~, 6o wave elevation, wave frequency w t, co2 natural frequencies of riser in two modes ~'1, ~'2 damping ratios in two specified modes phase differences 5, fl parameters in Newmark's relation r total time span in the dynamic analysis Subscripts and superscripts A all nodal horizontal translations apart from top surge B top surge degree of freedom H all horizontal translational degrees of freedom N all degrees of freedom besides horizontal translations overbar denotes complex amplitudes 0 values at time t = 0 184 Engng Struct., 1984, Vol. 6, July