Natural frequencies , mode shapes and modal damping ratios must be estimated to assess subsea pipeline spans for vortex -induced vibrations (VIV) and respons e to direct wave loading . Several approximate solutions exist f or a linearly elastic pipe under constant axial force supported by linearly elastic springs beyond the span’s shoul ders. An exact analytical solution has only recently been published. That s olution is used here in a Rayleigh -Ritz approximation to account for arch action arising from combined effect s of sag under gravity loads and axial restraint at the shoulders . The method allows survey data to be used directly to quantify arch action . Its accuracy is confirmed by finite element analysis. Further , the modal damping ratio is estimated based on the fractions of the potential energy in bending, the axial force, and the soil springs , all of which are determined analytically . Thus, it is found that the effective modal damping ratio increases without a bound as the axial load approaches the buckling load in compression. 1 Introduction Spans for subsea pipelines can aris e from laying a pipeline on an uneven seabed, from current -induced scour ing, wave action and/or other erosive phenomena [ 1] after the pipeline has been laid. This c an lead to fatigue damage from vortex induced vibrations . Mostly this is assess ed by modal analysis, together with empirical response functio ns, which DNVGL have periodically improved over the years as more experiment al data became available [2]. Finite element analysis can be used to calculate the required natural frequencies and mode shapes. Such analyses involve a nonlinear analysis under the static loads ( gravity , internal pressur e, temperature, pipe lay), followed by the linearized eig envalue analysis to obtain the natural frequencies and mode shapes [3-5]. Here a n isolated span is considered, with the pipe continuously supported by soil beyond the span shoulders as shown in Figure 1. The effect of axial force is included, as it can strongly affect natural frequenc ies, especially if , due to the combined effects of internal pressure and operating temperature, the axial load approaches the buckling load. An analytical solution for this isolated span problem , including the soil springs, but not the effective axial force [6-7] has b een given by Hobbs [8]. His results highlight the influence of the transverse soil stiffness on the natural frequency. In the same paper Hobbs also introduced the c oncept of a n effective span length, as the length of a span with pinned boundary conditions at the ends for which the frequency matches that calculated based on his beam -on-elastic foundation analysis. Park and Kim [9] provided a solution to the same problem without accounting for the axial force or the inertia in the soil -supported zone . Fur ther, Park and Kim’s solution neglects the off-diagonal term in the 2 by 2 stiffness matrix that represents the stiffness from the soil -supported part of the pipe. Kapuria et al. [10] do include the off -diagonal term, but still not the effect of inertia or axial forces in the soil - supported zone, though the effect of the axial force in the spanning zone is included . Finally , Sollund et al. [11] provide an analytical solution including both transverse soil springs and the effect of axial force within the soil supported zone as well as the span. Further, Vedel d and Sollund [12] provide a “semi -analytical” solution based on Fourier modes for the transverse and axial displacements . It leads to coupled equations for the amplitudes of these Fourier modes. The method not only accounts for a xial force and inertia everywhere, but also for dynamic changes in the axial load associated with arch action and is readily extended to multiple interacting spans [ 13]. However, for 1 Peek Solutions, Sant Andreu de Llavaneres, Spain 2 Crondall Energy Subsea, Ltd., Aberdeen, Scotland. 3 DNVGL, O slo, Norway. 2 accurate solutions a considerable number of Fourier modes need to be included g iving rise to coupled equations for the amplitudes of these modes. To develop practical approximate solutions including the effect of the axial force , Fyrileif and M ørk [14] also use the concept of an effective span length, but their effective length applies for fully -fixed instead of pinned boundary conditions. Their effective span length to match the natural frequency (from the solution by Hobbs [8]) and that to match the buckling load (from a solut ion for this from Hetényi [15]) are much the same. On this result , they build approximate expression s for the natural frequency and modal stress including the influence of axial force and soil stiffness. They calibrate their approximations by comparison with the results of a range of finite element analysis results . Their results are incorporated in the DNVGL span assessment guidance [2] and much used . Here , the solution of Sollund et al. [11] accounting for axial force everywhere as we ll as inertia in the soil-supported zone is presented in a compact normalized form (Sections 3.1-3.4), evaluated covering all cases of pipeline engineering interest (Section 3.5) and compared with the approximate formulations in [2,14] that are extensively used in engineering practice (Section 3.6). The above solution for the axially unconstrained case (i.e. constant axial force) is then extended to account for arch action that arises from static sag of the span together with the axial re straint at the shoulders. (With the sag, the pipe profile is like that of an upside -down arch . Although it carries gravity loads in axial tension rather than compression, the stiffening effect , which applies for dynamic upward and downward load increments , is the same.) With axial restraint at the shoulders, dynamic changes in axial force develop not only due to the sag under gravity loads , but also due to nonlinear effects associated with large -amplitude oscillation s. (Cross-flow VIV amplitude can easily reach or exceed the static sag. ) Both t hese effects have been considered by Bruschi and Vitali [ 16]. A difficulty with nonlinear m odal analysis, is that separable solutions (i.e. solutions that can be written as a product of a function of time t and a function of location x) typically do not exist , except for special cases, such as the simply -supported span. Standard modal analysis, either by finite elements [ 17] or the Fourier modes of [ 12,13] employ linearization. Thus , they include the arch action due to sag under gravity loads, but not the nonlinearity due to large - amplitude vibrations. Here also only linear dynamic response is considered. For a method to include nonlinearity in the dynamic response of the soil springs, see [ 18]. Although other authors did not use the term “arch action” for pipeline spans, its importance has been recognized and illustrated in several studies [3-21,14]. Further, Fyrileiv and M ørk [14] provide an approximate formula for the effect of arch action on natural frequency that is incorporated by DNVGL in [2]. However , their results only apply for a certain value of the axial stiffness at the span shoulders and they do not state what this value is. Clearly, from [ 21], it does not apply for other values. Here, for the first time, a general formula for the effect of arch action on the natural frequency of spans is given which account s for axial stiffness from soil support a t the shoulders as well the static equilibrium profile of the span. The method is quite gen eral. It applies for any seabed bathymetry provided the span is and remains a single span and the slopes remain small compared to 1 . It is based on a Rayleigh Ritz approximation (Section 4) that is found to be excellent when c ompared to a finite element solution for an example problem presented in Section 5 that does not use the small -slopes approximation . The analytical solution presented also provides a new and effective way to improve damping estimates , which is especially important for horizontal (in -line) vortex -induced vibrations. This is done by breaking down the modal potential energy into 4 components: (1) from flexural deformation of the pipe, (2) from the (constant) axial force , a compo nent which is positive for axial tension but negative for axial compression , (3) from arch action (varying component of axial force), and ( 4) from the transverse soil springs. Energy dissipation is best estimated as a fraction of the potential energy of each of these components and combined as described in Section 6. 2 Nomenclature A = amplitude of modal displacement y(t); corresponds to amplitude of vibrations at mids pan for the first mode a,b = dimensionless real parameters that de fine roots of the characteristic equation in the inner zone 3 D = outer diameter of pipe incl uding any coatings E = Young’s modulus of the pipe Ed = energy dissipated per cycle at the undamped natural frequency due to structural, soil, and axial -force dampi ng EI = effective flexural rigidity including contribution of pipe coatings, if applicable EIFJ = flexural rigidity at the field joints (not includi ng the contribution of concrete coating that is not present there) f = natural frequency f̅ = normalized na tural frequency, f̅ =ω/ω 1 = 2πf/ω1, where ω1 is from Eq. (3) f1, f2, f3, f4 = basic symmetric solutions that decay with distance from the span; f 1 and f 2 apply within the inner zone |x| ≤ L/2, f 3 and f 4 apply in the outer zone |x| > L/2 g = amount by which the pipe becomes longer due to the transverse (vertical and /or lateral) displacements Iϕϕ, Iw’ϕ’, Iϕ’ϕ’, Iϕ’’ϕ’’ = integrals over the mode shape, i.e. over x ∈ [-xmode, xmode], of the expression appearing as a subscript, e.g. Iϕϕ = ∫mode ϕ2 dx, Iw’ϕ’ = ∫mode w,x ϕ,x dx, etc. k = stiffness of soil springs (force per unit deflection per unit length along the pipeline) ks, kL, kV = stiffness of soil springs for static (vertical), dynamic lateral and dynamic vert ical loading respectively Ka, Kbending , K N, K soil = contributions to modal stiffness from arching, bending, the axial force and the soil springs, respectively. Defines dynamic potential energy for these components. For instance, the change in potential energy in bending is due to a displacement yϕ from the static equilibrium position is ½ K bending y2. Ktu = total modal stiffn ess under axially unconstrained conditions, K tu = K bending + K N + K soil Kt = total modal stiffness, including the effect of arch a ction arising from the combination of static deflection and axial restraint at the shoulders of the span, K t = K tu + K a L = span length L̅ = normalized span length, L̅ = L/x 1 Lea = effective length in regard to axial restraint: chosen such that a length L ea of pipe axially fixed at the ends but not otherwise axially restrained provides the correct relation between the amount Δg by which the pipe becomes longer due to dynamic displacements, and the corresponding change in axial force ΔN. See Section 4. m = effective mass per unit length, including pipe, coatings, content, and hydrodynamic added mass (further details after Eq. (2) in Section 3) N = effective axial force , positive for tension N̅ = normalized effective axial force, N̅ = N/N 1 N1 = unit value of axial force, given in Eq. (3) Nc = absolute value of the axial force when the pipe buckles in compression N̅c = normalized buckling force , N̅c = N c/N1, shown in Figure 2 and Figure 3 R = radius to point at which the bending stress is determined t = time w = transverse stati c equilibrium displacement; typicall y refers to the vertical displacement but the same form ulation can be used for the lateral direction if the pipe has lateral out -of- straightness, e.g. due to a route bend x = coordinate defining location on the pipeline; zero at midspan 4 x̅ = normalized coordinate, x̅ = x/x 1 where unit value x 1 is given in Eq. (3) X = distance of a point within the soil -supported zone from the shoulder o f the span, X = |x| - L/2. X̅ = normalized value of X, X̅=X/x 1, where x 1 is giv en in Eq. (3) x1 = unit value for axial coordinate, and span length , see Eq. (3) xmode = distance from midspan beyond w hich the (transverse) mode shape displacements become negligible (x mode → ∞ for the ana lytical solution) y=y(t) = modal displacement ; factor applied to mode shape ϕ to get the dynamic displacement y(t) ϕ(x) α,β = dimensionless real parameters that define the roots of the characteristic equation in the outer (i.e. soil -supported) zone βH = parameter “ β” introduced by Hobbs [8], related to normalized span length by βH = L̅4 βF105 = parameter denoted by “ β” in [2], given by βF105 = log 10(βH) = 4 log 10(L̅) Δg = amount by which the pipe beco mes longer (if positive) or shorte r (if neg ative) due to the dynamic transverse displacements y(t)ϕ(x) ΔN = dynamic change in effective axial force arising from arch action for axially constrained case (Section 4) ϕ = mode shap e displacement, ϕ = ϕ(x), normalized to a maximum displacement of 1 λ = root of characteristic equation (in general complex) Π(.) = change in total potential energy due to a displacement increment (.) from the static equilibrium position Πu(.) = change in total potential energy for axially unconstrained conditions due to a displacement increment (.) from the static equilibrium position ξbending , ξN, ξsoil, ξarch = component damping ratio for bending, the axial force, the soil , and arch action ; the energy d issipated per cycle in the component is 4π times the component damping ratio times the maximum potential energy change of the component over the cycle from the equilibrium position ζ = modal damping ratio ζu = modal damping ratio for axially unconstrained conditions, ignoring the effect of arch action for the cross -flow mode ζbending , ζN, ζsoil = contribution to the modal damping ratio under axially unconstrained conditions from damping in bending, of the axial force, and in the soil ω = undamped natural f requency in radians per unit time : up to and including Section 3.5, arch action is not included , thereafter ω is re -defined as the natural frequency including the effect of arch acti on ω1 = unit value for angular natural frequen cy given in Eq. (3) ωu = undamped natural frequency for axially unconstrained conditions ωarch = frequency parameter associated with arch action, ω2 = ω u2 + ω arch2 (.),x = derivative of (.) with respect to x ; (.), x = (.), x̅ /x1 (.),x̅ = derivative of (.) with respect to x̅ ∫mode integral over x ∈ [-xmode, xmode] with x mode → ∞ ∫soil integral over the part of the pipeline that is supported by the soil and mode shape displacements are significant ; if (.) is symmetric in x, ∫ soil (.) dx = 2 ∫ [L/2,∞) (.) dx 5 3 Solu tion for Axially Unconstrained Case A single span is considered, as shown in Figure 1. The axial force N is assumed to be constant. For small amplitude vibrations, this is the case not only for an a xially unconstrained pipeline, but also for the in -line direction if the stati c lateral displacements are zero. For vertical vibrations, the arch action associated with sag under gravity loads together with axial constraint at the shoulders of the span is considered in Section 4. The results presented in this section apply for vertical (cross -flow) or horizontal (in -line) vibrations, provided the appropriate soil spring stiffness k is used for the direction considered. The the ory used applie s for slopes that remain small compared to 1, but as long as that is the case the axially unconstrained mode shap es and natural frequencies do not depend on the static equilibrium profile nor any influence the seabed geometry may have on it . However , it is limited to a single span with continuous contact in the soil supported zone and therefore does not apply for seabed geometries leading to multiple interacting spans. Figure 1 Span geometry and coordinates 3.1 Formula tion Classical (Euler -Bernoulli) beam theory is used, except that the geometrically nonlinear Von Kármán term for the axial strain is included. Via the principle of virtual work, this contributes a term involving the axial force to the transverse equilibr ium equation. The effective axial force N is assumed to be an input to the analysis , determined from a prior static analysis including the effect of lay tension and considering the loading history since installation of the pipeline. The elastic pipe has flexural rigidity EI and effective mass per unit length m . It is unsupported in the region x ∈ (-L/2, L/2) where L denotes the span length and elsewhere it is supported by a bed of linear elastic Winkler springs of stiffness k (force per unit length alon g pipe per unit deflection of the springs). The governing equations for the mode sha pe ϕ=ϕ(x) and natural angular frequency ω is EI ϕ,xxxx - N ϕ,xx – ω2 m ϕ = 0 for x ∈ (-L/2, L/2) (1) EI ϕ,xxxx - N ϕ,xx + (k – ω2 m) ϕ = 0 |x| ≥ L/2 (2) Where N denotes the effective axial force in the pipeline (positive for tension), and (.), x denote s differentiation with respect to the axial coordinate x. For simplicity, the added mass outside the span is assumed to be the same as inside. In practice, this is not an important limitation . One can always select a frequency -dependent stiffness k so that k – ω2 m has the desired value . Ideally this value is determined as the real part of a complex impeda nce function from a pipe -soil-water interaction analysis. (The authors are not aware of any such analyses reported in the literature, though [ 22-25] are example s of impedance estimation not including the water and using [24] or [ 25] to obtain impedance requires the pipeline to be approximated as a strip footing .) The same eq uations apply for the horizontal and vertical directions; only the soil s tiffness k is different, and this leads to different solutions ϕ and ω for the lateral and vertical directions. The above equation is also valid when the pipe is axially unconstraine d or when the static deflection in the direction of vibrations considered is zero. (In the latter case any changes in axial load are of 2nd order in the dynamic displacements and therefore lost in the linearization involved in any 6 modal analysis.) For the vertical case with sag and axial constraint , an approximate solution building on this one is given in Section 4. Returning to this axially unconstrained case, o ne seeks values for ω for which non -trivial solutions ϕ=ϕ(x) are possible that become vanishingly small a s |x| → ∞. Further , ϕ and its first 3 derivatives must be continuous at x = ±L/2. To make ϕ unique , there is also a normalizatio n condition . For the first symmetric mode , this is taken to be ϕ=1 at x=0 (midspan). 3.2 Normaliz ation In order to simplify the solution process, and reduce the number of problem input parameters, it is convenient to rewrite the above in normalized (non -dimen sional) form. For thi s purpose, x, N, and ω are written as unit values multiplied by normalized values. The unit values are given by x1 = (EI/k)1/4 , N 1 = (EI k)1/2 , ω1 = (k/m)1/2 (3) and the corresponding normalized values are denoted by an o verbar, as in x̅ = x/x 1 , N̅ = N/N 1 , f̅ = ω/ω 1 (4) The unit value x 1 also applies to the span length, so that L = L̅ x1. In terms of the normalized variables, the formulation of the eigenvalue problem from Eqs. (1) and (2) becomes ϕ,x̅ x̅ x̅ x̅ - N̅ ϕ,x̅ x̅ – f̅2 ϕ = 0 for x̅ ∈ (-L̅/2, L̅/2) (5) ϕ,x̅ x̅ x̅ x̅ - N̅ ϕ,x̅ x̅ + (1 – f̅2) ϕ = 0 |x̅| ≥ L̅/2 (6) The boundary conditions, whether normalized or not , are that ϕ and its first 3 derivatives must be continuous at x̅ = ±L̅/2. A normalized axial force of N̅ = -2 corresponds to the axial force at which the pipe buckles if it is continuous ly supported by the soil (with n o span). Thus, N > -2. Further , f̅ = 1 is the frequency at which the pipe vibrate s as a rigid body if it is continuously supported by the soil. Since an unsupported length reduces the lowest natural frequency, 0 < f̅ < 1 for the first mode. It is clear from Eqs. (5) and (6) that only N̅ and L̅ affect the normalized solution. The normalized span length L̅ =L(k/EI)1/4 is related to the soil stiffness parameter “β” used in [8,14,11] according to βH = L̅4 = L4 k/EI (7) where βH denotes the parameter “β” introduced by Hobbs [8]. The term “normalized” rathe r than “dimensionless” is used here, since non -dimensionali zation may imply the simplifications that can be achieved by application of Buckingham’s Pi Theorem, but the normalization applied here goes beyond that, taking advantage of simplifications that th e approximate formulation permits. Here all normalized quantities are also dimensionless, but the converse is not necessarily true. 3.3 Solution The zones x̅ ∈ (-L̅/2, L̅/2) and elsewhere will be referred to as the inner and outer zones, respectively. In general, the solutions can be written as exp( λx̅), where λ can be complex and different values apply for the inner and outer zones. For both zones, the characteristi c equation is quadratic in λ2, so the standard formula for the solution of a quadratic can be applied to get λ2 = N̅/2 ± Δ1/2 if Δ > 0 or λ2 = N̅/2 ± i ( -Δ) 1/2 if Δ < 0 (8) 7 where the values of Δ for the i nner and outer zones are Δinner = (N̅/2)2 + f̅2, Δouter = (N̅/2)2 + f̅2 – 1 (9) If the point ( N̅/2,f̅) falls within a circle of radius 1 centered at the origin, one has that Δinner > 0 and Δ outer < 0 (10) For this case, the roots of the characteristic equation are of the form λinner = ±a or ±ib and λ outer = ±α ± iβ where a, b, α and β are positive real numbers that can be calculated by evaluating th e expressions below involving real n umbers only in the order stated: a = (Δinner1/2 + N̅/2)1/2, b = ( Δinner1/2 – N̅/2)1/2 (11) r2 = (1-f̅2)1/2 , α = (( r2 + N̅/2)/2)1/2 , β = (( r2 – N̅/2)/2)1/2 (12) Solutions that are either sym metric or antisymmetric and decay away from the span can then be written as ϕ = A 1 f1 + A 2 f2 for |x̅| < L̅/2 (13) ϕ = A 3 f3 + A 4 f4 for |x̅| > L̅/2 (14) where, A1 to A 4 are arbitrary constant s of integration to be determined from the continuity and mode normalization conditions . For the symmetric case, f1 = cosh(a x̅), f 2=cos(b x̅), f 3=exp( -αX̅) sin( βX̅), f4=exp( -αX̅) cos( βX̅), (15) and for antis ymmetric case, f1 = sinh(a x̅), f 2=sin(b x̅), f 3=sgn( x̅) exp( -αX̅) cos( βX̅), f4=sgn( x̅) exp( -αX̅) sin( βX̅), (16) where X̅ = |x̅| - L̅/2, and sgn( x̅)=1 for x̅>0 and -1 for x̅<0. The boundary conditions are that ϕ and its first 3 derivatives mus t be continuous at x̅=L̅/2. (If they are continuous there , they will also be continuous at x̅=-L̅/2, by virtue of the way these solutions have been constructed as symmetric or anti -symmetric. If one side is continuous, so is its reflection.) These boundar y conditions give rise to 4 simultaneous linear, homogeneous equations in the constants A1 to A 4, in the form fi1 A1 + fi2 A2 – fi3 A3 – fi4 A4 = 0 for i=0,1,2,3 (17) where f ij denotes the ith derivative of f j with respect to x̅, evaluated at x̅=L̅/2. Non-trivial solutions are possible only when the determinant of the coefficients vanishes. This determinant is den oted by det( F). Where F denotes the 4 by 4 matrix whose rows are [fi1 fi2 -fi3 -fi4] for i=0,1,2,3. The natural frequency corresponds to the value of f̅ for which det( F) vanishes. This is calculated numerically. T he constants A i are then determined from the homogeneous equation together with the mode normalization condition , ϕ=1 at x̅=0. The bending stress at vibr ation amplitude of one pipe diameter, D, can then be written as σb1 = -ERD ϕ,x̅x̅ (k/EI)1/2 (18) where σb1 denotes the bending stress due to displacements Dϕ and is positive for tension on the side towards which the displacemen ts increase; E denotes Young’s modulus for the pipeline steel, or other material in which the bending stress is being calculated, D is the pipe diameter including coatings, R is 8 the radius to the location where the bending stress is calculated (e.g. half the outer diameter of the pipe, excluding the coatings) . This bending stress does not include the possible stress that can arise due a stiff coating (like concret e) that is discontinuous at field joints. (For that one would need to apply a stress concent ration fac tor equal to EI/EI FJ where EI FJ is the flexural rigidity at the field joint where the concentrated stress is being determined. The underlying assumption for dealing with a local change in flexural rigidity in this manner is that the field joint is short enough that its lower flexural rigidity does not sign ificantly affect the mode shape and frequencies, but it does affect the local curvature at the field joint.) 3.4 Range of Validity of the Above Solution In compression, the axial force cannot be la rger than the buckling force. This buckling force has been given by Het ényi [15] and can also be obtained from the above solution by setting the frequency f̅ to zero. The lower line in Figure 2 shows t his buckling load , |N̅|=N̅c. (Both the Het ényi equation [15] and setting f̅=0 in the above for mulation yield the same value to within the numerical truncation errors.) The validity condition ( N̅/2)2 + f̅2 < 1 is satisfied at a nd above the buckling load in Figure 2, until the upper line in Figure 2, where (N̅/2)2 + f̅2 = 1. This line will be referred to as the high-tension limit. The solution at and above this high -tension limit is given in Appendix A, as well as in Sollund et al. [11]. However , this high -tension solution still requires that f̅ < 1. This condition is satisfied for the lowest frequency mode, but not necessarily for higher modes. Soluti ons for f̅ > 1 can also be obtained, but these involve standing waves extending to infinity at constant amplitude. Thus , normal modes are not realistic at such frequencies. Inst ead, energy is radiated away from the span in the form of elastic waves . Such pipe radiation damping can thus develop at frequencies f̅ = ω (m/k)1/2 > 1 in addition to soil radiation damping that can be represented as damping associated with the soil spri ngs [22]. (There can also be soil radiation dam ping, but this is included as the damping of the soil springs as described in Section 6, and not in the eigenvalue problem for the undamped modes considered here.) Figure 2 Range of validity of standar d solution presented in Section 3.3, between the high -tension limi t above, where (N̅/2)2 + f̅2 = 1, and buckling force below. See Appendix A for solution at and above the high -tension limit. 9 3.5 Parametric Study The solutions for the axially unconstrained case are shown in Figure 4 to Figure 6. These plots show the ratio of the solution parameter plotted divided by the corresponding value for a simply supported span at the same N/N c ratio. (Since N c differs from the value for a si mply supported span , as shown in Figure 3, the same N/N c ratio means different N values.) The simply supported span solution is shown in Table 1. Observations and remarks in regard to Figure 3 to Figure 6 are as foll ows: 1) Figure 3 to Figure 6 together with Table 1 provide a simple means to obtain key solution parameters of interest. However , for applications, it is more convenient and accurate to use a look - up table, which is provided in [26] together with a simple look -up and interpolation algorithm. An approximate response surface is also provided in [ 11]. 2) Calculating numerically the solu tion to Det( F)=0 requires a suitable initial guess . Otherwise, the solution may converge to a higher mode, or the iterations could go outside the domain of v alidity of the algorithm for Det( F) used . For instance, this can happen due to a transition from the standard solution of Section 3.3, to the high -tension solution of Appendix A. The results presented in Figure 4 to Figure 6 include results in the high -tension zone. 3) Large L̅ = L (k/EI)1/4 occurs for large soil stiff ness k . The buckling load then approaches the fixed -fixed value of 4 times the pinned -pinned value, as expected. The pinned -pinned value is matched at L̅ = π. Below that , the buckling load is lower than under p inned -pinned conditions. 4) For high tension compared to the buckling load and long spans, the span behaves more like a cable than a beam. Under such conditions: a) The pipe develops rotation over a short hin ge zone near the shoulder of the span . Here the curvature and associated bending stresses become very high compared to other parts of the span. b) The flexural stiffness has little effect on the natural frequency . Therefore , the natural frequency converges to the same value as for a simply supported span at t he same axial force, as illustrated in Figure 7. Such convergence is not seen in Figure 4, because for a given L, Figure 4 shows the ratio of frequencies at the same N/N c but diffe rent N. Descript ion Formula in normalized values Formula in actual values Natural Frequency f̅ = (π/L̅)2 (1 + N̅/N̅c)1/2 ω = (π/L)2 (EI/m)1/2 (1 + N/N c)1/2 Bucking Load N̅c = (π/L̅)2 Nc = EI (π/L)2 Maximum curvature associated with mode normalized to ϕ = 1 at midspan |ϕ,x̅x̅| = (π/L̅)2 at midspan |ϕ,xx| = (π/L)2 at midspan Table 1 Solution for simply -supported (i.e. pinned -pinned) span 10 Figure 3 Buckling load divided by value for simply -supported (i.e. pinned -pinned) pipe with the same span length and flexural rigidity. Figure 4 Natural frequency divided by natural frequency of a simply supported span of the same length and buckling load ratio N/N c. (Note that the same N/N c in general means different N values. See Figure 7 for a plot of frequ encies for the same N values in tension.) 11 Figure 5 Maximum Curvature or bending strain divided by that for a simply -supported span of the same length L and buckling load ratio N/N c Figure 6 Mode shape displacement ϕ at the shoulders of the span (for mode normalized to ϕ=1 at midspan) 12 Figure 7 Natural frequencies divided by value for simply supported span at the same axial force N. 3.6 Comparison with F105 Approximation Here “F105 Approximation” refers to the formulation given in Sections 6.8 and 6.9 of [2]. These come from Fyrileif and M ørk [14] when 3.16 ≤ L̅ ≤ 100 (2 ≤ βF105 ≤ 8) and from [ 27] when 0.18 ≤ L̅ ≤ 3.16 (-3 ≤ βF105 ≤ 2), respectively . The ratio of the exact solution divided by this F105 approximation is shown in Figure 8 and Figure 9. This applies to the in -line case, i.e. when arch action is negligible . The plot s only cover cases within the stated range of validity of the F105 approximation where these limits can be expressed in terms of the normalized parameters N̅ and L̅ that fully define the normalized problem. (The flexural rig idity including concrete coating contribution if any is denoted by EI here but written as “(1+CSF) EI” in [2]. However, where the validity limit in [ 2] is in terms of “EI” instead of “(1+CSF) EI”, “( 1+CSF)” is assumed to be 1 in enforcing the validity limit for the plot s. This leaves “L/D s ≤ 140” as the only validity limit in [2] that cannot be expressed in terms of the normalized problem parameters N̅ and L̅.) It is seen that the F105 approximation is good, especially for “very short” ( L̅ < 3.16) spans. For longer spans with larger N/N c values differences appear . However very large N/N c values can occur only for very long spans where the pipe hangs like a cable . This cannot occur for spans within the “L/D s < 140” validity limit of [2] without excessive axial stresses . Since N c is also approximated in [ 2] as shown in Figure 3, it does make a small difference whether the results are plotted for the same N, or for the same N/N c. In Figure 8 and Figure 9, they are plotted for the same N/N c value where L̅ > 3.16 and for the same N where L̅ < 3.16. 13 Figure 8 Natural frequencies divided by those from the F105 appr oximati on Figure 9 Maximum curvature or bending stress divided by the F105 approximation 4 Rayleigh -Ritz Approximation for Arch Action An arch relies on its shape and horizontal restraint at the abutments to carry vertical loads . Similar ly, a pipeline span with significant sag can act as an inverted arch if there is horizontal (i.e. axial) restraint at the shoulders of the span. The axial force remains constant when (a) deflections are only moderately large, (b) axial inertia is n egligible and (c) the pipeline is axially unconstrained. Here an approximate solution is developed in which the last of these assumptions is removed. It is based on a Rayleigh -Ritz approximation : The transverse displacements associated with the mode sha pe are assumed to be the same as for the axially unconstrained case. One simply includes the axial strain energy arising from the axial constraint in the Rayleigh quotient to obtain the natural frequency. Here a simple derivation of the result is presente d. A more rig orous one that also yields the damping associated with arch action in given in Appendix E. 14 Without the axial constraint at the shoulders, the problem is already linear . Thus , linearization is not required to define the eige nvalue problem for the mode shapes and frequencies , and the results do not depend on the static deflections . However, adding axial constraint makes the problem nonlinear . To obtain mode shap es and natural freque ncies, it is linearized about an equilibriu m state under the static loads . It i s therefore necessary to first determine the static deflection w=w(x). This is done in Appendix B for the same assumptions as for the dynamic case : the pipe is continuously supported by a flat seabed of linear springs except where it is spanning. However , w(x) may also be determined by other means, accounting for the actual seabed geometry. Even survey data for w(x) can be used directly in the formulation that follows. In this Section , all variables used refer to actu al (i.e. dimension al) values as opposed to normalized values, unless otherwise noted. The total (static + dynamic) transverse displacements assumed to be w(x)+y(t) ϕ(x), where ϕ(x) is the mode shape calculated for the axially uncon strained case. For give n transverse displ acements, neglecting axial inertia, the axial force and displacement quasi -static analysis yields an axial force that change s as a function of time , but it is still constant over the length of the span. Within the span this axial force can then be written as N + ΔN , where N now denotes the axial force in the static equilibrium configuration and ΔN denotes the change in axial force due to the dynamic transverse displacements y(t) ϕ(x). It is assumed that this quasi -static analysis has established a relatio nship between the amount Δg by which the pipe becomes longer due to dynamic lateral displacements and the corresponding change in effective axial force in the span ΔN. As in any modal analysis, this relationship must be linearized and is written in the fo rm ΔN = (EA/L ea) Δg (19) where L ea > L represents an effective length for axial constraint. A few examples can help with the estimation of this effective length Lea. A. If the soil provides no axial restraint at all until a dista nce Lu from the shoulders, and at that point the pipe is fully constrained axially, L ea = L + 2Lu. B. For axial restraint from distributed axial soil springs of stiffness k a, the dynamic axial force dies out exponentially with distance from the shoulder, with a characteristic length (EA/k a)1/2, leading to an effective length in regard to axial constraint of L ea = L + 2 (EA/k a)1/2. C. If frictional behavior is considered for the axial pipe -soil interaction, the ΔN -Δg relation becomes nonlinear and hysteretic . This case is not addressed here. ( To apply modal analysis , it is then necessary to linearize the nonlinear ΔN-Δg relationship , leading to amplitude - dependent stiffness and damping eff ects requiring an iterative solution process.) Based on the kinematic approximations for moderately large deflections, t he change in pipeline length due to the dynamic transverse displacement y(t)ϕ(x) is Δg = ∫ mode ½ {[ w .x + y ϕ.x]2 – w.x 2} dx (20) where ∫ mode includes not just the span, but also the supported portion of the pipe where the modal displacements are sign ificant . Using the analytical solutions for both ϕ(x) and w(x) , this integral is evaluated analyticall y over x ∈ (-∞,∞) , or, exploiting symmetry , as twice the integral over x ∈ (0,∞). By way of linearization , the higher order term in y is neglected, yielding Δg = y I w’ϕ’ where I w’ϕ’ = ∫mode w.x ϕ.x dx (21) For undamped free vibra tions y(t) = A sin(ωt), the maxim um value of y is A, and the maximum value of ẏ=dy/dt is ωA. When y is a maximum, ẏ and the kinetic energy are zero. Conversely, when ẏ is a maximum , y and the potential energy are zero. More precisely, these statement s apply to the change in potential energy from the equilibrium position . Thus , the potential energy when y is a maximu m must be equal to the kinetic energy when ẏ is a maximum. This is how Rayleigh’s quotient can be derived. Applying it first for the a xially unconstrained case and then for the axially constrained case leads to Πu(Aϕ) = ½ ( ωuA)2 m I ϕϕ where I ϕϕ = ∫mode ϕ2 dx (22) 15 Πu(Aϕ) + ½ (EA/L ea) Δg2 = ½ (ωA)2 m I ϕϕ (23) where Πu(Aϕ) represents the change in potential energy due to the displacements Aϕ under axially unconstrained conditions ; ωu represents the axially unconstrained natural frequency (denoted by ω in Sections 3.3 to 3.5), and ω now represents the natural frequency that accounts for the effect of axial constraint. Eliminating Π u(Aϕ) from the above two equations, substituting Δg = A I w’ϕ’ and solving for ω yields ω = (ω u2 + ω arch2)1/2 (24) where ωarch2 = EA Iw’ϕ’2/(L ea m I ϕϕ) (25) The same analysis may be applied for the in -line (i.e. horizontal) direction, but since the static deflection is zero for that direction (unless the sp an is part o f a route bend, for ins tance), there is no increase in the in-line frequency for arching action in the horizontal direction . Arching also has no effect on an antisymmetric mode, if the static deflection is symmetric. If the surveyed profile of the pipe provides accurate slopes, these can be used directly for the static slope w,x=w,x(x) in the expression for I w’ϕ’. In practice the survey data may contain some noise due to measurement errors. Differentiating to get slope tends to magnify this. Therefore , it is better to integrate the expression for Iw’ϕ’ in Eq. (21) by parts to obtain Iw’ϕ’ = – ∫mode w ϕ. xx dx (26) This is in a more suitable form to evaluate the integral numerically with survey data for w=w(x). In Eq. (26) it is envisioned that the integral is taken over a long enou gh length so that the modeshape displacement ϕ and its derivatives are negligible at the ends, and therefore it is not necessary to incl ude the boundary terms that arise from the integration by parts. It can readily be verified that adding c 1 + c 2 x to w , where c 1 and c 2 are constants , will not affect the integral in Eq. (26) or (21) if evaluated exactly, but it could affect it slightly if evaluated numerically. To minimize tha t effect, especially when w ma y be large because it is measured from the mean sea level, it is recommended to remove the mean value of w over the range of the integral before calculating the integral in Eq. (26) numerically . 5 Numerical Example In this section the application of the analytical solution is illustrated for a specific case and the results compared to those from finite element (FE) analysis , for the axially c onstrained as well as unconstrained cases. The example is chosen so that the sag has an impo rtant, but not dominant influence on the first cross -flow mode. For this a L=61m long span, of a 610mm (24”) diameter pipeline with a wall thickness of 17.1mm (0. 7”) is considered . Including corrosion and weight coatings the diameter increases to D c=702m m (27.64”). The pipe is assumed to remain elastic with a flexural stiffness of EI=345 MNm2 and an axial stiffness of EA=7844.86MN. (This also accounts for the stiffness of the concrete weight coating in the approximate, linearized way of [2].) During the operating conditions considered for the VIV assessment, the submerged weight of the pipe is q=880.9N/m, the effective mass per unit length is m=883.24kg/m, and the effective axial force is N= -1.5MN (i.e. in compression). Th e dynamic soil stiffness is assumed to be k L= 18 MN/m2 in the lateral direction, and k V = 28 MN/m2 in the vertical direction. The static soil stiffness is ks = 0.53 MN/m2 for the vertical direction. For simplicity, i t is assumed that the pipe is fully unrestrained axially a distance of one span length beyon d the shoulder on each side, and then is fully fixed axially , so that L ea=3L. (This corresponds to the same level of axial restraint as one would get with axial distributed soil springs of stiffness ka=7.42MN/m2.) Calculating the normalized values for the span length L and effective axial force N yields: 16 L̅s = 12.08, L̅L =29.15, L̅V = 32.56 N̅s = -0.111, N̅L =0.0190, N̅V = -0.0153 where the subscript “s”, “L”, or “V” indicates normalization is for the static (vertical), dynamic lateral , or dynamic ve rtical case, respectively. Solving the dimensionless problems, and converting the results back to dimensional values yields the vertical displacement and cross -flow mode shape (Figure 10), bending stress (Figure 11) and natural frequencies4 of f=f u=0.3563 Hz for the first in -line mode, and fu=0.3663 Hz for the first cross -flow mode if axially uncon strained . With the axial constraint, t he vertical frequen cy increases to f=0. 4506 Hz. Clearly at a midspan static deflection of 0. 813 D , the sag together with the axial constraint already have a significant influence on the natural frequency of the first cross -flow mode. A finite element model is constructed wi th Abaqus, version 6.11 [17]. The pipe is modeled with beam elements of type “B31H ”. This model includes the full geometric nonlinearities under the static loading rather than only the Von Karman nonlinear term in the axial st rain displacement relations. Further, the linearization for the eigenvalue analysis is performed for this fully nonlinear model . Thus , the finite element model does not rely on the small slopes approximation of the analytical model. However, in both cas es small vibrations are considered to justify the linear ization involved in the modal analysis. Based on the dynamic vertical soil spring stiffness, the characteristic length in the soil -supported zone is x 1=1.87m. To obtain good accuracy from the finite element solution, the element size must be small compare d to this characteristic length. Hence an element length of 0.25m is chosen. For simplicity, the same element length is used throughout the model, although it would suffice to use small elements in the soil supported zone up to 2 characteristic lengths f rom the shoulder. The soil is modeled by a discrete spring at each node, with the stiffness of each discrete spring being calculated by multiplying the distributed spring stiffness k by the element l ength, except that half the element length is used at shoulders and at the end of the model . The FE modeling procedure is as follows: 1) Start with soil springs active in the vertical ( with stiffness set to the static value ) and lateral (with dynamic stiffnes s) directions . The pipe is axially restrained only at the left side of the model (i.e. at x = -91.5m). Also , all 3 rotation components are restrained at both ends. 2) Apply the effective axial load at the right side of the mode l (i.e. at x= +91.5m), and sim ultaneously apply the gravity load q, which is assumed to be constant over the whole length of the pipe. 3) Add additional vertical springs (strain free) with stiffness equal to the dynamic vertical stiffness minus the static vertical stiffness. The combinati on of vertical springs provides the dynamic stiffness w ithout altering the pipe sag across the span. 4) Perform axially unconstrained modal analysis. 5) Fix the node at the right end of the model axially. 6) Perform axially constrained modal analysis. The axially c onstrained FEA results are compared with analytical res ults in Figure 10 and Figure 11, showing vertical displacement and bending momen t, respectively. The static results from the FEA show excellent agreeme nt with the analytic result . The dynamic FEA result demonstrates slight differences from the analytic result (seen in Figure 11) at the locations of maximum curvature; but overall also shows very good agreement. The axially unconst rained FEA results (not shown) are visually indistinguishable from the analytic results. 4 Natural frequencies have been converted from the circular frequencies ω in radians/s to frequencies f=ω/(2π) in Hz. 17 Figure 10 Analytical solution and FEA results for “static” deflections under gravity load, and “dynamic” cross -flow mode, when mid -span moda l displacement is normalized to mid -span static displacement Figure 11 Analytical solution and FEA results for static bending moment s under gravity loads, and dynamic bending moment associated with the first cross -flow mode when mid-span modal displacement is normalized to mid-span static displacement The effect of the static deflection at midspan is investigated by repeating the FEA for different effective pipe weights. The results are shown in Figure 12 for the effective length Lea of 3 times the span length, and in Figure 13 for an effective length of 5 span lengths. As can be seen , the Rayleigh Ritz approximation provides a good approximation to the FE results (shown as crosse s), and also agrees with the DNV res ult, when the effective axial restraint span length is 3 times the span length . However, DNV does not account for variation in frequency resulting from variation in axial restraint at the span shoulders. 18 Figure 12 Results of FE a nalysis with effective leng th three times the actual span length . Figure 13 Results of FE analysis with effective length five times the actual span length. 6 Damping The concept of modal damping applies only to classical modes. I. e. it applies only when the mode shape s determined for the undamped system are ort hogonal with respect to the damping matrix . In particular , it applies for Rayleigh or Caughey damping matrices. Even if this damping orthogonality does not hold, one can assume that it holds , as an approximation. Although not stated as such, t his is the approximation advocated in Section 4.6 of [28]. For single -mode response, this corresponds to a Rayleigh -Ritz approximation, but, if several modes are responding (as occurs not just in very long spans, but also in multiple spans that are sufficiently close to each other to interact), it can no longer be justified as a Rayleigh Ritz a pproximation. Unless otherwis e stated , the results in this section apply for classical modes, and under axially unconstrained conditions. The generalization to include effect s of arch action is given in Appendix F . The objective is to estimate equivalent viscous damping to match the energy dissipation per vibration cycle. The relationship between the modal damping ratio ζ and the energy dissipated per cycle E d is then: ζ = E d/(4π U total) (27) 19 where U total denotes maximum increase in total potential energy relative to the static equilibrium position. The energy dissipated, E d, includes 3 components: Ed,bending from bending deformation, Ed,N from axial deformation associated with the effective axial force N, and E d,soil from the soil. The potential energy change from the equilibrium position for each of the components may be written as ½ K bending A2, ½ K N A2, and ½ K soil A2, respectively, where Kbending = EI Iϕ’’ϕ’’, K N = N I ϕ’ϕ’ , K soil = k Iϕϕ (28) where Iϕ’’ϕ’’ = ∫mode ϕ,xx2 dx, I ϕ’ϕ’ = ∫mode ϕ,x2 dx, I ϕϕ = ∫mode ϕ2 dx, (29) The K values can be interpreted as contributions to the modal s tiffness. They are all positive, except for K N, which has the same sign as N. For each of the components, the energy dissipated per cycle is written as a fraction 4πξ of the change in potential energy of that component given above. Thus , the total energy dissipated can be written as Ed = ½ A2 4π (ξ bending Kbending + ξ N KN + ξ soil Ksoil) (30) where ξ bending , ξN and ξ soil are the component damping ratios for the pipe bending deformations, the axial force and the soil springs, respectively . The total potential energy can be written as ½ K tu A2 where Ktu = K bending + K N + K soil = ω2 m I ϕϕ (31) Substituting into Eq. (27) then leads to the following expression for the modal damping ratio ζ = [K bending ξbending + KN ξN + Ksoil ξsoil ]/Ktu (32) Here the subscript “u” in K tu is a re minder that the case considered here is either axially unconstrained, or the static deflections are zero. All K value s can be interpreted as a modal stiffness, defined such that the corresponding change potential energy due to dynamic transverse displacem ents yϕ from the static equilibrium position is ½ K y2. It is important to distinguish between component damping ratio s ξ and the contribution ζ of these components to modal damping. For the axially unconstrained case, the latter are given by ζbending = ξ bending Kbending /Ktu, ζN = ξ N KN/Ktu, ζ soil = ξ soil Ksoil/Ktu (33) so that ζ = ζ bending + ζ N + ζ soil (34) The ξ values are properties of each component whereas the ζ values depend on the system as a whole. The estimates given in [28] are ζ values. The “structural damping ” ratio th erein is ζ bending + ζ N. Hydrodynamic damping is significant but has not been included as one of the “components” above, because it is treated separately, together with the hydrodynamic excitation from the water flowing around the pipe. As in [ 2], hydrodynamic damping is therefore not included in ζ. It is shown in Appendix C that ξN=0 whenever energy dissipation is by viscous damping in the pipe and the soil. This also applies if the length of the model is finite and a dashpot is applied to any of the degrees of freedom at the ends of the model. To facilitate the estimation of damping , the values of K bending /Ktu, K N/Ktu, and Ksoil/Ktu have been obtained by calculating the integrals analytically, as described in Appendix D . The results are shown in Figure 14 to Figure 16 as a function of the buckling load ratio N/N c and the normalized span length 20 L̅ = L(k/EI)1/4. As above “potential energy” always refers to the change in potential energy from the static equilibrium position due to the dynamic displacemen t, rather than the total potential energy. Figure 14 Ratio Kbending /Ktu of bending potential to total potential energy Figure 15 Ratio KN/Ktu of potential energy of the axial force N to total potential energy 21 Figure 16 Ratio K soil/Ktu of potential energy in the soil s prings to total potential energy In practice , a look -up table provides a quicker and more accurate way of getting the results in Figure 4 to Figure 7 and Figure 14 to Figure 16 for any particular case. Such a table together with a simple look-up and interpolation algorithm is provided in [26]. For the numerical example of Section 5, considering the in -line response , one has (L̅,N̅) = (29.15 , -0.01903 ). The normalized buckling force is N̅c =0.038574 , so that N/N c = -0.493 5. The potential energy fractions are Kbend ing/Ktotal = 1.823 , K N/Ktotal = -0.955 2 and K soil/Ktotal = 0.1321 . These values have been calculated directly, bu t they could also be obtained from Figure 14 to Figure 16. If the component damping ratios are ξ bending = 1% for bending, ξ N=0 for the axial force, and ξ soil = 12% for the soil, the overall damping is ζ = 1.803 ξbending -0.9208 ξN + 0.1178 ξsoil = 3.4%. Performing the same calculation for different values of the axial force N gives the results in Figure 17. Clearly, with reasonable assumptions for component damping ratios, modal d amping varies strongly with axial force becoming very large as the buckling load is approached . This is especially important to the assessment of in -line VIV , since that is quite sensitive to damping. 22 A similar calculation can be done for the cross -flow direction, using the vertical soil stiffness k=k V. However, this does not include the effect of arch action, which is important for the numerical example considered here. Appendix F provides the potential energy fraction for arch action. Arching energy involves soil as well as the pipe, in both cases for the axial direction. Estimation of damping where arch action plays an important role is described in Appendix F . The simple methods presented in this paper are mainly recommended for in-line VIV where there is no arching action (unless the span is within a horizontal bend) . This is also the case of most practical interest, because the influence of damping on in -line VIV is much greater than for cross -flow. Figure 17 Modal dam ping ratio for the numerical example of Section 5 for the in -line direction as a function of the axial force N 7 Conclusion s Subsea pipeline spans must be assessed not only as part of the design, but also as part of maintenance activities that can become necessary as a result of span creation by seabed scour. Considering this and the simplicity of the exact analytical solution p rovided in [11] and this paper , it is surprising that this analytical soluti on has not been provided earlier. Here that solution has also been applied to account approximately for arch action and estimate damping by breaking down the dynamic change in potential energy in different “components” of the system: essentially the pipe and the soil. For the axially unconstrained case (no arch action), t he results agree with the approximate solution in [2] and clarify its range of validity. However, when arch action is involved , the results in [ 2] only apply for a tacitly assumed level of axial restraint at the shoulder, whereas the Rayleigh Ritz approximation provided accounts for the axial stiffness at the shoulder by means of the equivalent axially unrestrained length Lea, which v aries between L ea = L (the span length) for axially fully restrained conditions at the shoulder of the span, and L ea → ∞ for axially fully unrestrained conditions. Based on finite element calculations for a numerical example the accuracy of the Rayleigh Ritz approximation for arch action is excellent, but there has been no systematic exploration of the range of validity of this approximation with respect to all relevant parameters affecting it, except that the gravity load and associated sag has been varied, providing an excellent approximation in all cases. As observed in [ 23], not only arch action, but also nonlinear effects are important under cross -flow VIV, since the amplitude of vibrations can easily reach or exceed the static deflection . However, nonlinear effects are not addressed here nor in modal analysis by finite -element [ 17] or Fourier mode s [5,11]. There remains a need for an at least approximate , practical, and general way to address such nonlinearity. Time -domain VIV assessment methods (see e.g. [29] and