Dynamics of a Pipeline Span on an Elastic Seabed

Natural frequencies, mode shapes and modal damping ratios must be estimated to assess subsea pipeline spans for vortex-induced vibrations (VIV) and response to direct wave loading. Several approximate solutions exist for a linearly elastic pipe under constant axial force supported by linearly elastic springs beyond the span’s shoulders. An exact analytical solution has only recently been published. That solution is used here in a Rayleigh-Ritz approximation to account for arch action arising from combined effects 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.


Introduction
Spans for subsea pipelines can arise from laying a pipeline on an uneven seabed, from current-induced scouring, wave action and/or other erosive phenomena [1] after the pipeline has been laid.This can lead to fatigue damage from vortex induced vibrations.Mostly this is assessed by modal analysis, together with empirical response functions, which DNVGL have periodically improved over the years as more experimental 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 pressure, temperature, pipelay), followed by the linearized eigenvalue analysis to obtain the natural frequencies and mode shapes [3][4][5].
Here an 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 frequencies, 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 been 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 concept of an 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.Further, 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 soilsupported 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, Vedeld 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 axial 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 accurate solutions a considerable number of Fourier modes need to be included giving 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 solution for this from Hetényi [15]) are much the same.On this result, they build approximate expressions 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 well 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 restraint 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 oscillations.(Cross-flow VIV amplitude can easily reach or exceed the static sag.)Both these effects have been considered by Bruschi and Vitali [16].A difficulty with nonlinear modal 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 largeamplitude 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][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][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 accounts for axial stiffness from soil support at the shoulders as well the static equilibrium profile of the span.The method is quite general.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 compared 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 component 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.

Solution 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 axially unconstrained pipeline, but also for the in-line direction if the static 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 theory used applies for slopes that remain small compared to 1, but as long as that is the case the axially unconstrained mode shapes 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.

Formulation
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 equilibrium 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 along pipe per unit deflection of the springs).The governing equations for the mode shape ϕ=ϕ(x) and natural angular frequency ω is Where N denotes the effective axial force in the pipeline (positive for tension), and (.),x denotes 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 impedance function from a pipe-soil-water interaction analysis.(The authors are not aware of any such analyses reported in the literature, though [22][23][24][25] are examples 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 equations apply for the horizontal and vertical directions; only the soil stiffness 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 unconstrained or when the static deflection in the direction of vibrations considered is zero.(In the latter case any changes in axial load are of 2 nd order in the dynamic displacements and therefore lost in the linearization involved in any 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, one seeks values for ω for which non-trivial solutions ϕ=ϕ(x) are possible that become vanishingly small as |x| → ∞.Further, ϕ and its first 3 derivatives must be continuous at x = ±L/2.To make ϕ unique, there is also a normalization condition.For the first symmetric mode, this is taken to be ϕ=1 at x=0 (midspan).

Normalization
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-dimensional) form.For this purpose, x, N, and ω are written as unit values multiplied by normalized values.The unit values are given by (3) and the corresponding normalized values are denoted by an overbar, as in The unit value x1 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 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 continuously supported by the soil (with no span).Thus, N > -2.
Further, f̅ = 1 is the frequency at which the pipe vibrates 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.
The term "normalized" rather than "dimensionless" is used here, since non-dimensionalization 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 the approximate formulation permits.Here all normalized quantities are also dimensionless, but the converse is not necessarily true.

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 characteristic equation is quadratic in λ 2 , so the standard formula for the solution of a quadratic can be applied to get where the values of Δ for the inner and outer zones are 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 the expressions below involving real numbers only in the order stated: Solutions that are either symmetric or antisymmetric and decay away from the span can then be written as where, A1 to A4 are arbitrary constants of integration to be determined from the continuity and mode normalization conditions.For the symmetric case, and for antisymmetric case, f1 = sinh(ax ̅ ), f2=sin(bx ̅ ), f3=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 must 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 boundary conditions give rise to 4 simultaneous linear, homogeneous equations in the constants A1 to A4, in the form where fij denotes the i th derivative of fj 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 denoted 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.The constants Ai are then determined from the homogeneous equation together with the mode normalization condition, ϕ=1 at x ̅ =0.
The bending stress at vibration amplitude of one pipe diameter, D, can then be written as where σb1 denotes the bending stress due to displacements Dϕ and is positive for tension on the side towards which the displacements 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 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 concrete) that is discontinuous at field joints.(For that one would need to apply a stress concentration factor equal to EI/EIFJ where EIFJ 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 significantly affect the mode shape and frequencies, but it does affect the local curvature at the field joint.)

