S O M E A S P E C T S O F M A R I N E R I S E R A N A L Y S I S Mehernosh Boman Irani B.Ttch., Indian Institute of Technology, Bombay, 1980 M.A.Sc, University of British Columbia, 1982 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies Department of Mechanical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA May 1989 @ Mehernosh Boman Irani, 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Mehernosh Boman I r a n i Department of M e c h a n i c a l E n g i n e e r i n g The University of British Columbia Vancouver, Canada Date 11 May 1989 DE-6 (2/88) A B S T R A C T The problem of the static and dynamic analysis of marine risers is considered in this thesis. The equations of motion for a pipe conveying fluid are derived using a variational approach. The nonlinear expressions for the bending and torsional curvatures are obtained for a beam undergoing finite rotations. The elastic strain energy is ex-pressed in terms of these curvatures and Hamilton's principle used to derive the equations of motion valid for large deformations. The effect of internal flow is included in the variational formulation in a consistent manner. A finite element formulation is developed for the static analysis of compliant ris-ers. Large deformations are handled by using a convected coordinate system fixed to the elements. The nonlinear problem is solved using a combined incremental-iterative approach based on an instantaneous linearization of the Taylor series ex-pansion of the forces about the displaced configuration. The numerical results presented here show that the large displacement formulation can predict the static deflected shape of different types of compliant riser systems. The dynamic response of a typical marine riser system to excitations of the ocean waves, current and surface vessel motion is also studied. The hydrodynamic loading is represented by a general form of the Morison equation. The nonlinear drag force in the Morison equation is linearized using the method of equivalent linearization and a solution procedure in the frequency domain developed. Alternatively, the complete nonlinearity of the hydrodynamic loading is retained and the equations of motion are integrated in the time domain. The effect of this linearization on riser response predictions is evaluated. The method of equivalent linearization involves less computational effort and it approximates the response reasonably well. The effect of internal flow on the response of a typical production riser configuration is also studied. It is concluded that internal flow has a negligible effect on the dynamic response. C O N T E N T S ABSTRACT ii CONTENTS iii LIST OF TABLES v LIST OF FIGURES vi ACKNOWLEDGEMENT viii 1. INTRODUCTION 1 2. EQUATIONS OF MOTION 11 2.1 Equations of Motion of a Pipe Conveying Fluid 11 2.1.1 Kinematic relationships 13 2.1.2 Energy principle 19 2.2 Hydrodynamic Loading 25 2.3 Internal and External Fluid Pressures 29 3. STATIC ANALYSIS 33 3.1 Finite Element Formulation 33 3.1.1 Global degrees of freedom 34 3.1.2 Element degrees of freedom 41 3.1.3 Force-displacement relationship 43 3.1.4 Global system equilibrium 44 3.2 Numerical Solution Procedure 46 3.3 Numerical Examples 49 4. DYNAMIC ANALYSIS 63 • • • i l l 4.1 Equations of Motion 63 4.2 Solution in Frequency Domain 67 4.2.1 Linearization of the drag force 68 4.2.2 Numerical examples 70 4.3 Solution in Time Domain 81 4.3.1 Explicit time step integration 81 4.3.2 Numerical examples 82 4.4 Closure 86 5. CONCLUSIONS 91 REFERENCES 93 APPENDICES 98 I REVIEW OF LINEAR WAVE THEORY 98 II ELEMENT MATRICES FOR STATIC ANALYSIS 101 III ELEMENT MATRICES FOR DYNAMIC ANALYSIS 104 iv LIST O F T A B L E S 3.1 Data for "steep-S" riser configuration 54 3.2 Data for "steep-wave" riser configuration . 59 4.1 Data for API test case. , 71 V LIST O F FIGURES 1.1 A schematic of a drilling riser 2 1.2 Compliant catenary riser production system 4 1.3 Hybrid production riser concept 5 2.1 A schematic of riser configuration under study 12 2.2 A pipe conveying fluid 13 2.3 A schematic representation of the undeformed and deformed beam. 14 2.4 Eulerian angles xp,6 and <j> 16 2.5 Curvature components of the deformed beam 17 2.6 Equivalent representation of internal and external fluid pressures. 31 3.1 Elemental and global coordinate systems 35 3.2 Global degrees of freedom 38 3.3 Definition of incremental rotations 40 3.4 Element degrees of freedom 42 3.5 Cantilever beam subjected to moment at the free end 51 3.6 Cantilever beam subjected to torque and specified tip displacements. 52 3.7 A comparison of bending moments and torque 53 3.8 Typical "steep-S" riser configuration 55 3.9 Static solution strategy. 56 3.10 Static deflected shape as affected by the riser length 58 3.11 Typical "steep-wave" riser configuration 60 3.12 Static deflected shape of the "steep-wave" riser configuration. . . 61 4.1 Conventional twelve degree-of-freedom beam element . 65 4.2 Typical conventional riser configuration 72 v i 4.3 A comparison of the displacement with the API data 74 4.4 A comparison of the bending stress with the API data 75 4.5 Typical displacement results for a drilling riser in 453 m water depth. 76 4.6 Typical bending stress results for drilling riser in 453 m water depth. 77 4.7 Effect of internal flow on the riser response 79 4.8 Effect of structural damping on the riser response 80 4.9 A comparison between results obtained by the time stepping and equiv-alent linearization methods (H = 6.1 m) 84 4.10 A comparison between results obtained by the time stepping and equiv-alent linearization methods (H = 12 m) 85 4.11 Minimum-maximum displacement envelope for a compliant riser sub-jected to specified incident waves 87 4.12 Minimum-maximum displacement envelope for a compliant riser ex-posed to specified incident waves and current 88 4.13 Minimum-maximum displacement envelope of compliant riser due to specified incident waves and surface vessel motion 89 1.1 Definition sketch for linear wave theory 98 v i i A C K N O W L E G E M E N T The author wishes to thank his supervisor Dr. V.J. Modi for his interest and valuable suggestions. The author is indebted to John F. Baldwin, Dr. M. Isaacson, and Dr. N.D. Nathan for their help and encouragement. The author would also like to express his gratitude to the Hydraulics Laboratory of the National Research Council of Canada for its support and assistance. The financial support of the Natural Sciences and Engineering Research Coun-cil is gratefully acknowledged. • • • V l l l 1. INTRODUCTION Offshore technology has experienced a remarkable growth since the late 1940's, when drilling platforms were first used in the Gulf of Mexico. Today, fixed platforms have been installed in over 300 meters of water and drilling has extended to much greater depths. Presently, offshore structures are predominantly related to the petroleum industry but considerable developments in other fields such as ocean energy and mining have also taken place. As the search for oil, gas, minerals and alternate energy sources continues, offshore structures are being located in ever increasing depths and in harsher environmental conditions. This calls for greater efforts towards improvements in the existing techniques of design, analysis and construction. The marine riser is an important structural subsystem of any offshore platform used for oil and gas recovery. The riser in its simplest form is a pipe which connects the wellhead at the seabed to a fixed or floating platform or a vessel. It may be a production riser or a drilling riser. The riser must be able to sustain production or drilling operation while undergoing deflections due to platform motion, waves and currents. In the case of a simple drilling riser (Fig. 1.1), a steel pipe extends from the ball joint at the top of the blowout-preventer to the surface vessel. Ad-ditional components of the riser system consist of guide lines, choke-and-kill lines, and buoyancy and tensioning devices. When the riser forms part of the product transportation system it may comprise of several pipelines ('import', 'export', test and control lines). The system may consist of bundled or individual risers with separate tensioning devices. As production and exploration move to deeper waters, simple riser systems have been modified and alternative compliant production riser concepts have been pro-posed [l]. These concepts offer advantages in severe environmental conditions and 1 Drilling Platform 7 Tensioning System Lower Ball Joint , B l o w o u t Preventer S e a R o o r Well Head Figure 1.1 A schematic of a drilling riser. 3 deep water applications due to their ability to accommadate large motions. They also facilitate the development of smaller marginal reserves due to the relatively low installation and operation costs. The layout of one such compliant riser configuration is shown in Fig. 1.2. The use of flexible hoses (single or multitube) and subsurface buoys gives rise to a free hanging, catenary type profile. The flexible hose has a composite construction of interlocked steel and plastic layers [2]. Typically, the bending stiffness of a conventional steel pipe riser would be three orders of magnitude greater than that of a flexible pipe. Henceforth, the term 'conventional' riser would normally refer to a vertically aligned steel pipe system and a 'compliant' riser would mean a flexible hose system. Figure 1.3 shows a hybrid compliant riser concept developed by Mobil [3]. It consists of two parts: a lower rigid section, terminating about 100 meters below the sea surface, with a tensioning buoy assembly and a compliant upper section made up of a bundle of flexible hoses. Thus the lower section stays relatively undisturbed below the agitated surface zone and the upper section is compliant to waves, currents and surface vessel motions. Regardless of the type of riser configuration under consideration, marine riser analysis involves the study of a slender, flexible structure in the ocean environment. As such, this type of structural member forms major components in any offshore engineering application. Hence the problem of 'riser' analysis is equally applicable to flexible jacket platforms, individual tendons of a tension leg platform, mooring lines of a catenary platform, OTEC (Ocean Thermal Energy Conversion) cold water pipes, submarine pipelines, etc. The riser related fluid-structure interaction dynamics has received a great deal of attention in the past. Reviews of riser analysis include those by Chakrabarti and Frampton [4], and Narzul and Marion [5], while more general aspects of offshore 4 Figure 1.2 Compliant catenary riser production system. 5 Production Platform Compliant Section Rigid Section Sea Floor Figure 1.3 Hybrid production riser concept. 6 engineering are discussed in the texts by Sarpkaya and Isaacson [6] and Chakrabarti [7]. In what follows some of these aspects are touched upon. The main areas of interest in riser analysis are: (t) Mathematical modelling of the structure; (it) Evaluation of the hydrodynamic loading; (tit) Solution techniques. The first step in any riser analysis involves the mathematical modelling of the structure. This requires derivation of the equations of motion for a tensioned beam undergoing large three-dimensional deflections. In order to develop the equations of motion for large deformations of a beam, the general nonlinear theory of elas-ticity needs to be considered. The nonlinearities are both geometrical and physical in origin. Geometric nonlinearity is associated with the necessity to consider the deformed configuration in writing the equilibrium equations and the need to include nonlinear terms in the strain-displacement relations. Physical nonlinearity is due to the relations between the components of stress and strain being nonlinear. In riser applications the material remains elastic so that it is adequate to consider geometric nonlinearities only. The geometric nonlinearities can be considered up to varying degrees of ap-proximations, which in turn restrict the extent of rotations that can be adequately modelled. The majority of formulations for the analysis of 'conventional' risers are based on the assumptions usually employed in the Euler beam-column theory [8], where the rotations are considered neglible compared to unity and that the elonga-tions and shears are much smaller than the rotations (Patel et al. [9], Wang [10], Verbeek [11]). Riser problems involving large displacements can be handled by one of two pro-cedures. Firstly, the effects of finite rotations are included in deriving the equations of motion. Bernitsas [12] and Kim [13] use this approach and present equations 7 of motion for the large nonlinear three-dimensional deformations of risers. Unfor-tunately, the equations are not in a form amenable to numerical solution and the authors eventually resort to linearization. The second approach, which is more extensively used, is based on the finite element method. The structure is subdivided into a number of elements and the equations of motion are derived for an element based on the assumption of small rotations. The problem of large rotations is handled by the use of a moving co-ordinate system fixed to the elements. McNamara and Lane [14], Chucheepsakul [15] and McNamara et al. [16] present large displacement analyses of risers using the finite element method, but these are restricted to two-dimensional motion only. Chung et al. [17], and Felippa and Chung [18] show that in deep water applica-tions (typically beyond 3000 meters) of ocean mining pipes the deflected structure assumes a three-dimensional shape and large deflections produce highly nonlinear response. Recently, Martisen and Bergan [19], and Vogel and Natvig [20] have pre-sented general models for studying large three-dimensional motions of compliant risers. A production riser conveys oil or gas from the well-head to the surface ves-sel. It is known that internal flow affects the dynamic characteristics of a pipe and could even lead to buckling as well as oscillatory instabilities at sufficiently high flow velocities [21] - [23]. Paidoussis and Issid [24] have presented a comprehensive historical review of this subject and extensive experimental and theoretical inves-tigations of different aspects of the problem have been reported (Benjamin [25], [26], and Thurman and Mote [27]). Although the dynamics and stability of pipes conveying fluids have been studied for application in several fields, there is only one investigation the author is aware of which studies the effect of internal flow on marine riser dynamics. Wu [28] has investigated the effect of internal flow on the response of a conventional riser system. He develops the equations of motion for small deformations and solves them using a singular perturbation technique. This 8 restricts the application to specific configurations only. Hwang [29] has studied the effect of internal flow on the response of an OTEC cold water pipe using a modal superposition technique. A more general formulation to study the effect of internal flow on marine riser dynamics is needed. The second step in riser analysis involves the evaluation of the hydrodynamic loading on a slender structural member by assuming that flow separation dominates the fluid loading and that the incident wave kinematics does not change significantly in the vicinity of the structure. This aspect of the problem may be divided into two parts. The first part involves an adequate mathematical description of the water particle kinematics (i.e., velocities and accelerations) for a given sea state. The second part involves computation of the hydrodynamic loads from the known water particle kinematics. The problem of describing the water particle kinematics requires the use of an appropriate wave theory. Airy or linear wave theory is usually used, but it is necessary to consider a nonlinear wave theory for extreme waves. If a current is present, a simple superposition of the wave and current flow fields is generally used to calculate the total water particle kinematics. Ideally, wave-current interaction effects need to be taken into account [30]. Ismail [31] has shown that waves signifi-cantly modify the ambient current shear which affects the estimated hydrodynamic forces. Once the water particle kinematics are established, the in-line hydrodynamic forces are usually evaluated using a modified form of the Morison equation. This assumes that the in-line force on the structure is given as the sum of the inertia and drag forces. Due to the nature of the drag force, use of the Morison equation results in equations of motion which are nonlinear in both the response velocity and the fluid velocity. The nonlinearity can be handled by a numerical time domain analysis. Several linearization schemes may be used which result in simplified analyses. These are discussed in detail by Lipsett [32]. 9 Procedures used for the solution of the static or dynamic problem are numerous. However, they can be classified into two major categories: analytical and numerical. In general, analytical methods have proved to be successful for a riser which is uniform along its length [33]. A common procedure is to assume a mathematical model for deflection of the riser, satisfying geometric boundary conditions [34], [35]. Numerical methods for the spatial solution are more general in that they allow for changes in geometry. The two main approaches are: (i) direct solution of the partial differential equation by numerical approximations, and (ii) finite element methods. Examples of direct numerical methods are finite difference and numerical integration procedures [36]. The finite element method has found extensive use [10], especially for general three-dimensional, large deformation problems as discussed earlier. The solution of the dynamic problem with time dependent loading falls into two main categories: deterministic (time history or frequency domain) and non-deterministic (stochastic). A time history analysis is more general in its application and can handle nonlinearities without modification, but is computationally expen-sive. Frequency domain solutions are more suitable for fatigue analyses and involve shorter computer runs. However, the equations of motion have to be linearized. The stochastic approach involves determination of response statistics of interest for a given random sea-state. Generally, the equations are linearized so that the powerful spectral analysis methods can be used to obtain the response spectrum [37]. Sometimes repeated numerical simulations of the nonlinear equations are used to obtain statistics of the solution variables. The method, called Monte Carlo simulation [32], is computationally very expensive. The objective of this thesis is to develop tools that are useful for the design and analysis of marine risers. The development is also applicable to a wide range of structural systems. A variational approach is used to derive the equations of motion of a pipe conveying fluid and undergoing large three-dimensional motions (Chapter 2). The hydrodynamic loading due to the ocean waves and currents is 10 evaluated using the Morison equation and the effect of internal and external static fluid pressures is included in the formulation. Chapter 3 presents the finite element formulation for the large deflection anal-ysis of a framed structure. A combined incremental-iterative solution procedure for the nonlinear equations is discussed. Results are compared with some examples from the available literature and static deflected shapes of typical compliant risers are obtained. Dynamic response of the riser system is studied in Chapter 4. Two proce-dures are explored here. First, the drag force in the Morison equation is linearized and a frequency domain analysis performed. Secondly, the full nonlinearity of the hydrodynamic loading is retained and a solution procedure in the time domain is presented. Finally, Chapter 5 concludes this thesis by summarizing the present work. 2. EQUATIONS OF MOTION The study of a marine riser in general involves the analysis of a pipe conveying fluid in the ocean environment. Hence the first step in riser analysis would be the derivation of the equations of motion for a pipe conveying fluid. As discussed in the introductory chapter this aspect has been studied extensively. But a thorough literature review reveals that the formulations developed so far are restricted to the assumptions inherent in the Euler-Bernoulli beam theory, i.e., small deformations. Also, for a pipe conveying fluid there is mass transport across the boundaries which leads to nonconservative forces. This must be accounted for and care must be taken when using a variational approach to derive the equations of motion. Unfortunately, this point has been generally overlooked (e.g., Housner [23]) leading to improper boundary conditions. Here a general derivation is presented, for the first time, which clarifies the various aspects involved in the formulation of the equations of motion for a pipe conveying fluid and undergoing large deformations. To achieve this the mathemat-ical modelling of the system shown in Fig. 2.1 is considered. First, the equations of motion for large three-dimensional deflections of a pipe conveying fluid are derived using a variational approach. The hydrodynamic loading due to ocean waves and current is evaluated using a modified form of the Morison equation. Finally, the effect of internal and external static fluid pressures is included in the derivation using'the concept of effective tension and weight. 2.1 Equations of Motion of a Pipe Conveying Fluid In this section the equations of motion of an element of the riser oriented arbi-trarily in three-dimensional space (Fig. 2.1) are considered. To start with, expres-sions for the bending and torsional curvatures of the deformed beam are obtained in 11 12 Heave K±J Yaw Waves Element Coordinate System Sea Floor Figure 2.1 A schematic of riser configuration under study. 13 terms of the elastic displacements. The strain energy of the system is expressed in terms of these curvatures and the governing nonlinear, coupled equations of motion are derived using Hamilton's principle. The effect of internal flow is included in the derivation by a consistent energy approach. 2.1.1 Kinematic relationships The system consists of a pipe of length L containing a fluid flowing at a constant velocity Ui as shown in Fig. 2.2. y Ui (internal flow) Figure 2.2 A pipe conveying fluid. The governing equations are derived by considering this as a straight, slender, pris-matic beam which can undergo bending in two directions, torsion, and extension. A schematic representation of the undeformed and deformed beam geometries is shown in Fig. 2.3. All deformations of the beam are referred to the inertial coordinate system 14 Figure 2.3 A schematic representation of the undeformed and deformed beam. 15 x,y, z. The z axis is taken to be aligned with the centroidal axis of the undeformed beam. When the beam deforms, the centroidal axis at an arbitrary section translates by amounts u,v,tw in the x,y,z directions, respectively, and the section twists by an amount <f> about the centroidal axis. The goal is to arrive at expressions for the curvature components of the deformed beam in terms of the elastic deformations u,v,w and <j>. The coordinate axes xu, j/o, ZQ are fixed at an arbitrary point on the centroidal axis of the undeformed beam. Before deformation, x0,y0,z0 a r e parallel to x,y,z, respectively. Deformations u, v, w and <p displace the xo,yo, ZQ triad to x\, t/i, z\ and rotate £1,2/1,21 to 14,1/4,24, where the axis 24 is tangential to the deformed axis. The orientation of £4,1/4,24 with respect to £1,3/1,2! maybe expressed in terms of Eulerian angles. Figure 2.4 shows the Euler angles %p,6 and <p where the rotations are taken in the following order: a positive rotation rp about the yi axis resulting in £2,2/2522; a positive rotation 6 about the £2 axis resulting in £3,1/3,23; a positive rotation <f> about the 23 axis resulting in £4,1/4,24. Assume that the origin of the £4,2/4,24 coordinate system moves along the deformed centroidal axis of the beam with unit velocity such that the £4,1/4 and 24 axes are always directed along the principal torsion-flexure axes of the beam (Love [38]). Kirchoff's kinetic analogy states that the instantaneous angular velocity components of the £4,2/4,24 system with respect to the x,y,z system correspond to the components of curvature of the deformed axes [38]. Figure 2.5 illustrates this where the space derivative with respect to the curvilinear coordinate, s, along the deformed centroidal axis corresponds to the time derivative for the angular velocities. Thus the bending curvatures kx, ky and the torsional curvature k^ are obtained by projecting dxp/ds, 80/ds and d<j>/ds along the £4, t/4, and 24 axes and are given Figure 2.4 Eulerian angles xp, 6 and 4>. 17 2/1,2/2 2/3 ^3, *4 Figure 2.5 Curvature components of deformed beam. 18 as: 30 dtp kx = — cos <f> + —— cos 0 sin ^ ; O S O S K„ = —— cos & cos <p — —— sin © ; From the definition of Green's strain tensor [39], the derivative with respect to 5 can be related to the derivative with respect to z as, » - ( ! + * , - / • » . ...(2.2) where e is the nonlinear extensional strain of the centroidal axis, e = w' + -{u'2 + v'2 + w'2) , •••(2.3) and the prime denotes derivative with respect to z . Hence equation (2.1) can be rewritten as: kx = (l + 2e)~1/2[0'cos<£ +V>'cos0sin^] ; ky = {l + 2e)~1 / 2[^'cos0cos<£ - 0'sin<£] ; * * = (l + 2 £ ) - 1 / V - ^ ' t » h i 0 ] . •••(2.4) Also the following expressions for the trigonometric functions in terms of the displacements can be obtained [39]: 1 + w' costp (1 + 2e - v'*)1/2 ' (l + 2e-t;' 2) 1 / 2 (l + 2e) 1 / 2 19 sinV> = u [l + 2e - v f2) 1/2 ' sin 0 = —v (1 + 26) 1/2 (2.5) It should be noted that the angles of rotation are finite and transformations corresponding to finite angles of rotations are not commutative, so that different forms of the curvature expressions will result if the order of the Eulerian rotations is altered. The exact curvature expressions (Eq. 2.4), which are highly nonlinear, may be approximated to any desired degree by applying the binomial theorem. Retaining nonlinearities upto the second degree in displacements and their derivatives, the curvature expressions reduce to: kx = -v" + 2w'v" + v'w" + u"<f> ; ky = u" - u'w" - 2w'u" + v"<t> ; k4=<p'- w'<p' + u'V . • • • (2.6) 2.1.2 Energy principle A variational approach is used to derive the equations of motion of an elastic pipe containing a fluid flowing at a steady velocity. Note that for the problem on hand there is mass transport across the boundaries of the system. This leads to nonconservative forces interacting with the system so that care must be taken in applying Hamilton's principle to this type of problem [40]. The extended Hamilton's principle [41] equates the first variation of the Hamiltonian to the virtual work done by the external nonconservative forces as the system goes from a configuration at 20 time t\ to a configuration at time ti, S f h{Tu-V)dt=- I'2 {6WNF + 6Wne)dt, .-.(2.7) Jtx Jt1 where: Tu = kinetic energy of the system; V = potential energy of the system; 6Wnf = virtual work done by fluid forces at the boundaries; 6WNE = virtual work done by other external nonconservative forces. The potential energy of the pipe of length L, due to bending, torsion and extension can be expressed as, Li where : M x = EIxkx ; My = EIyky ; Mf = GJk.fi ; T = EAe; E = modulus of elasticity ; Ix^Iy — principal area moments of inertia about z and y axes, respectively ; A = cross-sectional area of pipe ; GJ = torsional rigidity of pipe . The kinetic energy of the system, Tu, consists of contributions from the kinetic energy of the pipe, Tp, and the kinetic energy of the fluid flow within the pipe, Tj, Tu = Tp-rTf. • • • (2.9) Neglecting the effects of rotary inertia and shear deformations the kinetic energy of the pipe can be written as, TP = J{im(«2 + v 2 + u,2) + \lrf}dz , • • • (2.10) 21 where m is the mass per unit length of the pipe; Ip, mass moment of inertia of the pipe about its axis; and the dot denotes derivative with respect to time. The kinetic energy of the fluid flow in the pipe is written as, Tf = j [\piAiUa-Ua]dz, -..(2.11) L where p,- is the mass density of the internal fluid; Ai, the internal cross sectional area of the pipe; and U_0, the absolute velocity of the fluid which can be expressed as, U a = (u + Uiu')x +{v + UiV')y + U i Z . • • • (2.12) Here Ui is the internal flow velocity and x, y and z are unit vectors in the x, y and z directions, respectively. The virtual work done by the fluid forces on the pipe at the boundaries can be written as, 6Wnf = [Fx6u + Fy6v + Fz6w] , • • • (2.13) z=0,L where F x , F y , and Fg are the fluid forces acting on the pipe in the x,y and z directions, respectively. These fluid forces act on the pipe due to the net rate of change of momentum through a section [42] and are given by the vector expression, £ = - / H a P i l L a - d A , ..-(2.14) A where the integral is evaluated over the cross-sectional area. Using equations (2.13) and (2.14), the extended Hamilton's principle (Eq. 2.7) is rewritten as, ; / (ru - v)dt 22 / 2 PiAiUi[{u + Uiu')6u + {i> + Uiv')6v + U{6w] dt Jt + /"*2 / piAiUi\{u + Uiu')6u + (t> + Uiv')8v + UiSw] dt •••(2.15) Note that the second and third integrals on the left hand side of this equation (Eq. 2.15) represent the net work that the fluid imparts to the pipe at the boundaries z = 0 and z = L, respectively. For pinned or fixed boundary conditions, the virtual displacements at the boundaries will be zero and these additional terms will vanish. The expansions for the kinetic and potential energies along with the nonlinear curvature expressions (Eq. 2.6) are substituted for the functional given by equation (2.15). Finally, taking the variation of the functional the following four Euler equations along with the appropriate boundary conditions are obtained: tt-coordinate (m + PiAi)u + 2piAiUiu' + PiAiUfu" + EIy{uiv - {u'w")" - 2{w'u")" + {v"4>) II + {u"w"Y - (u'w")' - 2{w'u")' + {v"<pw"Y - 2(u'V)" + 2 ( u W ) " + 4(u/V')" - 2{v"<f>w')"} + EIx{-(v"<p)" + 2{w'v"<p)" + {v'w"<j>)" + (uV)"} + GJ{{<p"v')" - {w'<p'v')" + (uV2)"} - EA[u'{w' + \{u' 2 + v'2 + w' 2)}}' = /, ; with boundary conditions: EIyu"'-EAu'{w' + \{u'2 + v'2 + w' 2)} = 0, or 6u = 0 at 2 = 0, L; and EIyu"' = 0, or <5u' = 0 at 2 = 0, L. v-coordinate (m + PiAi)v + 2PiAiUiv' + PiAiUfv" - EIx{-viv + 2{w'v")" + (t>V)" + («"<£)" + 2{v"w')" - 4{w'v"w')" - 2{v'w"w')" - 2{u"<pw') - {v"w")' + 2(w'v"w")' + {v'w" 2)' + {u"<pw")'} + EIy{{u"<p)" - {u'w"<f>)" - 2{w'u"<p)" + {v"<j>2)"} - GJ{{<p'u")' - {w'<p'u"Y + (u"V)'} -EA{v'{w' + hu' 2 + v'2 + w' 2)}}' = fy ; 24 with boundary conditions: EIxv'"-EAv'{w' + i(u' 2 + v' 2 + w' 2)} = 0, or 6v = 0 at 2 = 0,1; • • • (2.20) and EIxv"' = 0, or 6v' = 0 at 2 = 0,1. • • • (2.21) w-coordinate mw - EA{w' + i(u' 2 + v' 2 + u;'2)}.' + + \{u' 2 + t>'2 + w' 2)}]' - EIx{-2{v"2)' + 4{w'v" 2)' + 2{v'w"v")' + 2(u"2<f>)' - {v" 2)' + 2{w'v" 2)' + {v'w"v")' + (u'WY + («V')" - 2{w'v"v'Y' - ( t / V ) " - (tt>')"} - EIy{{u'u"Y' - (u'2w"Y' - 2f> W ) " + {y"<t>v'Y' - 2(u"2)' + 2(uW')' + 4(w'u"2Y - 2{v"<t>u"Y) + 2GJ{(<f>'2Y - (w'</>'2)' + (uW)'} = /, ; ... (2.22) with boundary conditions: PiAiU? - EA{w' + \(u' 2 + v' 2 + w' 2)} = 0, or 6w = 0 at z = 0; • • • (2.23) and p t^l7? + JE;>l{«;,+ i(u , 24-t; / 2 + Ti;'2)} = 0, or tfu; = 0 at z = L. • • • (2.24) ^-coordinate. 25 IP4> - GJ{4>" - {w'4>')' + (uV)' - {4>'w')' + (u/V)' - (u'Vu/)'} + EIx{-v" + 2w'v" + v'w" + u"(f>}u" + EIy{u" - u'w" - 2w'u" + v"4>}v" = U, ••• (2.25) with boundary conditions: GJ(f>' = 0, or ^ = 0 at z = 0, L. ••• (2.26) Here fx,fy,fz and are the components of external forces acting on the system. The mathematical model derived in this section is valid for large deflections where the rotations are less than unity. This restriction comes through the use of the binomial expansion for approximating the curvature expressions (section 2.1.1). The solution of the resulting nonlinear, coupled equations is not an easy task and no attempt is made to solve them in their present form. However, the equations do provide a perspective and some appreciation as to the complex character of the problem. It also helps in assessing the significance of the errors involved in the various simplified models which appear in the literature. Note, the proper boundary conditions are derived by including the effect of the nonconservative fluid forces on the system boundaries. Neglecting their contributions could lead to incorrect results in certain class of problems. 2.2 Hydrodynamic Loading The in-line wave force on a slender offshore structure is usually computed using the Morison equation. In their landmark paper, Morison et al. [43] assumed that the wave force on a structure is given as the sum of an inertia force, // and a drag force, fr). The original form of the Morison equation applies to a fixed structure 26 and for a typical circular section of diameter D, the force per unit length is / = / / + / D = PojDtCmU" + \ P o D C d U w \ U w \ , ••• (2.27) 4 2 where / v i s the fluid density; Uw, the undisturbed fluid velocity; and, Cd and Cm are the drag and inertia coefficients, respectively. The Morison equation is empirical and thus its application to complex time dependent flows can be questioned. However various attempts at more precise rep-resentation of the fluid forces have not generally been satisfactory. As such, given the uncertainties in modelling the ocean environment itself, the Morison equation has proved reliable in predicting the wave forces on slender members and is exten-sively used in the design of offshore structures. It is obvious that appropriate values of Cd and C m need to be used for accu-rately predicting the wave forces. For the fundamental case of a two-dimensional uniform sinusoidal flow past a circular section the force coefficients have been found to depend on the Reynolds number, the Keulegan-Carpenter number and the rela-tive roughness. The Reynolds number is the ratio of the inertia force to the viscous force and is written as Re = UmD/v, where u is the kinematic viscosity of the fluid and Um is the maximum flow velocity. The Keulegan-Carpenter number, K, is pro-portional to the ratio of the maximum fluid particle displacement, to the diameter of the cylinder and is given by K — UmT/D where T is the period of the flow. Numerous investigators have addressed this problem and a vast library of data on Cm and Cd are available from extensive laboratory and field tests. A comprehensive summary of the dependence of C m and Cd on these parameters is given by Sarpkaya and Isaacson [6]. The Morison equation in its original form applies to a fixed structure. If the structure is moving with a velocity Up which is colinear with Uw, a modified form 27 of the Morison equation can be used to predict the wave force [7], f = fi + fD = PojD 2CmUt> - p0-D 2[Cm - 1)U 4 4 w •••(2.28) The first two terms on the right hand side represent the inertia force on a fixed body in an accelerating inviscid fluid plus the contribution due to an accelerating body in an otherwise quiescent inviscid fluid [44]. The term p0irD 2(Cm — l)/4 is commonly called the added mass. The drag force is taken to be proportional to the square of the relative velocity between the fluid and the structure. In order to evaluate the hydrodynamic loading due to ocean waves and currents on a structural member arbitrarily inclined in space, consider the pipe element shown in Fig. 2.1. The fluid velocity relative to the element in the global XYZ coordinates can be expressed as: U r = U rJ + U ryj + Urzk; • • • (2.29) where: XT =U W + U C -U p; Uy = fluid particle velocity due to waves ; U_ c = current velocity ; ZJf = velocity of the pipe ; A A A *)i>^ = u1 1^ vectors in X, Y, Z directions. The components of the relative velocity in the x and y directions are given by the dot products: 28 U ry = Ur-y . •••(2.30) Where x, y are unit vectors in the x and y directions, respectively: A A A x. = xxi + yxj + zxk ; A A A y = xyi + yyj + Zyk . • • • (2.31) Similarly, the components of the fluid particle acceleration are: U? = U w-x; = U W -y . •••(2.32) The component of the hydrodynamic forces tangential to the pipe axis is ne-glected and the Morison equation is used to evaluate fx and fy, the components of the forces normal to the pipe axis, in the x and y directions, respectively: fx = ±PoDeCdUrxy/(Uy + U?) + Po-D\Cjj™ ; / „ = \p0DtCdUrvyJ{U? + U**) + Po^DlCmU? ; • • • (2.33) where De is the effective hydrodynamic diameter of the pipe section. Finally, the forces per unit length in the global X, Y, Z directions may be obtained as: fx = fx xx H" fy xy j fv = fxy* + fyyy; Sz = SxZx + Syzy. ••• (2.34) As before, it should be noted that there is no theoretical basis for this modified form of the Morison equation. The force coefficients are empirical and have to be 29 determined for a complete range of flow parameters. It is beyond the scope of the present study to consider the effect of variable force coefficients and, as is usually the case, representative values of Cm and are used. 2.3 Internal and External Fluid Pressures The riser pipe is subjected to external hydrostatic pressure from the surround-ing sea water and internal pressure due to the presence of the drilling mud or oil and gas. In this section it is shown how the effects of the internal and external static fluid pressures are included in the equations of motion. For conciseness, the nonlinear terms in equations (2.16), (2.19), (2.22) and (2.25) are neglected and the equations of motion are rewritten in the form: ,92u d 2u od 2 u (m + + 2 * ^ — + P^U 2^ d r„dui d 2 r d V -*[T*] + a ? [ E I ' a ? l m f - i ,d 2v d 2v ,d2v (m + P i A i ) — + 2 * A ^ — + HAtf^ ~il[ Ta^} +d^[ EIza^\- fv] d 2w dT t d 24> d 2d> where the axial tensile force T is now given by, T = EA[w' + \{u' 2 + v' 2)} . ••• (2.36) These equations are consistent with the level of approximations of Euler beam-column theory, namely, the rotations are negligible compared to unity and the 30 strains are much smaller than the rotations. The force due to the fluid pressures can be evaluated by integrating the pressure force along the wetted perimeter of the pipe section. This method is tedious and may not give the exact force for complicated sectional shapes. An alternative method, which is extensively used, includes the effect of external and internal fluid pressures as statically equivalent forces acting on the riser. The details of this approach have been very lucidly explained by Sparks [45]. The procedure is illustrated in Fig. 2.6 which shows the free body diagram of a pipe segment of length dL acted upon by internal and external fluid pressures and the tension on the pipe walls (Fig. 2.6a). This force and the pressure system acting on the pipe segment are decomposed into three parts (Fig. 2.6b). The external and internal pressure fields are separately closed by the addition of "missing pressures" as if the element were disconnected from the rest of the pipe and submerged by itself. The equivalent forces acting on the segment due to this "closed" pressure field is simply the buoyancy upthrust (p0gA0dL) on the segment and the weight of the contained fluid (pigAidL). In order not to modify the total force system it is necessary to add to the true wall tension (T), forces equal and opposite to those created artificially by the addition of the above mentioned "missing pressures". The sum of these forces, as shown in Fig. 2.6(c), is called the effective tension. Hence the influence of the external and internal fluid pressures, pa and p,-, is achieved by replacing the tension T by an effective tension, Te = T + PoA0 - p^ , ••• (2.37) where A0 and A,- are the external and internal sectional areas, respectively. Similarly the weight per unit length of the riser in air, W, is replaced by an effective weight, 31 T + dT T / (true wall tension) Pi (internal fluid pressure) ' p0 (external fluid pressure) WdL (true weight) p0gA0dL (upthrust) PigAidL (weight of internal fluid) T + dT (Po + dp0)Ac T - piAi + p0A0 + (dT - dpiAi + dp0A0) (W + PigAi - PogA0)dL PiAi + p0A0 Te + dTe WedL (effective weight of riser and contents) Te (effective tension) 2.6 Equivalent representation of internal and external fluid pressures. 32 WE = W - P0QAo + pigAi. ••• (2.38) For a riser generally the tension at the top is known. Therefore, a simple procedure would be to calculate the effective tension (Te) at any point from the top tension and the effective weight (WE) of the intervening length of the riser (plus contents). Thus the effects of internal and external fluid pressures are represented by these statically equivalent forces which produce identical effects as far as the bending of the beam-column is concerned. 3. STATIC ANALYSIS The first step in the analysis of compliant riser.systems involves a definition of the static equilibrium state of the riser due to the action of the 'dead loads' (i.e., self weight and static fluid pressures). This static equilibrium state is used as the reference configuration on which the various 'live loads' (i.e., ocean currents, waves and surface vessel motions) are imposed. The definition of the static equilibrium configuration is not a trivial problem as large displacements occur due to the ex-treme flexibility of the system. The strains may be small, but the relative rotations are large, i.e., line segments will not change significantly in length but they may translate and rotate appreciably. Hence it is necessary to determine the nonlin-ear solution for the static equilibrium configuration using incremental or iterative procedures. In this chapter a formulation for the analysis of the riser undergoing large deformations is presented. This formulation is based on the finite element method and was originally derived by Nathan [46] and Radomske [47] in the context of structural stability problems involving large deformations. It is extended here to the analysis of compliant riser systems. First, the finite element discretization procedure for the large displacement problem is presented. An incremental-iterative scheme for the solution of the nonlinear equations is described. Then some numerical results involving large deformations are presented. For conciseness, basic finite element procedures are not discussed and reference can be made to standard texts for details (e.g., Bathe [48], Zienkiewicz [49]). 3.1 Finite Element Formulation The nonlinear equations of motion derived in the previous chapter include 33 34 the effects of finite rotations within the element boundaries. Discretizing these equations would involve the use of sophisticated higher order elements to derive the element stiffness properties. Here a more efficient approach is used wherein the structure is subdivided into a number of elements so as to reduce the rotations and translations of the elements within acceptable bounds for the linearized form of the equations to be applicable. A less refined element stiffness is needed and the problem of large deformations is handled by the use of a moving element coordinate system fixed to the elements. To start with the structural deformations are defined with respect to a global inertial reference frame. Next, the local moving coordinate system associated with each element is defined and the element stiffness matrices are formulated. Finally, the system response in the global coordinates is deduced from the element equilibrium equations. 3.1.1 Global degrees of freedom Consider an element JK, connecting nodes j to k, arbitrarily oriented in space with respect to a global inertial reference frame X, Y, Z (Fig. 3.1). All geometry (nodal coordinates) and global forces are ultimately related to this global reference frame. Each element has twelve global degrees of freedom (f), three translations (r*) and three rotations (r*) at each node. The translations, E* = ( r i , r 2 , r 3 , r 7 , r 8 , r 9 )T , • • • (3.1) are measured in the global X, Y, Z axes direction. The rotations, E1" = {r4, ^ 5. ^65 rich ^ n , ri2) T , • • • (3.2) being finite are expressed in terms of Eulerian angles. At nodes j and k of element JK, there exist local reference frames whose coordinate axes are defined by the orientation of the member in space. The local 35 Figure 3.1 Elemental and global coordinate system. 36 reference frame x 3,y 3 ,z 3 (Fig 3.1) has its origin at node j. The x 3 axis is tangent to the centroidal axis of the element. The y 3 axis is coincident with the major principal axis of the element section while z 3 is coincident with the minor principal axis. Similar observations can be made for the local reference frame x k,y k,z k at node A;. Initially, before any deformations are applied to the structure, the element is straight and the x 3,y 3,z 3 and x k,y k,z k frames are aligned in space. The initial orientation of the element with respect to the global reference frame is defined by the finite angles ( r ° , rg, r|?)T. Hence the initial orientation of the local reference frame can be related to the global reference frame by an orthonormal transformation: (;) -4} -"m where transformation matrix xp is a function of the angles ( r ° , r^, r^) . As the structure (and the element JK) deforms, the local reference frames move with respect to the global reference frame. The current orientation of the local frames are again related to the global X, Y, Z frame by the orthonormal trans-formations: and, where tp 3 is a matrix function of the large global rotations (rj, r|j, r £ ) T which correspond to (r^, r^, TQ) and define the current (deformed) orientation of the element tangent at node j. Similarly tp is a matrix function of the large rotations (r10> rl H r12) a t n o d e k-In order to evaluate r*, the deformations need to be added to the angles defining 37 the initial orientation of the element. Hence at node j of element JK, fr°4 + r4\ = rg + r 8 . •••(3.6) The finite rotations (r^, r]j, rg) are taken in a specified order as shown in Fig. 3.2. To start with, the x\,yi,z\ axes are parallel to the global coordinates X, Y, Z respectively. First, the rotation r\ about the x\ axis results in x2,y2,Z2- Secondly, the rotation r£ about the y2 axis results in £3,1/3,23. Finally, the rotation rg about the 23 axis results in x }, y J, zK Corresponding to this rotational sequence the transformation matrix is: ( cos(r|) cos(rg) cos^) sin(rg) - cos^) sin(r£) cos(rg) ^\ + sin^J) sin(rg) cos(rg) + sin^) sin(rg) rp> = — cos(r£) sin(rg') cos^) cos(r )^ cos^) sin(r£) sin(rg') - sin(f4) sin(r|) sin(r|) + sin(r^) cos(rg) s i n(r5 ) - sin(rj) cos(rj) cos(r^) cos(r£) J •••(3.7) The tj)_ transformation composed of the large rotations (rJo, r^ , r^) at node k is similar in form to xp 3. The configuration of the structure is completely determined by the structural degrees of freedom r which correspond to the elemental global degrees of freedom r as shown in Fig. 3.2. The global degrees of freedom of an element can be related to the structural degrees of freedom as, ¥=br, • • • (3.8) where 6 involves the large Eulerian rotations r*. But r* defines the orientation of the element in the deformed state and therefore is not known a priori. This relationship is simplified for application to the incremental solution procedure. It is assumed 38 Figure 3.2 Global degrees of freedom. 39 that the structure deforms gradually as the loads are applied in incremental steps and the corresponding displacements and rotations are small. Thus equation (3.8) is written as an incremental relationship 6r = b_6r, (3.9) where the incremental rotations of 6f are taken to act about axes whose orientations are assumed not to change during the current incremental step, i.e., r* is known. Figure 3.3 shows the incremental rotations at node j from which the relationship between the incremental deformations 6f and the structure incremental degrees of freedom 6r follows. The b matrix can be written as where: *1 = b = o ki Q o o o / o Vo o o bzJ o cos(rj) sin(rj) cos(r5r) 0 sin(r*) cos(r4) co s ^ ) (3.10) 1 0 0 \ 0 cos(rJo) sin(r^ 0) .s i n(r* i ) -s i n(ri o )c o s(ri i ) c o s(ri o ) c o s(ru ) / and I = 0 = Figure 3.3 Definition of incremental rotations. 41 3.1.2 Element degrees of freedom There are six local or element degrees of freedom per element. The six degrees of freedom 5 are measured with respect to the local moving coordinate system x 3 ,y 3 ,z 3 as shown in Fig. 3.4. The 5 are kept small with respect to the length of the element by suitably increasing the number of elements required to model a given structure under a large rotation or translation. The local nodal deformation vector s with respect to the convected system is illustrated in Fig. 3.4 and is defined as s = (si, S2, s 3, s 4, s 5, s 6 )T. • • • (3.11) The translations (si, s2, S3) act along the local x3,y 3,z 3 axes. As the 5 are kept small with respect to the length of the element, the rotations (54, 55, 56) are small and are assumed to be about the same local x 3,y 3 ,z 3 axes. The local degrees of freedom s are functions of the global degrees of freedom r defined previously. Figure 3.4 shows the element JK before and after the finite displacements r, in which the element deformations are exaggerated for clarity. It is assumed that initially node j is at the origin of the X, Y, Z axes. The element has an initial length L with projections li,l2,l3 in the X,Y,Z directions, respectively. Due to the global displacements r* the element translates as a rigid body and also deforms as shown. Using appropriate transformations, the local degrees of freedom are related to the global displacements by the following expressions: si = (/1 + r 7 - r^U, 1) + (/2 + r 8 - r2)VJ'(l, 2) + {l3 + r 9 - r 3)V y(l, 3) - L ; s 2 = (/1 + r 7 - n)xp 3(2,1) + (/2 + r 8 - r2)^'(2,2) + (/3 + r 9 - r 3)0 y(2,3) ; s 3 = {h + r7~ ri)fp> (3,1) + {l2 + r 8 - r2)VJ(3,2) + (/8 + r 9 - r3)VJ'(3,3) ; 54 = -^(3,1)VJ'(2,1) - V*(3,2)VJ'(2,2) - Vfc(3,3)V'(2,3); z * node j (initial) Figure 3.4 Element degrees of freedom. 43 s6 = V*(3,1)^(1,1) + ^(3,2)^(1,2) + V*(3,3)V y (l,3) ; 5 6 = -V f c(2,1)^(1,1) - tpk(2,2)^(1,2) - V*(2,3)VJ'(2,3) . •••(3.12) 3.1.3 Force-displacement relationship The local force displacement relationship determines the element stiffness ma-trices. The element forces Si are first determined as nonlinear functions of the displacements s: Si = Si(s). Next, they are expanded in a Taylor series about some instantaneous element displacements s° By making the displacement increments s suitably small, second and higher order terms can be discarded, and the equations linearized as Si{s) - Si{s°) = Si -Sf = ASi l 5 • • • (3.13) where or A S = k As , •••(3.14) where k = dSjjs) dsj The matrix k can be written as k — kg -\- k_i , •(3.15) 44 where the stiffness matrix kg is independent of SQ and the matrix kx is linear in S Q . Use of these matrices is based on the assumption that the displacements must be small compared to the element dimensions. The stiffness matrices are presented in Appendix II. These were given by Nathan [46] but have been transformed to conform to the present element definition. 3.1.4 Global system equilibrium Now the system response in global coordinates is deduced from the element equilibrium equations. The element or local degrees of freedom s have been ex-pressed previously as functions of the coordinates f , Si = f{¥). -..(3.16) The element incremental displacements can be expressed as ds Ssi = ^6r} = a^(r)bjk(f)6rk , • • • (3.17) or in matrix notation, 6s = a*b6r = g_6r , •••(3.18) where the transformation matrix a is a function of the global coordinates f. By the principle of virtual work, the external work done by the R forces moving through the corresponding virtual displacement Br must be equal to the internal work of the element forces 5 moving through the compatible virtual displacements 6s. Thus, 6r TR = 6s TS_. •••(3.19) This compatibility condition (Eq. 3.19) gives, 45 If equation (3.20) is true for an arbitrary virtual displacement 5r, then, R = a S . (3.21) The can be expanded in a Taylor series about an initial global displacement r° on the assumption that there are no singularities, dRi(r) d 2Ri(r) AryAr f c (3.22) + If the increments r are restricted such that the second and subsequent terms of the expansion can be neglected, then or in matrix form, A j R i = V ^ i ( r ) A r , ^ dry V J r 0 1 ' J — AR = KAr , (3.23) (3.24) where K = . dr3 J r O The jth column of K_ is given by dR(r) dr, _ d_ ro ~ drj 46 r m = 9? [ko + ki 9L*J+ ^{¥jk+a*c4)Tbij '-t=i (3.25) where m is the total number of structural degrees of freedom, ay is a column vector of zeros with a one in the jth. place, and: c. = drdri db (3.26) Thus the global stiffness K_ is given by K. — Ko + K.i + K.2 > (3.27) where: and the columns of Kj = k_i(§°)a ; »=i Thus given the global displacements r, f and the element forces 5 at any stage, the global instantaneous or tangent stiffness of the structure K_ can be evaluated. This stiffness defines the tangent plane to the load-displacement surface at the instantaneous displacement configuration. 3.2 Numerical Solution Procedure The formulation for the large displacement problem involves, in principle, a set of nonlinear algebraic equations which can be expressed as K{r)r_=R, (3.28) 47 where: K = overall stiffness matrix ; r = vector of structural displacements; R = external force vector. A solution of these nonlinear equations (Eq. 3.24) can be obtained using one of three approaches: an incremental method; an iterative method; or some sort of combined incremental-iterative approach. In the incremental method of solution, the loads are applied in small increments and the structure is assumed to respond linearly during each load increment. The solution at a particular value of the load is taken as the solution at the end of the previous increment plus the solution to the locally linearized equations. Thus the set of linear algebraic equations solved at the tth increment can be expressed as KiAri = ARi, •••(3.29) where A r are the incremental displacements corresponding to load increments AR, and K_ the global instantaneous or tangent stiffness of the structure depends on the global displacements and is re-evaluated at every load increment. The advantage of this method is that the complete load-displacement history of the structure can be followed even for highly nonlinear problems. But in this method the errors that are generated within each increment accumulate making the solution progressively less accurate. The use of many small increments for acceptable accuracy may require extensive computational effort in some problems. Alternatively, there are several iteration methods which may be used for the solution of Eq. (3.24). The Newton-Raphson method has found extensive use in structural mechanics. This iterative scheme corresponds to the principle that if a structure is in an equilibrium configuration, the externally applied loads balance the internal forces. An imbalance between the external and internal loads would 48 produce a displacement increment and equilibrium will not prevail. Thus an equi-librium configuration can be obtained by the iterative process: K_j A r j = ARj ; r y+i = r y + A r y ; • • • (3.30) where K_j i S the tangent stiffness at the j'th iterate, and A r are the displacement increments produced by the load imbalance Ai2y = R e x t — Rint- Here R e x t are the external forces acting on the structure and R_int are the internal elemental forces transformed to the global coordinate system. The iteration is carried on till the external and internal forces sum to zero, i.e., when the structure attains equilibrium. The Newton-Raphson scheme is numerically stable and efficient and its main strength is that the accuracy of the solution at a particular value of the forcing function is not dependent on the accuracy obtained at any other value. But to guarantee convergence, an initial trial solution is required in the neighbourhood of the exact solution [50]. If the nonlinear equations are a slight perturbation of a set of equivalent linear equations, the initial solution can be taken as the solution to the equivalent linear set. Unfortunately, for certain type of structures the nonlinear equations are not slight perturbations of a linear problem. To overcome the drawbacks of the above mentioned methods, a combined scheme employing incremental load steps with subsequent iterations is used here. In this way the load-displacement path of a structure can be followed efficiently using large load steps for the incremental method. At the same time an accurate solution is assured by performing the Newton-Raphson iterations for convergence to equilibrium. This approach combines the advantages of both methods and can be applied to a wide class of problems. However certain types of flexible structures (e.g., cable type structures) do not possess any initial stiffness. This leads to unsta-ble equations when the first load increment is applied to the unstressed structure. 49 In such cases the above solution procedure cannot be initiated and special methods have to be used to start the solution process. Effectively an initial stability has to be provided to the solution process by some means and several methods can be used to achieve this (e.g., viscous relaxation and iterative procedures). In the viscous relaxation method an artificial dynamic problem is solved in which some inertial or viscous resistance provides stability to the solution process. The iterative method involves imposing an artificial initial stiffness to help initiate a stable incremental solution process. This artificial initial stiffness is gradually removed as the structure stiffens under the increasing loading. Subsequently, Newton-Raphson iterations are carried out to converge to the correct configuration under the given loading. 3.3 Numerical Examples The solution procedure described above is implemented in a general purpose computer program. The Fortran program can be used to analyse any three-dimensional framed structure subjected to an arbitrary loading. The user can specify the num-ber of load increments to be used and the sequence in which the Newton-Raphson iterations are to be performed. Hence for a given problem, the total load may be applied in a number of small increments and the iterations performed after each load increment. Alternatively, if the deformations are not excessive, it might be more efficient to apply the loads in a single step and then perform the iterations. How-ever, such a decision is purely problem dependent and, in general, the magnitude of deformations in a large displacement problem cannot be predicted beforehand. As discussed before, the convergence criterion for the Newton-Raphson iterations can be defined in terms of the imbalance between the internal and external loads. But a criterion based on displacement is preferable. For each structural degree of freedom r*, the ratio = Ar^/r^ is computed after each iteration. Convergence is achieved when the largest is less than a specified value of the tolerance, Et, which depends on the accuracy required. A value of Et = 10~ 4 is specified for all 50 the numerical examples presented here. The program is verified for a number of pertinent examples involving large deflections. In this section some typical results are presented. First, a cantilever subjected to a moment at the free end is considered. Next, the results are validated for a cantilever undergoing large three-dimensional deformations due to a specified bidirectional tip loading and axial torque. Finally, the static deflected shapes of typical compliant riser systems are obtained. Cantilever Beam with End Moment To begin with, a cantilever undergoing large rotations in one plane is consid-ered. Fig. 3.5 shows the problem specifications and a comparison of the computed tip deflections with values obtained independently using the ANSYS [51] program. Both results were obtained using ten elements of equal length. The results from the two finite element programs agree reasonably well. The discrepancy in the displacements predicted by the two programs decreased as the number of elements was increased. Cantilever Beam Undergoing Three-Dimensional Deformations Next, a vertical cantilever subjected to an axial torque and tip loading is con-sidered. The physical constants for the cantilever are presented in Fig. 3.6. It is subjected to a torque of 70 KNm and displacements of 3 m in both the Y and Z directions as shown. The results obtained for the displacements and rotations are shown in Fig. 3.6 while Fig. 3.7 presents the forces along the beam. Values obtained by O'Brien et al. [52] are also presented. The agreement between the two independent approaches is quite good. Here again, ten equal length elements were used. 51 f~7\ Moment = 105 kg - cm A = 10 cm 2 Ix — Iy = 5 cm4 E = 2.0x 106 kg/cm 2 L = 100 cm '///////////A TIP DEFLECTION X-Displacement (cm) Z-Displacement (cm) Y-Rotation (radians) LINEAR 0.00 50.00 1.00 ANSYS (Nonlinear) -15.82 45.99 1.00 PRESENT (Nonlinear) -15.72 46.07 1.00 Figure 3.5 Cantilever beam subjected to moment at the free end. 52 T 70 KNm DISTANCE ALONG X-AXIS (m) Figure 3.6 Cantilever beam subjected to torque and specified tip displacements. 600.0-, D I S T A N C E A L O N G X - A X I S (m) Figure 3 . 7 A comparison of bending moments and torque. 54 Riser Static Analysis Finally, the static analysis of typical compliant riser systems is undertaken. The static equilibrium configuration of a riser, under the dead loading, i.e., self-weight and buoyancy, is obtained. First, a "steep-S" type riser configuration is analysed. The riser consists of a flexible hose hanging from a surface vessel and a subsurface buoy as shown in Fig. 3.8. The problem involves predicting the deflected shape of the riser of any given length L when it is suspended between the surface vessel and the subsurface buoy. The physical properties of the riser are given in Table 3.1. Table 3.1 Data for "steep-S" riser configuration Water depth 150 m Distance from sea floor to 50 m subsurface buoy Riser pipe outer diameter 0.353 m Riser pipe inner diameter 0.254 m Modulus of elasticity of riser pipe 8.0 x 108 N/m? Weight per unit length of riser in air 208 kg/m Density of sea water 1025.18 kg/m3 Density of oil 919.5 kg/m3 The subsurface buoy is situated 50 m above the seafloor and the surface vessel is offset 100 m from the subsurface buoy in a water depth of 150 m. The computation involves three main steps as shown in Fig. 3.9. First, the initial configuration of the riser is taken as a straight line of length L. The starting point of the straight line is taken at the subsurface buoy while the Figure 3.8 Typical "steep-S" riser configuration. surface vessel STEP III end coordinates moved to correct position final initial position position subsurface buoy Figure 3.9 Static solution strategy. 57 terminating point is offset horizontally from the location of the surface vessel. Next, the deflected shape of this initial straight line configuration is obtained by subjecting it to the specified weight and buoyancy (keeping the positions of the two end points fixed). An artificial stiffness in the form of a fictitious tensile force is imposed on the straight line configuration to start the incremental loading procedure (see section 3.2). This stiffness is gradually removed over a number of load increments. Finally, the position of the terminating point is shifted to the correct location (at the surface vessel) in a number of equal increments. Newton-Raphson iterative procedure is applied after each increment to converge to the equilibrium position. Figure 3.10 shows the deflected shape of three different riser lengths (L = 145,155,165 m) obtained using the above procedure. It is noted that the computational effort involved can be drastically reduced by assuming that the initial shape is a catenary (instead of a straight line). Owen and Qin [53] obtained results for the same set of data using a simple analytic approach. They neglected the bending stiffness and solved the catenary problem. Their results match the present curves within plotting accuracy. The results show that the bending stiffness has negligible effect on the static deflected shape for the riser configuration under consideration. However, effect of the fiexural rigidity cannot be ignored in general. This is particularly true for some types of proposed compliant riser systems. Such problems are not amenable to analytic solution and are ideally suited to numerical treatment using the finite element method presented here. In addition, certain riser configurations also rely on variable buoyancy for their structural integrity. An example of this type is presented here as shown in Fig. 3.11. This is referred to as a "steep-wave" type compliant riser configuration. The 58 surface vessel 100.0 80.0 -20.0 H .20.0 _^-4U0 60.0 80.0 X C O O R D I N A T E (m) 'subsurface buoy 100.0 -40.OH 7Z^V Sea Floor V77^ 777^ 777^T 777^ -60.0 J Figure 3.10 Static deflected shape as affected by the riser length. 59 Surface Vessel Figure 3.11 Typical "steep-wave" riser configuration. 60 riser is located in 289 m depth of water and the surface vessel is offset 225 m from the point of attachment at the sea floor (as shown in Fig. 3.11). The riser is 420 m long and has two sections with additional mass and buoyancy. Buoyancy section I starts 30 m from the lower end and is 60 m long while buoyancy section II extends from this for another 30 m. The additional mass and buoyancy for these sections are listed in Table 3.2 together with the riser properties. Table 3.2 Data for "steep-wave" riser configuration Water depth 289 m Riser pipe outer diameter 0.375 m Riser pipe inner diameter 0.275 m Modulus of elasticity of riser pipe 8.0 x 108 N/m 2 Mass per unit length of riser in air 160 kg/m Additional mass in section I 222.4 kg/m Additional mass in section II 129.3 kg/m Buoyancy per unit length 599.4 N/m Additional buoyancy in section I 4970.3 N/m Additional buoyancy in section II 2889.7 N/m The static deflected shape for this configuration is shown in Fig. 3.12. The procedure used to arrive at the solution is similar to the one outlined above. Note that the additional buoyancy in the lower section of the riser keeps it clear off" the sea floor. The numerical examples presented here illustrate the large displacement capa-bility of the present computational procedure. In addition, the inherent advantages of the finite element method makes the present formulation general enough to be applicable to a large class of structural systems including configurations of com-300.0 surface vessel 250.0-^ 200.0 2 Q o o o ISI 150.OH 100.0 50.0 0.0 50.0 100.0 150.0 200.0 250.0 X C O O R D I N A T E (m) Figure 3.12 Static deflected shape of the "steep-wave" riser configurati 62 pliant riser systems proposed for deep water production facilities. In particular, complex geometries with variable or lumped mass, stiffness and buoyancy can eas-ily be handled and the effects of soil-stucture interaction can be incorporated using equivalent stiffness representation of the soil properties. 4. DYNAMIC ANALYSIS This chapter studies the dynamic response of typical riser systems when sub-jected to excitations due to ocean currents, waves and surface vessel motions. The dynamic analysis is necessary to ascertain the safe operational envelope of the drilling or production facility. It is used to calculate the maximum excursions of the riser displacements and stresses for a given sea state. This decides the limits of the environmental conditions beyond which the operations are stopped and the system is disconnected. Even if the extreme wave conditions are never encountered, offshore structures have to be safeguarded against fatigue failure. The peak stresses may never be reached but many reversing stress cycles imposed on a relatively high stress level can cause failure. This is particularly of concern for marine risers due to the continuous flexing of the pipe in the natural ocean environment. Hence dynamic analysis is needed to ensure the structural integrity of the system for the duration of its design life. The equations of motion used for the dynamic analysis are presented and a nu-merical solution procedure based on the finite element method developed. The drag force in the Morison equation representing the fluid loading is linearized and a fre-quency domain analysis performed. Next, the full nonlinearity of the hydrodynamic loading is retained and the equations are solved in the time domain. 4.1 Equations of Motion It is assumed that the dynamic loading is such that the displacements about the static equilibrium position are small. Hence, the linearized form of the equations of motion derived in Chapter 2 are adequate for the analysis. The equations used 63 64 for the finite element discretization are written as: .d 2u d 2u (m + PiAi)—-^ + 2piAiUi dt 2 rx 1 ldtdz o «d 2u d 4u + (piAiU? ~Te)-^ + EIy— = fx ; xd2v d 2v (m + PiAi)—-¥ + 2piAiUi dt 2 yx 1 xdtdz d 2v _ _ 6 dz^ + EIxdz 4 + (PiAiU?-Te)d^ + EIx^ = fy; d 2w 3Te _ 171 dt 2 dz ~ fz 5 d 26 d 26 In addition to the usual terms associated with the inertia forces, and the flexural, axial and torsional stiffnesses, effects of the static fluid pressures and internal flow are included in the equations (Eq. 4.1). Influence of the internal and external static fluid pressures is included by using the effective tension, Te (section 2.3). Internal flow gives rise to additional terms in the equations of motion in the transverse directions. These correspond to damping type forces associated with 2p,AtZ7,^^ and 2piAiUi4fejz- terms and Coriolis type forces which effect the stiffness of the 2 2 system through piAiU 2^^ and piAiU 2^^ terms. It is noted that although the equations are linearized in the sense that a fixed geometry is considered, some nonlinear effects are included by retaining the "geometric" stiffness. The finite element discretization of equations (4.1) is a relatively straightfor-ward procedure. The conventional beam element possessing six degrees of freedom per node (Fig. 4.1) is used. The axial displacement and the twist are linearly interpolated within the el-ement while a cubic shape function approximates the transverse displacements. 65 Figure 4.1 Conventional twelve degree-of-freedom beam element. 66 The discretization procedure leads to the matrix representation of equation (4.1) in terms of the nodal degrees of freedom in the local (x, y, z) coordinates. The ele-mental mass, damping and stiffness matrices so derived are presented in Appendix III. Using standard finite element techniques [54] the elemental matrices are trans-formed to the global (XYZ) coordinates and assembled in the usual manner to give a system of equations which can be represented in the form, Mr + Cr+Kr = / , •••(4.2) where: r = vector of nodal variables; / = external force vector; M = mass matrix including added mass effect; C = overall damping matrix, C_ ,• + C_ 8 ; Qi = internal flow damping matrix; cB = structural damping matrix; K = overall stiffness matrix, K r„ + K g\ KL = linear stiffness matrix; KG = geometric stiffness matrix. The structural damping is due to the friction within the riser material or at the connections between the elements of the system. The resulting damping forces are generally nonlinear functions of the strain or deflections. It is customary to replace the structural damping forces in a system by an equivalent viscous damping causing the same amount of energy dissipation [55]. Rayleigh damping [56] is often used in structural dynamics where the damping matrix is represented in the form: CB = 01M + P2K. •••(4.3) 67 The coefficients (5\ and are obtained by specifying the damping ratios ft and & in any two modes, where w,- is the natural frequency in the t th mode which is obtained from an eigen-value analysis. The actual level of structural damping present in marine riser systems is rather unclear. No experimental work in this area has been reported and typically a damp-ing ratio of anywhere upto 3% of the critical damping is generally used. However, the choice of an appropriate value for the damping ratio is not very crucial due to the substantial amount of hydrodynamic damping present. The hyrodynamic damping is introduced through the force vector / (Eq. 4.2) which is evaluated us-ing the Morison equation. The loading is a function of both the structural velocity and the fluid velocity due to the nonlinear drag force component. The presence of the nonlinear fluid loading is a major consideration in the choice of a solution technique for equation 4.2. The equations of motion can be numerically integrated through discrete time-steps. In this way the nonlinearities can be handled without any modification but this is computationally expensive. A solution in the frequency domain involves relatively shorter computer runs but the equations of motion have to be linearized. Both these approaches are presented here. 4.2 Solution in Frequency Domain Because of the random nature of the ocean waves a sea state is often specified in terms of known statistical properties (e.g., a significant wave height and a peak period). Then a probabilistic approach is used to obtain the structural response statistics of interest. Generally, the equations of motion are solved in the frequency 68 domain to use the powerful spectral analysis methods of random vibrations. Obvi-ously, a frequency domain solution is also an efficient tool to evaluate the structural response due to a deterministic wave of specified height and period. But a frequency domain solution technique requires the equations of motion to be linear. In what follows a linearized solution of the equations of motion in the frequency domain is presented. 4.2.1 Linearization of the drag force The equations of motion (Eq. 4.2) are nonlinear through the hydrodynamic loading which is evaluated using the Morison equation. Different linearization tech-niques can be used [32] to facilitate the solution of the equations of motion in the frequency domain. The procedure described here is often called the method of equivalent linearization. The nonlinear drag force in the Morison equation can be written as where A = \p0DeCd, and the relative velocity U r = U w — U p, in the absence of a current. The drag force is expressed in a linear form, where a and b are constants. These constants are chosen such that the resulting linear equation, in some sense, represents the original relation as closely as possible. One way to do this is to minimize the mean square error of the linearization. The mean square error, c, is fD=AUr\U r\, • ••(4.5) fD = a + bUr, • • • (4.6) e = {AU r\U r\-{a + bU r)} 2, • ••(4.7) where the bar denotes the average over a period. The linearization constants a and 69 b are obtained by minimizing e with respect to a and 6: da The relative velocity can be expressed in the harmonic form U r = Bsin(wf). • • • (4.8) In the case of regular waves, the linearization constants have the value: a = 0 AB. • • • (4.9) Note that the values a and 6 depend on the response velocity through the parameter B. Therefore an iterative procedure is required starting with an initial choice for the coefficients and the iterations continued until their values converge. As mentioned before, the analysis presented here is for a deterministic wave of specified height and period. For the case of random forcing the resulting expressions for the linearization constants will be different [6]. The linearization procedure described above was presented by Lipsett [32] for a single degree of freedom system. Here it is extended to a continuous system through the finite element discretization process. The linearization of the drag force given by equation 4.6 leads to an additional damping term in the discretized equations of motion (Eq. 4.2). Once the nonlinear drag force is linearized a solution can be easily obtained in the frequency domain. The excitation is due to a harmonic wave of frequency w. Linear wave theory is used to obtain the wave kinematics needed to evaluate the forcing terms (Appendix I). The periodic forcing function so obtained can be expressed as, The solution of the linear problem can also be represented in the periodic form f = fe 70 giving a set of algebraic equations as u 2M + iuC + iuCea + K f = / . •••(4.10) This system of complex equations is banded but not symmetric due to the nature of the internal flow damping matrix (Appendix III). The hydrodynamic damping that arises from the equivalent linearization procedure is represented as C_eq. This equivalent damping matrix, C_eq, depends on the value of the linearization con-stant b at each node. Hence the equations are solved iteratively till the values of the linearization constants at every node do not change. To initiate the iterative procedure it can be assumed that the structural response velocity is zero. But to facilitate rapid convergence to the final solution a reasonably accurate initial esti-mate of the velocities is obtained by setting the structural damping ratio to 15% for the first pass. This reduces the computational effort and generally leads to con-vergence within only two or three iterations for forcing frequencies away from the structure resonant frequencies. Up to ten iterations may be necessary at or around the resonant frequencies. , 4.2.2 Numerical examples The frequency domain solution procedure described above is implemented through a computer program which is capable of analysing the three-dimensional response of any general riser configuration. The loading can be given by a specified incident wave and noncolinear current. The surface vessel motion is assumed to be known (i.e., the amplitude and phase lag with respect to the incident waves are specified). In this section some typical numerical results are presented. First, the computational procedure is validated by comparison with known results for a conventional riser configuration. Next, the response of a typical drilling riser is presented. Finally, the effect of internal flow and structural damping on riser 71 response is studied. Computational Validation The American Petroleum Institute (API) Committee on the Standardization of Offshore Structures has defined a set of standard test cases as a basis of comparing the performance of riser analysis methods. API specified the input data for these test cases and eight participants were invited to submit solutions [57]. The purpose being to present data which can be used to validate other computational procedures. Figure 4.2 shows a schematic of the vertically aligned, drilling riser configuration of the API cases. The present formulation was validated by comparison with these API test cases. Results for one of the test cases are presented here. Table 4.1 presents the input data specifications for the test case designated as 500-20-1-D in the API bullettin [57]. Table 4.1 Data for API test case Distance from mean sea level to riser support ring Distance from seafloor to LBJ Water depth Riser pipe outer diameter Riser pipe inner diameter Modulus of elasticity of riser pipe Weight per unit lenght of riser joint Density of sea water Density of mud Drag coefficient Mass coefficient Effective hydrodynamic diameter Current at surface* Surface vessel static offset Wave height Wave period Vessel surge amplitude Vessel surge phase angle 9.144 m 152.4 m 0.4064 m 0.3747 m 2.07 x 10u N/m 2 256.59 kg/m 1025.18 kg/m 3 1438.46 kg/m 3 0.7 1.5 0.6604 m 0.2574 m/s 4.572 m 6.1 m 9s 0.61 m 15 deg 15.24 m * (varying linearly to zero at LBJ) F i g u r e 4.2 Typical conventional riser configuration. 73 In this case a tension of 54.43 t was applied at top of the riser. Figures 4.3 and 4.4 show the envelope of the maximum and minimum of the deflected shape and bending stress distribution. The API values plotted are the average of the combined results as obtained by the eight independent investigators. The comparison gives an indication of the agreement between the other methods and the analysis presented here. Next, the response of a similar drilling riser system in 453 m water depth is considered. The input data is the same as for the API case considered except the top tension is now 1311 and the surface vessel has a static offset of 13.5 m. The large top tension is needed to maintain the structural integrity of the riser. Figures 4.5 and 4.6 show the predicted extremes of the deflected shape and the bending stress distribution. The bending stresses near the sea floor are much greater than those near the surface due to the greater tension in the top part of the riser. The large bending stresses in the lower part of the riser can be reduced by using intermediate buoyancy modules which would increase the axial tension and thus the bending stiffness of the riser. Effect of Internal Flow The purpose of a production riser is to convey oil or gas from the well head to the surface vessel. It is known that internal flow affects the dynamic characteristics of a pipe and could even lead to buckling as well to as oscillatory instabilities at sufficiently high flow velocities. Internal flow has been incorporated in the present formulation so that its effect on marine riser dynamics can be evaluated. Here the effect of internal flow on a typical production riser is analysed. For simplicity the data specified in Table 4.1 is used again but the riser is subjected to wave action only (no current or surface vessel motion). The top tension is taken as 45.36 t and the internal fluid density pt=919.5 kg/m3. Figure 4.7 shows the DISPLACEMENT (m) Figure 4.3 A comparison of the displacement with the API data. 75 BENDING STRESS (N/mm2) Figure 4.4 A comparison of bending stress with the API data. 7 6 500.0-1 Figure 4.5 Typical displacement results for a drilling riser in 453 m water depth. 77 500.0 -i Figure 4.6 Typical bending stress results for drilling riser in 453 m water depth. 78 maximum deflection of the riser subjected to the specified incident waves (wave height=6.1 m). The response is shown as a function of wave period for different internal flow velocities (Z7,-=0, 7.32 and 14.63 m/s). It is seen that increasing the flow rate increases the response. As discussed in section 4.1, internal flow affects the system by decreasing the stiffness and also provides a type of negative damping. This leads to a shift in the response curves. The high flow velocities chosen here illustrate the effect of increasing flow rates on the system response. Furthermore, the present example is such that the funda-mental period of the riser lies within the range of typical wave periods. Ofcourse the effect of internal flow is more pronounced near the resonance peak. If the system natural period does not lie within the range of exciting wave periods the effect of internal flow on the system response will be much smaller. Hence it appears that internal flow has a negligible effect on the dynamical response of typical production risers. Effect of Structural Damping A damping ratio of 2% of the critical damping in the first two modes is specified for all the numerical examples presented so far. As discussed in section 4.1 the exact level of structural damping is not known, and is not expected to effect the riser response prediction by an appreciable amount. This effect is illustrated in Fig. 4.8 which shows the response of the API drilling riser subjected to wave action. The response is evaluated for two damping ratios (f = 0% and f = 2%). It is seen that the structural damping has only a small effect on the response at the resonance peaks. This is to be expected due to the substantial amount of hydrodynamic damping acting on the system. 79 2.0 - i w o a, x 1.5-• A • WAVE HT=6.1 m — NO INTERNAL FLOW A = FLOWRATE=7.32 m/s • = FL0WRATE=14.63 m/s 1.0 10.0 15.0 - i 1 r 20.0 1 25.0 PERIOD (s) F i g u r e . 4.7 Effect of internal flow on the riser response. 80 w u DM X 3.On 2.5 A 2.0 H 1.5 H 1.0H W A V E H T . = 6.1 m 2% 0 % 0.5 H 0.0 0.0 5.0 10.0 i 1 1 r 15.0 I 1 r — ' — I 20.0 25.0 P E R I O D (s) Figure 4.8 Effect of structural damping on the riser response. 81 4.3 Solution in Time Domain The frequency domain solution can only be obtained if the nonlinear drag force component of the Morison equation is linearized. The nonlinearity of the hydro-dynamic loading can be handled without modification by numerically integrating equation (4.2) in the time domain. This also allows the input of irregular wave sequences. But the solution in the time domain is computationally expensive and is generally used for evaluating the response of a structure when it is subjected to a deterministic loading only (as opposed to stochastic). Numerical time step integration procedures are broadly classified as explicit or implicit. In explicit methods, the solution at the end of the time-step is obtained by approximating equation (4.2) at the start of the time-step. On the contrary, in implicit methods, the solution corresponding to the end of the time-step is obtained by evaluating the equations also at the end of the step. Explicit time integration methods are easier to program, however, one has to be careful to choose a small enough time-step size, to ensure numerical stability. Whereas in implicit methods, the maximum time-step size that can be used is governed by the accuracy of the so-lution and not by the stability of the integration procedure [58]. Krieg and Key [59] have studied different methods to be used in conjunction with finite element models. They concluded that explicit, conditionally stable, central difference time integra-tion used with a diagonal mass matrix provides a practical means of computing the dynamic response. This method is chosen for the present analysis. 4.3.1 Explicit time step integration In the central difference method, it is assumed that the acceleration at any time t is given by It = [r t_ A t - 2r t + rt+At]/At2 , • • • (4.11) 82 and the velocity is written as, t t = [-Et-At + Lt+At]/2A* , •••(4.12) where At is the time-step used. If the displacements and velocities at time t are known, the accelerations can be evaluated using equation (4.2). Then, the displacements and velocities at time (t + At) can be obtained by using equations (4.11) and (4.12), respectively. Hence starting with some initial conditions for the problem this procedure can be continued till the desired time limit. The general expression that is evaluated at each time-step can be expressed in the form Note, no matrix inversion is required if the mass matrix M is diagonal. Moreover, there is no need to assemble the global matrices which leads to substantial reduction in computer memory requirements. 4.3.2 Numerical examples Here some results obtained using the time step integration procedure are pre-sented. First, the effect of retaining the full nonlinearity of the hydrodynamic drag force is assessed by comparing results with the frequency domain solution. Next, the effect of different combinations of loading on the response of a typical compliant riser configuration is studied. The choice of an appropriate time-step size for the numerical integration is rt+At = At2M- 1[f_-Krt] -AtM^Clrt-rt-At] + 2 r t - r t _ A t • • • • (4.13) 83 important. As mentioned above, the central difference method is conditionally stable, i.e., At must be kept below a critical time-step or the solution becomes unstable [58]. For this, At must be less than the time it takes for an elastic wave (flexural, torsional or axial) to travel across an individual element. A typical value of the time step-size used for the numerical results presented here is At = 0.001 s. Also, for the following results the structural damping is set to zero (f = 0). Effect of Linearization of Drag Force Effect of linearization of the drag force has been discussed thoroughly by Lipsett [32] for a single degree of freedom system. Inherent in the work is the assumption that his conclusions can be extended to a continuous system. To study the effect of linearization on the accuracy of the response predictions for a continuous system, the example used in section 4.2.2 to evaluate the effect of structural damping is considered again. Results are obtained using the time-step integration and frequency domain solution procedures. The numerical integration solution is used to assess the accuracy of the method of equivalent linearization. Figures 4.9 and 4.10 show the maximum riser displacement for two different wave heights (H = 6.1 m and H = 12 m). The wave height represents the magnitude of the loading and the importance of the nonlinearity. Note that the method of equivalent linearization approximates the response reasonably well except for the larger wave periods. For larger wave periods it underpredicts the response. This discrepancy increases with increasing nonlinearity. These trends in the results are analogous to the conclusions reached by Lipsett [32] for the single degree of freedom system. Effect of Different Load Combinations The "steep-S" compliant riser configuration defined in Chapter 3 is used to illustrate the effect of different combinations of loading on riser response. The 84 P E R I O D (s) Figure 4.9 A comparison between results obtained by time stepping and equivalent linearization methods (H = 6.1m). P E R I O D (s) Figure 4.10 A comparison between results obtained by time stepping and equivalent linearization methods (H = 12 m). 86 static deflected shape of the riser of length L — 155 m (see section 3.3) is taken as the reference configuration about which the dynamic analysis is performed. The nodal coordinates and forces obtained from the static analysis are used as input for the dynamic program and the system analysed in the time domain. The riser is subjected to the action of waves, current and surface vessel motion. All the results presented here are for a wave height of 12 m and a period of 18 s. Figure 4.11 shows the dynamic response of the riser subjected to wave action only. The plot shows the envelope of the maximum and minimum X-displacements of the riser about its static position. Figure 4.12 is a similar plot of the riser response under the combined action of waves and a current profile varying linearly from 1 m/s at the surface to zero at the sea floor. Next, the riser was subjected to the specified wave condition and a prescribed surface vessel motion (amplitude = 2 m, phase lag = 15 deg). The response of the riser is plotted in Fig. 4.13. Note that in these plots the dynamic displacements are magnified for clarity. 4.4 Closure In this chapter the equations of motion for the dynamic analysis of marine risers including the effect of internal flow are presented. The implementation of the computational procedure is illustrated by typical numerical examples. Although results are presented here for riser systems, it is readily seen that the method is also applicable to a wider variety of ocean based structures, for example, jacket platforms, articulated towers, marine pipelines, etc. In particular, OTEC cold water pipes and ocean mining systems also include internal flow. Also, for such systems the interaction between the pipe and the surface vessel cannot be ignored. 87 X C O O R D I N A T E (m) Figure 4.11 Minimum-maximum displacement envelope of compliant riser subjected to specified incident waves. 100,0 -i 80.0-w 6 0 0 H—t 2 40.0 O o t s : 20.0 Static deflected shape Envelope of min. and max. disp. Disp. enlarged by a factor of 5 Wave height = 12 m Wave period = 18 s Linearly varying current profile from 0 to 1 m/s 0.0 0.0 20.0 40.0 l 60.0 80.0 l 100.0 X C O O R D I N A T E (m) Figure 4.12 Minimum-maximum displacement envelope of compliant riser exposed to specified incident waves and current. 89 100.0 80.0-2 •—I S 40. O o o ISl 60.0-0 20.0-0.0 Static deflected shape Envelope of min. and max. disp. Disp. enlarged by a factor of 5 Wave height = 1 2 m Wave period = 18 s Surface vessel surge amp. = 2 m Phase lag = 15 deg 0.0 20.0 40.0 60.0 80.0 —I 100.0 X C O O R D I N A T E (m) Figure 4.13 Minimum-maximum displacement envelope of compliant riser due to specified incident waves and surface vessel motion. 90 This interaction can be easily included in the analysis if the response amplitude operators for the surface vessel are known. 5. CONCLUSIONS Some aspects of marine riser analysis have been looked at. Specifically, the problem of a pipe conveying fluid and undergoing large deformations has been addressed. The emphasis has been to look at the fundamental aspects involved in the analysis and in the process to develop tools which are applicable to a large class of problems. The equations of motion for a pipe conveying fluid and undergoing large three-dimensional deformations are derived. The nonlinear curvature expressions for a beam are obtained for the case of finite rotations. These are used to express the elastic strain energy of the beam and Hamilton's principle is applied to derive the equations of motion. The effect of internal flow is included in the variational formulation in a consistent manner. The nonconservative forces arising due to the fluid flow across the system boundaries are properly accounted for. It is shown that neglecting the contributions of the external fluid forces leads to improper boundary conditions. The problem of predicting the static deflected shape of compliant riser sys-tems has been addressed next and a finite element formulation for the analysis of a structure undergoing large deformations presented. It is based on the method developed by Nathan [56] which involves the use of a convected coordinate sys-tem moving with the elements. The nonlinear problem is solved using a combined incremental-iterative scheme. Typical numerical results illustrating the large dis-placement capability of the formulation are presented and it is concluded that the static equilibrium configuration of different types of riser systems can be predicted very accurately. Finally, the dynamic response of typical marine riser systems due to excitations of the ocean waves, currents and surface vessel motions is studied. The hydrody-91 92 namic loading is represented by a general form of the Morison equation. The nonlin-ear drag force in the Morison equation is linearized using the method of equivalent linearization and a solution procedure in the frequency domain developed. Alter-nately, the complete nonlinearity of the hydrodynamic loading is retained and the equations of motion are numerically integrated using an explicit time-step integra-tion scheme. The effect of the linearization on riser response predictions is evaluated. The method of equivalent linearization approximates the response reasonably well except for low frequency excitations. The effect of internal flow on the response of typical riser systems is also studied. Internal flow tends to reduce the stiffness of the system and gives rise to a negative damping mechanism. Hence it tends to have a destabilizing effect on the system response. However, realistic internal flow rates have a marginal effect on riser response. 93 REFERENCES 1. Chateau, G.M., "Oil and Gas Production Facilities for Very Deep Water," Pro-ceedings of the Third International Behaviour of Offshore Structures Confer-ence, BOSS'82, Boston, Massachusetts, 1982, Vol. I, pp. 50-70. 2. de Oliveira, J.G., Goto, Y., and Okamoto, T., "Theoretical and Methodolog-ical Approaches to Flexible Pipe Design and Application," Proceedings of the Seventeenth Offshore Technology Conference, Houston, Texas, 1985, OTC 5021, Vol. 3, pp. 517-526. 3. Panicker, N.N., and Yancey, I.R., "Deepwater Production Riser," Proceedings of the Fifteenth Offshore Technology Conference, Houston, Texas, 1983, OTC 4512, Vol. 2, pp. 9-17. 4. Chakrabarti, S.K., and Frampton, R.E., "Review of Riser Analysis Techniques," Applied Ocean Research, 1982, Vol. 4, No. 2, pp. 73-90. 5. Narzul, P., and Marion, A., "Static and Dynamic Behaviour of Flexible Cate-nary Risers," Proceedings of the Fifth International Offshore Mechanics and Arctic Engineering Symposium, Tokyo, Japan, 1986, Vol. Ill, pp. 378-386. 6. Sarpkaya, T., and Isaacson, M., "Mechanics of Wave Forces on Offshore Struc-tures," Van Nostrand Reinhold, New York, 1981. 7. Chakrabarti, S.K., "Hydrodynamics of Offshore Structures," Springer-Verlag, New York, 1987. 8. Dym, C.L., and Shames, I.H., "Solid Mechanics," McGraw-Hill, New York, 1973. 9. Patel, M.H., Sarohia, S., and Ng, K.F., "Finite-Element Analysis of the Marine Riser," Engineering Structures, 1984, Vol. 6, pp. 175-184. 10. Wang, E., "Analysis of Two 13,200 ft. Riser Systems Using a Three Dimen-sional Program," Proceedings of the Fifteenth Offshore Technology Conference, Houston, Texas, 1983, OTC 4563, Vol. 3, pp. 435-444. 11. Verbeek, P.H.J., "Analysis of Riser Measurement in the North Sea," Proceed-ings of the Fifteenth Offshore Technology Conference, Houston, Texas, 1983, OTC 4562, Vol. 2, pp. 425-434. 12. Bernitsas, M.M., "Contribution Towards the Solution of Marine Riser Design Problem," Ph.D. Thesis, 1977, Massachusetts Institute of Technology, Boston, U.S.A. 13. Kim, Y.C., "Nonlinear Vibrations of Long Slender Beams," Ph.D. Thesis, 94 1983, Massachusetts Institute of Technology, Boston, U.S.A. 14. McNamara, J.F., and Lane, M., "Practical Modeling for Articulated Risers and Loading Columns," ASME Journal of Energy Resources Technology, 1984, Vol. 106, pp. 444-450. 15. Chucheepsakul, S., "Large Displacement Analysis of a Marine Riser," Ph.D. Thesis, 1983, University of Texas at Arlington, Arlington, U.S.A. 16. McNamara, J.F., O'Brien, P.J., and Gilroy, S.G., "Nonlinear Analysis of Flexi-ble Risers Using Hybrid Finite Elements," ASME Journal of Offshore Mechan-ics and Arctic Engineering, 1988, Vol. 110, No. 3, pp. 197-204. 17. Chung, J.S., Whitney, A.K., and Loden, W.A., "Nonlinear Transient Motion of Deep Ocean Mining Pipe," ASME Journal of Energy Resources Technology, 1981, Vol. 103, pp. 2-10. 18. Felippa, C.A., and Chung, J.S., "Nonlinear Static Analysis of Deep Ocean Mining Pipe - Part I: Modeling and Formulation," ASME Journal of Energy Resources Technology, 1981, Vol. 103, pp. 11-15. 19. Mathisen, K.M., and Bergan, P.G., "Nonlinear Static Analysis of Flexible Ris-ers," Proceedings of the Fifth International Offshore Mechanics and Arctic En-gineering Symposium, Tokyo, Japan, 1986, Vol. Ill, pp. 337-353. 20. Vogel, H., and Natvig, B.J., "Dynamics of Flexible Riser Systems," Proceedings of the Fifth International Offshore Mechanics and Arctic Engineering Sympo-sium, Tokyo, Japan, 1986, Vol. HI, pp. 363-370. 21. Blevins, R.D., "Flow Induced Vibration," Van Nostrand Reinhold, New York, 1977. 22. Naguleswaran, S., and Williams, C.J.H., "Lateral Vibrations of a Pipe Convey-ing a Fluid," Journal of Mechanical Engineering Science, 1968, Vol. 10, No. 3, pp. 228-238. 23. Housner, G.W., "Bending Vibrations of a Pipe Line Containing Flowing Fluid," ASME Journal of Applied Mechanics, 1952, pp. 205-208. 24. Paidoussis, M.P., and Issid, N.T., "Dynamic Stability of Pipes Conveying Fluid," Journal of Sound and Vibration, 1974, Vol. 33, No. 3, pp. 267-294. 25. Benjamin, T.B., "Dynamics of a System of Articulated Pipes Conveying Fluid - Theory," Proceedings of the Royal Society (London), 1961, Series A, Vol. 261, pp. 457-486. 26. Benjamin, T.B., "Dynamics of a System of Articulated Pipes Conveying Fluid - Experiments," Proceedings of the Royal Society (London), 1961, Series A, Vol. 95 261, pp. 487-499. 27. Thurman, A.L., and Mote, CD., "Nonlinear Oscillations of a Cylinder Con-taining a Flowing Fluid," ASME Journal of Engineering for Industry, 1969, Vol. 91, pp. 1147-1155. 28. Wu, M-C, "Effects of Rigidity and Internal Flow Velocity on Marine Riser Dy-namics," Ph.D. Thesis, 1985, Texas A & M University, College Station, U.S.A. 29. Hwang, Y-L., "Normal Mode Analysis of the Dynamic Response of the OTEC Cold Water Pipe and Platform System to a Random Seaway," Ph.D. Thesis, 1979, University of California, Berkeley, U.S.A. 30. Isaacson, M., "Wave and Current Forces on Fixed Offshore Structures," Cana-dian Journal of Civil Engineering, 1988, Vol. 15, No. 6, pp. 937-947. 31. Ismail, N.M., "Effects of Wave-Current Interaction on the Design of Marine Structures," Proceedings of the Fifteenth Offshore Technology Conference, Hous-ton, Texas, 1983, OTC 4615, Vol. 3, pp. 307-316. 32. Lipsett, A.W., "Nonlinear Response of Structures in Regular and Random Waves," Ph.D. Thesis, 1985, University of British Columbia, Vancouver, Canada. 33. Brouwers, J.J.H., "Analytical Methods for Predicting the Response of Marine Risers," Proceedings of the Royal Netherland Academy of Arts and Sciences, 1982, Series B, Vol. 85, pp. 381-400. 34. Kirk, C.L., Etok, E.U., and Cooper, M.T., "Dynamic and Static Analysis of a Marine Riser," Applied Ocean Research, 1979, Vol. 1, No. 3, pp. 125-135. 35. Kim, Y.C., and Triantafyllou, M.S., "The Nonlinear Dynamics of Long, Slen-der Cylinders," Proceedings of the Third International Offshore Mechanics and Arctic Engineering Symposium, 1984, Vol. I, pp. 477-488. 36. Spanos, D. P-T., and Chen, T. N., "Vibrations of Marine Riser Systems," ASME Journal of Energy Resources Technology, 1980, Vol. 102, pp. 203-213. 37. Kirk, C.L., "Dynamic Response of Marine Risers by Single Wave and Spectral Analysis Methods," Applied Ocean Research, 1985, Vol. 7, No. 1, pp. 2-13. 38. Love, A.E.H., "A Treatise on the Mathematical Theory of Elasticity," Fourth Edition, Dover Publications, New York, 1944. 39. Novozhilov, V.V., "Foundations of the Nonlinear Theory of Elasticity," Gray-lock Press, Rochester, New York, 1953. 40. Finlayson, B.A., "The Method of Weighted Residuals and Variational Princi-96 pies," Academic Press, New York, 1972. 41. Sommerfeld, A., "Mechanics: Lectures on Theoretical Physics," Academic Press, New York, 1942. 42. Fox, R.W., and McDonald, A.T., "Introduction to Fluid Mechanics," John Wi-ley & Sons, New York, 1985. 43. Morison, J.R., O'Brien, M.P., Johnson, J.W., and Schaff, S.A., "The Force Ex-erted by Surface Waves on Piles," Petroleum Transactions of the American In-stitute of Mining and Metallurgical Engineers, 1950, Vol. 189, pp. 149-154. 44. Isaacson, M., "The Forces on Circular Cylinders in Waves," Ph.D. Thesis, 1974, Department of Engineering, University of Cambridge, Cambridge, U.K. 45. Sparks, CP., "The Influence of Tension, Pressure and Weight on Pipe and Riser deformations and Stresses," ASME Journal of Energy Resources Tech-nology, 1984, Vol. 106, pp. 46-54. 46. Nathan, N.D., "Finite Element Formulation of Geometrically Nonlinear Prob-lems of Elasticity," Ph.D. Thesis, 1969, University of Washington, Seattle, U.S.A. 47. Radomske, B.R., "Large Deflection Analysis of Shallow Framed Structures," M.A.Sc. Thesis, 1972, University of British Columbia, Vancouver, Canada. 48. Bathe, K.-J., "Finite Element Procedures in Engineering Analysis," Prentice-Hall, New Jersey, 1982. 49. Zienkiewicz, O.C., "The Finite Element Method," McGraw-Hill, London, 1977. 50. Leonard, J.W., "Tension Structures," McGraw-Hill, New York, 1988. 51. Peter, K., "ANSYS: Engineering Analysis System Theoretical Manual," Swan-son Analysis Systems, Inc., Houston, Pensylvania, 1986. 52. O'Brien, P.J., McNamara, J.F., and Dunne, F.P.E., "Three-Dimensional Non-linear Motions of Risers and Offshore Loading Towers," Proceedings of the Sixth International Offshore Mechanics and Arctic Engineering Symposium, Houston, Texas, 1987, Vol. I, pp. 171-176. 53. Owen, D.G., and Qin, J.J., "Non-Linear Dynamics of Flexible Risers by the Finite Element Method," Proceedings of the Sixth International Offshore Me-chanics and Arctic Engineering Symposium, Houston, Texas, 1987, Vol. I, pp. 163-170. 54. Cook, R.D., "Concepts and Applications of Finite Element Analysis," Second 97 Edition, John Wiley & Sons, New York, 1981. 55. Hurty, W.C., and Rubinstein, M.F., "Dynamics of Structures," Prentice-Hall, New Jersey, 1964. 56. Clough, R.W., and Penzien, J., "Dynamics of Structures," McGraw-Hill, New York, 1975. 57. American Petroleum Institute, "Comparison of Marine Riser Analysis," API Bulletin 23, First Edition, 1977. 58. Key, S.W., "Transient Response by Time Integration: Review of Implicit and Explicit Operators," Advanced Structural Dynamics, Edited by J. Donea, 1978, pp. 71-95. 59. Krieg, R.D., and Key, S.W., "Transient Shell Response by Numerical Integra-tion," Advances in Computational Methods in Structural Mechanics and De-sign, Edited by Oden, T.T., Clough, R.W., and Yamanoto, Y., 1972, pp. 237-257. APPENDIX I REVIEW OF LINEAR WAVE THEORY Linear wave theory forms the basis for computing the water particle kinematics for a given wave. Figure 1.1 shows the general form of a wave train propagating in the X-direction. The wave height H is the vertical distance from trough to crest, the wavelength L is the distance between successive crests, the wave period T is the time interval between successive crests passing a particular point and c is the speed of the wave travelling through the fluid (c = L/T). The governing equation is the Laplace equation in two dimensions, d2<j> d26 , , a ? + a ? = ° - - < " > where cf> is the velocity potential. 99 The boundary conditions that apply in linear wave theory are: dj>_ dz 86 dtj ~dz~ ~di d6 -dl + 9r} 0 at z = —d 0 at z = 0; 0 at z = 0: (7.2) (7.3) (7.4) where n is the free surface elevation. The boundary conditions expressed by equa-tions 7.3 and 7.4 are linearized in the sense that these boundary conditions are applied at the still water level (z = 0) instead of the actual water surface (z = n) and both these equations neglect nonlinear terms. The two free-surface boundary conditions are combined to give: d 26 d6 , , -o£ + 9 £ = 0 at z = 0; .-.(7.5) ld6 r/ = at z = 0. ••• (7.6) g dz Another condition that must be satisfied is: 6{x,z,t) =6{x-ct,z) ; •••(^•7) which describes the periodic nature of the wave train. The solution to this boundary value problem may be written in the form: ,/ x irHcosh{k(z+ d)\ . . , where the wave crest tj = 77/2 passes x = 0 at time t = 0, k is the wave number related to the wavelength L by k = 2n/L and to is the wave frequency related to the wave period T by CJ = 2%/T. The boundary condition expressed by equation 7.2 results in the dispersion re-lationship between the wave frequency and the wave number. Substituting equation 100 1.8 into equation 1.2 gives, u 2d = kd tanh kd , • • • (1.9) 9 which is a transcendental equation for the parameter kd. Once the solution for the velocity potential 6 is obtained the other variables of interest are easily determined. The horizontal and vertical velocity components, Ux and , of the fluid particles follow from the appropriate derivatives of 6: 86 nH cosh{k[z + d)} Ux = TT~ = . W , ,. COSIKX - ut) ; x dz T sinh(fcd) v ' „ „ . 86 ir H smh.{ k[z + d)\ ... . . . z dx T sinh(A:d) v ; v 7 And the acceleration components, and Z7]j, follow from the temporal derivatives of the velocities: • 27r 2Hcosh{k{z + d)} . ^ = -7^ sinh(M S 1" (f c* - " ' ) ; ^ _ 2 ^ s i n h W z + d)} ^ _ u t ) ... ( z T2 sinh(A:<f) K ] K ' APPENDIX II ELEMENT MATRICES FOR STATIC ANALYSIS The element stiffness matrices £ Q and kx defined in Chapter 3 are presented here: AE L 0 0 0 0 0 0 0 0 0 0 12EIV 0 0 0 0 GJ L 0 0 0 6EIV if 0 AEIy L 0 6EIZ 0 0 0 0 0 0 0 4EIZ 0 0 0 0 0 0 0 6AE BL 2 0 0 0 AE 10L 0 0 6AE 0 AE 10L 0 0 0 0 Ely 0 0 0 0 AE 10L 0 6AE 0 0 AE 10L 0 0 0 6AE TIT 101 k.l{ s2,) — 52 0 6AE 0 0 AE 10L 6AE TU 0 0 0 0 0 0 0 0 6EIV 0 0 0 0 —if 0 21 Ely 0 0 0 0 2 lEh -JiX 0 0 AE 10L 0 0 0 0 0 0 0 6AE 0 AE 10L 0 6AE 5L* 0 0 &EIy if 0 0 6AE 5L2 0 0 0 0 0 0 0 0 0 21 Ely -nf AE lOL 0 0 21 Ely nf 0 0 0 0 0 21 Ely ~nf 0 0 i(*2) = «! o AE 10L 0 0 2 1SAE 0 0 0 21 Ely ~nf o 0 0 AE 10L 6EIV V 0 0 0 21 Ely 0 0 0 IQEIy BL 15AE 0 0 0 0 0 0 0 0 16EIy 5L 0 0 103 0 AE 10L 0 0 0 2 15AE AE 10L 0 0 0 0 0 0 0 0 2 IE I. 0 0 y 0 0 21 Ely 0 16gJy 5L 0 0 0 0 16EIy SL 0 0 1 5 A £ 0 0 0 0 0 Where: Iz is the moment of inertia about the element z3 axis; Iy is the moment of inertia about the y 3 axis (Iy > Iz); J is the torsional moment of inertia; A is the element cross-sectional area; L is the length of the element; E is the modulus of elasticity; and, G is the torsional modulus. APPENDIX III ELEMENT MATRICES FOR DYNAMIC ANALYSIS The element matrices for the conventional twelve-degree-of-freedom beam ele-ment used in Chapter 4 are presented here. The degrees of freedom are as shown in Fig. 4.1. Linear Stiffness Matrix The contributions to the linear stiffness matrix corresponding to the different degrees of freedom follow. Axial deformation: kT = EA 1 -1 -1 -1 Torsional deformation: kT. = GJ 1 -1 -1 -1 Bending in x — z plane: 12 6L -12 6L Ely 6L 4L2 -6L 2L 2 L 3 -12 -6L 12 -6L 6L 2L 2 -6L L 2 104 105 Bending in y — z plane: 12 -6L -12 -6L -6L 4L2 6L 2L2 -12 6L 12 6L -6L 2L2 6L L 2 Inserting these stiffness contributions in to the appropriate locations gives the 12 x 12 element stiffness matrix: 1 0 0 0 0 0 -7 0 0 0 0 0 0 12a 0 0 0 6La 0 -12a 0 0 0 6La 0 0 12/3 0 -6L/3 0 0 0 -12)3 0 -6L/30 0 0 0 0 6 0 0 0 0 0 -6 0 0 0 0 -6L0 0 4L2/3 0 0 0 6L/3 0 2L20 0 0 6La 0 0 0 4L2a 0 -6La 0 0 0 2L2a -1 0 0 0 0 0 1 0 0 0 0 0 0 -12a 0 0 0 -6La 0 12a 0 0 0 -6La 0 0 -12/3 0 6L/3 0 0 0 12/? 0 6L/3 0 0 0 0 -6 0 0 0 0 0 6 0 0 0 0 -6L/3 0 2L2p 0 0 0 6L/3 0 4L20 0 0 6La 0 0 0 2L2a 0 -6La 0 0 0 4L2a where: EIy a EIX EA . GJ u EI* Geometric Stiffness Matrix Similarly, the contributions to the geometric stiffness matrix are expressed as. 106 Bending in x — z plane: 36 ZL -36 ZL r ZL 4L 2 -ZL -L 2 30L -36 -ZL 36 -ZL ZL —L 2 -ZL AL 2 Bending in y — z plane: 36 -ZL -36 -ZL -ZL 4L 2 ZL —L 2 -36 ZL 36 ZL -ZL —L2 ZL 4L 2 where T'= [piAiU?-Te] . These sub-matrices are assembled as before in to a 12 x 12 element geometric stiffness matrix. Internal Flow Damping Matrix Besides contributing to the geometric stiffness the internal flow also gives rise to a damping type matrix as discussed in Chapter 4. The contribution to the damping matrix is presented here. Bending in x — z plane: 2ptAiUi l 2 • A * _1 2 _1_ 10 A * 0 _ j _ 10 1 T 2 60 1 2 10 1 2 J. 10 1_ 10 6 0 ^ x 10 - A 0 107 Bending in y — z plane: c. 2PIAIUI _ l 2 _1_ 10 _ 1_ 2 __1_ 10 0 _i_ 10 6 0 ^ i_ 2 _1_ " 10 I 2 _1 10 10 6 0 ^ 1_ 10 iTT 0 Mass Matrix Finally, the mass matrices used in the dynamic analysis are also presented. Axial degrees of freedom: mL m = 2 1 1 2 where m is mass per unit length of the element. Torsional degrees of freedom: m _ IPL 2 1 1 2 where Ip is the mass moment of inertia per unit length of the element about the z axis. Bending degrees of freedom: m = mL 420 156 22L 54 -13L 22L 4L2 13L -3L2 54 13L 156 -22L -13L -3L2 -22L 4L2 This consistent mass matrix is used for the frequency domain analysis, how-ever as discussed in section 4.3, a diagonal mass matrix is used for the time step 108 integration procedure. The diagonal terms of this mass matrix are: m = / mL/2 \ 39mL/78 IpL/2 mL3/78 roL3/78 mL/2 39mL/78 IpL/2 mL3/78 mL3/78 J \
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Some aspects of marine riser analysis
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Some aspects of marine riser analysis Irani, Mehernosh Boman 1989
pdf
Page Metadata
Item Metadata
Title | Some aspects of marine riser analysis |
Creator |
Irani, Mehernosh Boman |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | The problem of the static and dynamic analysis of marine risers is considered in this thesis. The equations of motion for a pipe conveying fluid are derived using a variational approach. The nonlinear expressions for the bending and torsional curvatures are obtained for a beam undergoing finite rotations. The elastic strain energy is expressed in terms of these curvatures and Hamilton's principle used to derive the equations of motion valid for large deformations. The effect of internal flow is included in the variational formulation in a consistent manner. A finite element formulation is developed for the static analysis of compliant risers. Large deformations are handled by using a convected coordinate system fixed to the elements. The nonlinear problem is solved using a combined incremental-iterative approach based on an instantaneous linearization of the Taylor series expansion of the forces about the displaced configuration. The numerical results presented here show that the large displacement formulation can predict the static deflected shape of different types of compliant riser systems. The dynamic response of a typical marine riser system to excitations of the ocean waves, current and surface vessel motion is also studied. The hydrodynamic loading is represented by a general form of the Morison equation. The nonlinear drag force in the Morison equation is linearized using the method of equivalent linearization and a solution procedure in the frequency domain developed. Alternatively, the complete nonlinearity of the hydrodynamic loading is retained and the equations of motion are integrated in the time domain. The effect of this linearization on riser response predictions is evaluated. The method of equivalent linearization involves less computational effort and it approximates the response reasonably well. The effect of internal flow on the response of a typical production riser configuration is also studied. It is concluded that internal flow has a negligible effect on the dynamic response. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0098209 |
URI | http://hdl.handle.net/2429/29117 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1989_A1 I72.pdf [ 4.33MB ]
- Metadata
- JSON: 831-1.0098209.json
- JSON-LD: 831-1.0098209-ld.json
- RDF/XML (Pretty): 831-1.0098209-rdf.xml
- RDF/JSON: 831-1.0098209-rdf.json
- Turtle: 831-1.0098209-turtle.txt
- N-Triples: 831-1.0098209-rdf-ntriples.txt
- Original Record: 831-1.0098209-source.json
- Full Text
- 831-1.0098209-fulltext.txt
- Citation
- 831-1.0098209.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0098209/manifest