Range of Validity of the Above Solution
In compression, the axial force cannot be larger 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 this buckling load, |N ̅ |=N ̅ c. (Both the Hetényi equation [15] and setting f̅ =0 in the above formulation yield the same value to within the numerical truncation errors.) The validity condition (N ̅ /2) 2 + f̅ 2 < 1 is satisfied at and 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.
Solutions 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.Instead, 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 springs [22].(There can also be soil radiation damping, 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 standard solution presented in Section 3.3, between the high-tension limit above, where (N ̅ /2) 2 + f̅ 2 = 1, and buckling force below.See Appendix A for solution at and above the high-tension limit.

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/Nc ratio.(Since Nc differs from the value for a simply supported span, as shown in Figure 3, the same N/Nc 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 follows: 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 lookup 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 solution 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 validity 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 stiffness 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 pinned-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 hinge 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 the 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/Nc but different N.
The plots 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 rigidity 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 plots.This leaves "L/Ds ≤ 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/Nc values differences appear.However very large N/Nc values can occur only for very long spans where the pipe hangs like a cable.This cannot occur for spans within the "L/Ds < 140" validity limit of [2] without excessive axial stresses.
Since Nc 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/Nc.In Figure 8 and Figure 9, they are plotted for the same N/Nc value where L̅ > 3.16 and for the same N where L̅ < 3.16.The transverse displacements associated with the mode shape 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 presented.A more rigorous one that also yields the damping associated with arch action in given in Appendix E.
Without the axial constraint at the shoulders, the problem is already linear.Thus, linearization is not required to define the eigenvalue 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 shapes and natural frequencies, it is linearized about an equilibrium state under the static loads.It is 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 actual (i.e.dimensional) 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 unconstrained case.For given transverse displacements, neglecting axial inertia, the axial force and displacement quasi-static analysis yields an axial force that changes 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 relationship 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 form ΔN = (EA/Lea) Δg (19) where Lea > 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 distance Lu from the shoulders, and at that point the pipe is fully constrained axially, Lea = L +2Lu.B. For axial restraint from distributed axial soil springs of stiffness ka, the dynamic axial force dies out exponentially with distance from the shoulder, with a characteristic length (EA/ka) 1/2 , leading to an effective length in regard to axial constraint of Lea = L + 2 (EA/ka) 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 amplitudedependent stiffness and damping effects requiring an iterative solution process.) Based on the kinematic approximations for moderately large deflections, the change in pipeline length due to the dynamic transverse displacement y(t)ϕ(x) is Δg = ∫mode ½ {[ w.x + y ϕ.x] 2w.x 2 } dx (20) where ∫mode includes not just the span, but also the supported portion of the pipe where the modal displacements are significant.Using the analytical solutions for both ϕ(x) and w(x), this integral is evaluated analytically 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 Iw'ϕ' where Iw'ϕ' = ∫mode w.x ϕ.x dx (21) For undamped free vibrations y(t) = A sin(ωt), the maximum value of y is A, and the maximum value of ẏ=dy/dt is ωA.When y is a maximum, ẏ and the kinetic energy are zero.Conversely, when ẏ is a maximum, y and the potential energy are zero.More precisely, these statements apply to the change in potential energy from the equilibrium position.Thus, the potential energy when y is a maximum must be equal to the kinetic energy when ẏ is a maximum.This is how Rayleigh's quotient can be derived.Applying it first for the axially unconstrained case and then for the axially constrained case leads to Πu(Aϕ) = ½ (ωuA) 2 m Iϕϕ where Iϕϕ = ∫mode ϕ 2 dx (22) Πu(Aϕ) + ½ (EA/Lea) Δg 2 = ½ (ω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 Iw'ϕ' and solving for ω yields ω = (ωu 2 + ωarch 2 ) 1/2 (24) where 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 span is part of a route bend, for instance), 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 Iw'ϕ'.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 enough length so that the modeshape displacement ϕ and its derivatives are negligible at the ends, and therefore it is not necessary to include the boundary terms that arise from the integration by parts.
It can readily be verified that adding c1 + c2 x to w, where c1 and c2 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 that effect, especially when w may 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.

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 constrained as well as unconstrained cases.The example is chosen so that the sag has an important, 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 Dc=702mm (27.64").The pipe is assumed to remain elastic with a flexural stiffness of EI=345 MNm 2 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).The dynamic soil stiffness is assumed to be kL= 18 MN/m 2 in the lateral direction, and kV = 28 MN/m 2 in the vertical direction.The static soil stiffness is ks = 0.53 MN/m 2 for the vertical direction.For simplicity, it is assumed that the pipe is fully unrestrained axially a distance of one span length beyond the shoulder on each side, and then is fully fixed axially, so that Lea=3L.(This corresponds to the same level of axial restraint as one would get with axial distributed soil springs of stiffness ka=7.42MN/m 2 .) Calculating the normalized values for the span length L and effective axial force N yields: L̅ s = 12.08, L̅ L =29.15, L̅ V = 32.56N ̅ s = -0.111,N ̅ L =0.0190, N ̅ V = -0.0153where the subscript "s", "L", or "V" indicates normalization is for the static (vertical), dynamic lateral, or dynamic vertical 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=fu=0.3563Hz for the first in-line mode, and fu=0.3663Hz for the first cross-flow mode if axially unconstrained.With the axial constraint, the vertical frequency 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 with 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 strain 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 cases small vibrations are considered to justify the linearization involved in the modal analysis.
Based on the dynamic vertical soil spring stiffness, the characteristic length in the soil-supported zone is x1=1.87m.To obtain good accuracy from the finite element solution, the element size must be small compared 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 from 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 length, 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 stiffness) 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 model (i.e. at x= +91.5m), and simultaneously 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 combination of vertical springs provides the dynamic stiffness without 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 constrained FEA results are compared with analytical results in Figure 10 and Figure 11, showing vertical displacement and bending moment, respectively.The static results from the FEA show excellent agreement 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 unconstrained FEA results (not shown) are visually indistinguishable from the analytic results.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 crosses), and also agrees with the DNV result, 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.

Figure 12
Results of FE analysis with effective length three times the actual span length.

Figure 13
Results of FE analysis with effective length five times the actual span length.

Damping
The concept of modal damping applies only to classical modes.I.e. it applies only when the mode shapes determined for the undamped system are orthogonal 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, this 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 approximation.
Unless otherwise stated, the results in this section apply for classical modes, and under axially unconstrained conditions.The generalization to include effects 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 Ed is then: where Utotal denotes maximum increase in total potential energy relative to the static equilibrium position.
The energy dissipated, Ed, includes 3 components: Ed,bending from bending deformation, Ed,N from axial deformation associated with the effective axial force N, and Ed,soil from the soil.The potential energy change from the equilibrium position for each of the components may be written as ½ Kbending A 2 , ½ KN A 2 , and ½ Ksoil A 2 , respectively, where Kbending = EI Iϕ''ϕ'', KN = N Iϕ'ϕ' , Ksoil = k Iϕϕ (28) where Iϕ''ϕ'' = ∫mode ϕ,xx 2 dx, Iϕ'ϕ' = ∫mode ϕ,x 2 dx, Iϕϕ = ∫mode ϕ 2 dx, The K values can be interpreted as contributions to the modal stiffness.They are all positive, except for KN, 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 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 ½ Ktu A 2 where Substituting into Eq.( 27) then leads to the following expression for the modal damping ratio ζ = [Kbending ξbending + KN ξN + Ksoil ξsoil ]/Ktu (32) Here the subscript "u" in Ktu is a reminder that the case considered here is either axially unconstrained, or the static deflections are zero.All K values can be interpreted as a modal stiffness, defined such that the corresponding change potential energy due to dynamic transverse displacements yϕ from the static equilibrium position is ½ K y 2 .
It is important to distinguish between component damping ratios ξ 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 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 therein 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 Kbending/Ktu, KN/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/Nc and the normalized span length 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 displacement, rather than the 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/Nc = -0.4935.The potential energy fractions are Kbending/Ktotal = 1.823,KN/Ktotal = -0.9552and Ksoil/Ktotal = 0.1321.These values have been calculated directly, but 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 damping 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.
A similar calculation can be done for the cross-flow direction, using the vertical soil stiffness k=kV.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 damping ratio for the numerical example of Section 5 for the in-line direction as a function of the axial force N

Conclusions
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 provided in [11] and this paper, it is surprising that this analytical solution 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), the 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 varies between Lea = L (the span length) for axially fully restrained conditions at the shoulder of the span, and Lea → ∞ 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 modes [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 references therein) do provide this possibility, but they involve much more computational effort without any guarantee that for the linear case the results agree with the established practice.
the static slope w' becomes negligible.Normalizing, proceeds as before, for the dynamic case, except that in this case the deflections w and axial elongation g also need to be normalized.For this the unit values are w1=q/k and g1=w1 2 /x1.In terms of normalized values and derivatives, the problem can be written as: Although the equations for unit values N1 and x1 are the same as for the dynamic case, i.e.Eq. ( 3), the effective soil stiffness k is generally much lower for the static case, due to nonlinear inelastic soil response for which the linear stiffnesses are used as an approximation.Therefore, in general the dimensionless span length, and the dimensionless axial force will be different for the dynamic and static solutions.
The symmetric solutions can be written as for N ̅ =0 = (cosh(|N ̅ | 1/2 x ̅ ) -1)/N ̅ for N>0 and f3 and f4 are calculated as for the symmetric mode in the dynamic case setting f̅ =0 for the calculation of the coefficients α and β.Indeed, the homogeneous part of the solution above is the same as the homogeneous solution for the mode shape equations when f̅ =0, except that in this case positive and zero values of N ̅ are considered, as well as negative values.Although different valid expressions could be picked for f2, the choice above has the advantage that all three expressions converge to f2 = ½ x ̅ 2 as |N ̅ | →0.
As before , the boundary conditions are that w ̅ and its first 3 derivatives must be continuous at x ̅ =L̅ /2.These boundary conditions give rise to 4 linear equations in the constants A1 to A4.In this case, the equations are not homogeneous (because of the particular solution) and can be solved for A1 to A4 unless the axial force is equal to the buckling load.

Appendix C -Conditions under which the axial force does not contribute to the energy dissipation
This appendix is self-contained in regard to notation.
Consider the damping matrix of a discretized system (e.g. by a finite element discretization) in which ui denote the nodal displacement and/or rotation components.The following assumptions are made: 1) The strain displacements relations are quadratic, so that the relationship between strains at any point within the system can be written as ϵp = Σi Bpi ui + ½ Σi Σj Gpij ui uj , where ϵp denotes the p th component of strain, and Bpi and Gpij are functions of location within that are derived from the displacement discretization and the strain-displacement relations, and Σ indicates summation over all values of the index appearing as a subscript.This applies for the theory of moderately large deflections, since there the strain displacements relations are linear, except for the expression for the axial strain which is quadratic.2) The materials involved are viscoelastic.I.e. it is assumed that the energy-conjugate stresses to the strains ϵp can be written as σp = Σq { Lpq ϵq + Dpq ϵq}, where Lpq and Dpq are material stiffness and viscous damping coefficients.For the system considered this applies when the pipe material, and the soil springs are viscoelastic.
Under these assumptions, applying standard methods to derive and linearize the governing equations leads to The model analyzed is assumed to be long enough so that the boundary conditions at the ends no longer influence the results.Those boundary conditions are assumed to be fully fixed for displacements and rotation, except that axial displacements at one of the two ends are allowed for the axially unconstrained case only.The integral ∫(.)dx is over this entire model.Virtual displacements are zero unless otherwise noted.
The theory of "moderately large displacements" means that the axial strain is approximated by ϵ = u' + ½ w' 2 , and the linearized expression for the curvature κ = -w'' is used.
Inertia is included as a transverse inertial force -m wd,tt.Axial inertia is neglected.with wd = ϕ sin(ωut), it is readily verified that this leads to the governing equations of Section 3.
(Recall that in Section 3, "ω" means the axially unconstrained natural frequency ωu.)Since those equations are satisfied exactly, the corresponding PVW equation is also satisfied exactly.I.e., ∫ (Ns ϕ'δw' + EI ϕ'' δw'' + kϕ δw) dx -∫ mωu 2 ϕ δw dx = 0 ∀ (δu,δw) (E.5) Return now to the dynamic portion of the PVW (E.3) that includes axial constraint, apply the Rayleigh-Ritz approximation wd=ϕy, and δw=ϕδy, and simplify using (E.4) to obtain ∫ (Nd (δu' + (ws'+ϕ'y) ϕ' δy) +rxd δu) dx -∫{qxd δu+(qzdmw,tt -ωu 2 mwd) ϕδy }dx=0 ∀ (δu,δy) (E.6) Setting δy=0, integrating by parts where necessary, applying the lemma of calculus of variations leads to Nd' + qxdrxd = 0. Further, Nd = EA ϵd, ϵd = ud' + ws' ϕ' y + ½ ϕ' 2 y 2 , rxd = ka ud.This means that Nd and ud are the axial force and displacement one would get for a perfectly straight pipe undergoing axial deformations only under a thermal contraction ws' ϕ' y + ½ ϕ' 2 y 2 .Now it is consistent with the linearization of all modal analysis to drop higher order terms in y.Further, there is no axial dynamic applied load, i.e. qxd=0.Thus, one obtains Neglecting the change in Ndy where the transverse dynamic displacements are still significant, the expression for Karch simplifies to Karch = Ndy,span Iw'ϕ', where Ndy,span denotes the constant value of Ndy within the span.This leads to free vibrations at a natural frequency of ω = (Kt/M) 1/2 .With these approximations the results in this appendix are seen to be the same as those in Section 4 and Appendix F, with a modal stiffness for arch action of Karch = (EA/Lea) Iw'ϕ' 2 (E.12) in both cases.Thus, this appendix serves to confirm the results, clarify the assumptions and approximations made and it provides the consistent way to calculate Lea.

Appendix F -Damping Associated with Arch Action
The appendix provides the potential energy fraction associated with arch action, and describes how the approach of Section 6 can be extended to account for arch action.The nomenclature of Section 2 applies here as well.A derivation of the same results that highlights the assumptions and approximations from a different perspective is given in Appendix E. An application of the methods to estimate damping associated with arch action for the numerical example of Section 5 is also included in this appendix.
For undamped vibrations y = A sin(ωt), the maximum speed is |ẏ|max = ωA.If the frequency increases from ωu to ω due to arch action, this means that the kinetic energy increases by a factor (ω/ωu) 2 , because the kinetic energy is proportional to velocity squared.What is true for the maximum kinetic energy must also be true for the maximum change in potential energy from the equilibrium position.Thus, the potential energy also increases a factor (ω/ωu) 2 due to arch action.Thus, the maximum change in total potential energy from the equilibrium position is of modal displacement y(t); corresponds to amplitude of vibrations at midspan for the first mode a,b = dimensionless real parameters that define roots of the characteristic equation in the inner zone D = outer diameter of pipe including 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 damping EI = effective flexural rigidity including contribution of pipe coatings, if applicable EIFJ = flexural rigidity at the field joints (not including the contribution of concrete coating that is not present there) f = natural frequency f̅ = normalized natural 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; f1 and f2 apply within the inner zone |x| ≤ L/2, f3 and f4 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 vertical loading respectively Ka, Kbending, KN, Ksoil = 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 ½ Kbending y 2 .Ktu = total modal stiffness under axially unconstrained conditions, Ktu = Kbending + KN + Ksoil Kt = total modal stiffness, including the effect of arch action arising from the combination of static deflection and axial restraint at the shoulders of the span, Kt = Ktu + Ka L = span length L̅ = normalized span length, L̅ = L/x1 Lea = effective length in regard to axial restraint: chosen such that a length Lea 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/N1 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 = Nc/N1, shown in Figure 2 and Figure 3 R = radius to point at which the bending stress is determined t = time w = transverse static equilibrium displacement; typically refers to the vertical displacement but the same formulation can be used for the lateral direction if the pipe has lateral out-ofstraightness, e.g.due to a route bend x = coordinate defining location on the pipeline; zero at midspan x ̅ = normalized coordinate, x ̅ = x/x1 where unit value x1 is given in Eq. (3) X = distance of a point within the soil-supported zone from the shoulder of the span, X = |x| -L/2.X ̅ = normalized value of X, X ̅ =X/x1, where x1 is given in Eq. (3) x1 = unit value for axial coordinate, and span length, see Eq. (3) xmode = distance from midspan beyond which the (transverse) mode shape displacements become negligible (xmode → ∞ for the analytical 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 = log10(βH) = 4 log10(L̅ ) Δg = amount by which the pipe becomes longer (if positive) or shorter (if negative) 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 shape 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 dissipated 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 frequency 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 action ω1 = unit value for angular natural frequency given in Eq. (3) ωu = undamped natural frequency for axially unconstrained conditions ωarch = frequency parameter associated with arch action, ω 2 = ωu 2 + ωarch 2

Figure 3
Figure 3Buckling load divided by value for simply-supported (i.e.pinned-pinned) pipe with the same span length and flexural rigidity.

Figure 4 Figure 5 Figure 6 Figure 7
Figure 4 Natural frequency divided by natural frequency of a simply supported span of the same length and buckling load ratio N/Nc.(Note that the same N/Nc in general means different N values.See Figure 7 for a plot of frequencies for the same N values in tension.)

Figure 8 Figure 9
Figure 8Natural frequencies divided by those from the F105 approximation

Figure 10
Figure 10Analytical solution and FEA results for "static" deflections under gravity load, and "dynamic" cross-flow mode, when mid-span modal displacement is normalized to mid-span static displacement

Table 1 Solution for simply-supported (i.e. pinned-pinned) span
denotes total axial force, instead of static portion only which now is denoted by Ns