Vortex-Induced Vibration Response Bifurcation Analysis of Top-Tensioned Riser Based on the Model of Variable Lift Coefficient

Abstract

Top-tensioned riser (TTR) is one of the most frequently used pieces of equipment in offshore petroleum engineering. At present, the research of its vortex-induced vibration (VIV) is mainly focused on numerical simulation with few analytical approaches. In this paper, the nonlinear dynamic equation of VIV about TTR is developed by introducing the third-order variable lift coefficient hydrodynamic model. The amplitude-frequency response equation under 1:1 primary resonance excitation is deduced based on the single mode discrete results. The nonlinear dynamic responses characteristics of the TTR vibration under the excitation of vortex-induced resonance are discussed. Bifurcation theory is used to study the singularity of the analytical solution of the riser system. The hysteresis and solitary solutions from the results of singularity analysis are found. Such results can provide guidance for the design and optimization of riser structural parameters.

1. Introduction

With the increasing demand for oil and gas resources globally, the exploitation of marine oil and gas resources gradually extends from shallow water to deep water. As the deep-sea risers are slender with large aspect ratio, the natural frequency and structure stability of the riser are reduced, and the vortex-induced vibration (VIV) of the deep-sea riser is more complicated with respect to the water depth. Under this trend, a lot of new problems and new phenomenon emerged [1]. At present, the dynamic studies of deep-sea riser are mainly concentrated in the field of VIV. VIV is essentially nonlinear, self-excited and self-limiting multi-degree of freedom resonance response, which generally involves strong nonlinear structural motions. Typical nonlinear characters such as lock-in, hysteresis, displacement jumps, bifurcation, and chaos can be found in VIV of deep-sea risers. Despite rich dynamics of the VIV of deep-sea risers, researchers have not been able to find exact analytical solutions of the VIV response. This hampers the understanding of the underlying mechanisms of some observed nonlinear phenomenon [2].

With the increasing number of VIV experiments and CFD simulations for deep-sea conditions in recent years, deeper understandings of the riser VIV are gained and thereafter lead to progress in accurate modeling of vortex lift. Sarpkaya measured the Fourier average of hydrodynamic force over many cycles of vibration. The lift force is decomposed into two parts, the drag part and the inertia part, which are, respectively, related to the velocity and acceleration of the vibrating cylinder. Sarpkaya pointed out that, for practical Reynolds numbers, the nonlinear expression with respect to structural motion could capture the hydrodynamic feature better than the linear expression [3]. Gopalkrishnan and Govardhan implemented plenty of VIV experiments and presented the lift coefficient in terms of structural motion [4, 5]. Vandiver suggested that a piecewise parabola function of structural amplitude could be used for industrial model of lift force to calculate riser displacement by using wake oscillator model [6]. Based on the experimental results of Gopalkrishnan and Vandiver et al., a third-order polynomial of the structure velocity was proposed to model the lift coefficient so as to take account of nonlinear interaction between structural and fluid dynamics [7], which provides a new model of lift coefficient for theoretical calculation and analysis of VIV about marine risers.

At present, there are two main methods to research VIV, namely, the model experiment method and numerical simulation method. In order to make a regular description of the riser movement, analyze the characteristics of the frequency domain, and express the influence of the parameters of the system on its motion, in this paper, the nonlinear vibration theory and the singularity theory are used to obtain the analytical solution and analyze the top-tensioned riser (TTR) system. The influence imposed by lift coefficient on the amplitude of the system is studied quantitatively. Such quantification can provide predictive and new ideas for the VIV problem about the TTR.

2. Dynamical Model and Method of Solution

For the deep-sea risers with large aspect ratio, we take the normal TTR as an example. It is assumed that the TTR is vertically supported by the platform. The body of TTR is a uniform circular cross section. The lower end of the riser is hinged on the universal joint, while the top end of the riser has top tension applied by the tensioner. Thus, TTR can be modeled as a simply supported beam with top tension. Assume the VIV caused by the current; the schematic diagram and simplified force diagram of the TTR are shown in Figures 1(a) and 1(b), respectively [8–10].

(a) VIV of a riser subject to current

(b) Simplified force diagram of riser

Figure 1 

Schematic diagram of vortex-induced vibration of the TTR system.

Based on the bending vibration theory of the Euler Bernoulli beam, the motion equation that governs the VIV at transverse flow direction (y direction) reads as follows: where EI is the bending stiffness; is the transverse displacement of the riser. T is the top tension, which can be expressed as , in which is the riser cross-sectional area, is the acceleration of gravity, and are the riser material density and seawater density, respectively, is the top tension coefficient, and is the length of the top tension riser [11]; is the riser mass per unit length, is the riser fluid mass per unit length, is the additional mass per unit length (, in which is the coefficient of the added mass), C is structural damping, D is the diameter of the riser, U is the velocity of the ocean current, C_L is the lift coefficient, and C_D is the drag coefficient. In this paper, the third-order variable lift coefficient model is used, which models the lift coefficient as [12], in which is vortex shedding frequency (as , and is the Strouhal number). The values of coefficients in the expression can be found from reference [7, 12], such as C_L0=0.5, C₁=1.82, C₂=-1.29, C₃=-0.71 or C_L0= 0.22, C₁= 1.62, C₂=-2.31, C₃= 0.75. Note that these coefficients C_L0, C₁, C₂, and C₃ are dimensionless.

Riser boundary conditions are as follows:

2.1. Dimensionless Equation

We introduce the following dimensionless quantities: , , and ( is time scale and has ). Let , by substituting the dimensionless transformation into (1), we get the dimensionless equation as follows:where (dimensionless frequency).

2.2. Galerkin Discretization

Let , is generalized coordinates, and let the th order mode function be the basis function. Note that the basis function satisfies the boundary condition specified in (2) and the orthogonality condition. This paper considers the first-order modal function, as N = 1. Multiply both side of (3) by the mode function and integrate both sides over the interval . Using the orthogonality of the mode function, the ordinary differential equation with the single mode vibration of the TTR can be expressed as [13]in which , , , , , .

2.3. Periodic Solution at Primary Resonance

In this paper, we solve the first-order approximation solution of the system by nonlinear vibration averaging method [14–16] when [17–19]. Assume the solution of (4) as follows:where with being the small quantity that characterizes the perturbation and .

Note the fact that ; by substituting this expression into (4) and (5), we have

In order to determine the amplitude-frequency response equations in steady state, let us set . This leads to the amplitude-frequency response equation in steady state shown in (7). The final form of (7) takes advantage of the trigonometric function relation .

With given dimensionless frequency , (7) enables the calculation of steady state amplitude via solving nonlinear algebraic equation, when is close to , which implies the occurrence of near resonance vibration. As a result, one can obtain the fixed point in state plane by solving (6) and (7) simultaneously.

2.4. Stability Analysis

Stability analysis of fixed points that belonged to the nonlinear system defined in (6) is carried out herein. The Jacobian matrix can be obtained from (6) with Y₁ and Z₁ representing the right-hand side of the differential equations. For the first-order approximation solution of the averaging method, , the characteristic equation can be written as follows:

Stability criterion of the system at each of its fixed points can be determined via its eigenvalues and computed from (8). The eigenvalues of stable solution should satisfy

3. Analytical Results and Discussion

At first, to validate the proposed riser dynamic model, we take the VIV experimental results from Ocean University of China as baseline [20]. The riser model studied in [20] is a 20 mm diameter copper pipe with a length of 6.2 meters. For validation purposes, the displacement envelope of the riser VIV model in [20] has been computed at flow velocity =0.85m/s for comparison.

The black dashed line in Figure 2 is the experimental data measured in [20]; the blue dash-dotted line is the calculated result from [20]; and the red line is the numerical result calculated by the proposed method in Section 2 of this paper. From Figure 2, it can be seen that the proposed method used in this paper results in closer dynamical response to the experimental data. The improved agreement between numerical and experimental results can lay the foundation for later numerical and theoretical analysis in this paper.

Figure 2 

Comparisons of the displacement envelopes between the presented numerical results and the existing experimental results of the riser model VIV at =0.85 m/s.

In this section, we take the 1500m length TTR with API X65 steel as the material under study [21]. The key parameters of the VIV governing equation are listed in Table 1.

Table 1 

Parameters of TTR.

3.1. Vibration Influence of Riser with External Excitation Frequency

The change of external excitation frequency is correlated with the change of flow velocity, which leads to the VIV response change. Figure 3 shows the bifurcation diagram of the TTR vortex-induced vibration response. Amplitude A varies with the external excitation frequency near the first-order natural frequency of structure. In Figure 3 the red dot line indicates the numerical solution and the black line is the analytical solution obtained from (7). It can be seen that when the vibration depicts stable periodic pattern, the analytical solution is consistent with the numerical solution.

Figure 3 

Bifurcation diagram of the TTR VIV response amplitude A varies with the external excitation frequency and the reduced velocity .

The reduced velocity is also an important reference for studying the VIV of the cylinder. In this section the explicit expression of the reduced velocity is , in which is the first-order natural frequency of the riser. The transverse vibration of VIV of the cylindrical pipeline occurs mainly in the range of reduced velocity from 3 to 7, and the range of 4 to 6 is of particular significance [22]. The upper axis in Figure 3 shows the relationship between the amplitude and the reduced velocity of the VIV. Compared with the results obtained in experimental study results, the variation range of VIV amplitude in the reduced velocity in this paper is relatively consistent. It is shown that the model used in this paper and its approximate solution can describe the VIV of the deep-sea TTR accurately, and it is feasible to study the influence of riser parameters based on the model.

From Figure 3 we can find that when the VIV occurs, the vertical cross-flow amplitude increases gradually with the increasing of the external excitation frequency (means that the flow velocity increases), and it will reach the maximum value (about 0.23 meters) of vibration near the natural frequency of the riser. Then the amplitude of vibration gradually reduces and remains basically unchanged in a certain range. When the external excitation frequency is close to the next-order natural frequency of the riser, the riser amplitude will significantly increase again with the increase of external excitation. In addition, along with the varying of external excitation frequency, the riser vibration response would undergo the process from almost periodic to periodic and then almost periodic again.

3.2. Parameter Study of Lift Coefficient

Note the lift coefficient in this study is modeled as ; for the calculation of rigid riser and flexible riser, the choice of parameters in C_L is different. Because the deep-sea riser aspect ratio is relatively large, with the length increases, the stiffness of riser system decreases and presents more characteristics of flexible riser. Sensitivity study indicates that the changes of parameters C₂ and C₃ have very little effects on C_L. The influences brought by C₂ and C₃ are orders of magnitude smaller than that of C_L0 and C₁. Therefore, in this paper, we only study the impact of varying C_L0 and C₁ on the riser nonlinear dynamics via singularity theory. C₂ and C₃ are taken as constants.

According to the singularity theory [23], the critical parameter condition qualitatively changes the topology of amplitude-frequency response; such condition defined in parameter space is called the transition set. The transition set includes the bifurcation set , the hysteresis set H, and the double limit set DL [24]. The bifurcation equation based on (7) can be expressed as where and are treated as the state variable and the bifurcation parameter, respectively; and are unfolding parameters. The solutions of the transition set can be, respectively, obtained from (11). Each transition set represents a typical topological pattern of the nonlinear dynamics; within each subregion of the parameter space, the dynamical system behaves similarly. In other words, (11) governs the shapes of subregion boundaries in (, ) plane that cateogize different solution tolopologies.

by replacing the dimensionless frequency ν and the dimensionless amplitude with vortex shedding frequency and amplitude and substituting them into (11), three types of transition set can be obtained in parameter space. It is found that the double limit set is empty; results of the other two transition sets are shown as follows.

From the above formula we can draw the boundaries of transition set as shown in Figure 4. In Figure 4, red solid line defined the boundary of the hysteresis set, and the blue dotted line is the bifurcation set. The transition set can divide C_L0 -C₁ plane into three regions (①, ②, ③).

Figure 4 

Transition set of parameters C_L0 and C₁ in lift coefficient.

In order to explore the difference of the amplitude-frequency response of the system in each region, we take three representative pairs of C_L0 and C₁ (0.06, 8), (0.015, 8), (0.01, 8) from regions ①, ②, and ③ in Figure 4 to calculate the amplitude-frequency response of the TTR. Analytical solutions under these parameter sets are compared by numerical calculations. Also, stability analyses are carried out within these subregions. Comparisons are displayed in Figure 5, in which the black solid line is the stable part of the analytical solution, the black dotted line is the unstable part, and the red dotted line is the numerical solution.

(a) Amplitude-frequency response of general single solution (region ①)

(b) Amplitude-frequency response of multiple solutions as hysteresis phenomenon (region ②)

(c) Amplitude-frequency response of isolated solution as bifurcation phenomenon (region ③)

Figure 5 

Amplitude-frequency response of the riser VIV calculated with parameters C_L0 and C₁ in different regions.

From Figure 5 we find that when C₁ stays unchanged, the amplitude will reduce with C_L0 gradually decreasing. In addition, with the unfolding parameters changing in different areas, the riser system has typical hysteresis phenomenon in region ②, and an isolated solution appears in region ③. This brings us to the conclusion that the system will have multiple solutions when parameters drop in regions ② and ③. In other words, this means that the same external excitation frequency may cause different values of amplitude. For large C₁/C_L0 ratio, which corresponds to region ③, the amplitude-frequency response may have multiple solutions. It can be found from the lift coefficient expression that C_L0 represents forced excitation due to vortex shedding and C₁ represents the viscous damping of the fluid. Thus the multiple solution case occurs when vortex shedding excitation force is small and the fluid damping is large.

In the case when the numerical solutions are periodic solutions, the stability part of the analytical solutions coincides quite well with the numerical solutions. As the other parts of analytical solutions are unstable in the frequency range of periodic solutions, they cannot be matched with the numerical solutions. When the numerical solutions are almost periodic solutions, the analytical solutions can only reflect the periodic solutions when the system is stable. Therefore, the analytical solutions in the almost periodic range cannot coincide with the almost periodic solutions.

3.3. Analysis of Transition Process

As shown in Figure 5, all three types of responses contain only one stable solution. In order to explore the difference of physical significance about the VIV between the cases of multiple solutions and the single solution, we study the transient process of the VIV from different initial conditions in this section. In order to calculate conveniently and seek the rule, we study the transient process of the VIV in the case of the isolated solution (one case of multiple solutions) and the general single solution with similar amplitude in steady state.

3.3.1. Evolution of Vector Field

Under dimensionless condition two parameters set with C_L0=0.01, C₁=8 and C_L0=0.06, C₁=2.1 are studied, respectively. Figure 6 shows the amplitude response under the first-order primary resonance. It can be seen that Figure 6(a) reveals the isolated solution in the amplitude-frequency response diagram, and Figure 6(b) has the general single solution. The maximum amplitudes of these two solutions are similar with the value around 0.265. Figures 6(c) and 6(d) display the vector fields with the fixed external excitation of the dimensionless frequency at =825. Drastic topological difference can be observed from Figures 6(c) and 6(d), which validates the prediction of parameter space boundary via singularity theory.

(a) Amplitude-frequency response of isolated solution (C_L0=0.01, C₁=8, and =825)

(b) Amplitude-frequency response of single solution (C_L0=0.06, C₁=2, and =825)

(c) Vector field diagram of isolated solution (C_L0=0.01, C₁=8, and =825)

(d) Vector field diagram of general single solution (C=0.06, C₁=2, and =825)

Figure 6 

Comparisons of amplitude-frequency response and vector field diagram between the isolated solution and single solution at dimensionless frequency =825.

We now investigate the case of isolated solutions in detail. Substitute into (7); the amplitude corresponding to the fixed points A, B, and C in Figure 6(c) can be calculated. The coordinates of fixed point A , B , and C are marked in Figure 6(c) with red letters.

The stability of the fixed points A, B, and C could be calculated by using the method in Section 2.4. Eigenvalues at each fixed point under the isolated solution case are listed as follows: A, B, and C. From the signs of eigenvalue it can be found that fixed point A is stable and the remaining two are unstable.

Likewise, for the case of single solution shown in Figure 6(d), it is found that the fixed point A is with eigenvalue , which is stable. By comparing the vector fields of Figures 6(c) and 6(d), the following phenomena are found:

Isolated solution includes special fixed points such as stable node A, saddle point B, and unstable node C, while the vector field diagram of single solution has only stable node A.

Under (0-2π) of phase θ, any initial positions will eventually reach the stable node A, but the transition process from different initial positions to stable node A is not the same. In the case of isolated solution, due to the existence of more special points, the transition process becomes more complex.

From vector field diagram, Figure 6(c), it can also be found that the state starting near the unstable node C will flow outward. The amplitude of the system may appear to decrease, increase, or remain the same temporarily and so on. However, the system amplitude will gradually increase until reaching to the stable node A.

When the system arrives near the saddle point B in the vector field diagram, Figure 6(c), it may move above or below fixed point B. And the system vibration amplitude will eventually reach to the stable node A or the stable node in the adjacent phase period.

3.3.2. Comparison of Transition Processes with Different Initial Conditions

In order to compare the changes and difference of VIV transition process response between the multiple solutions and the single solution under different initial conditions, we select some different initial points under dimensionless condition and set dimensionless frequency =825. We study two parameter cases with C_L0=0.01, C₁=8 and C_L0=0.06, C₁=2.1. In Figure 7, the blue curves show the displacement responses of multiple solutions condition and black curves show the displacement responses of single solution condition.

(a) =0.046, θ=0, multiple solution case

(b) =0.046, θ=0, single solution case

(c) =0.135, θ=0, multiple solution case

(d) =0.135, θ=0, single solution case

(e) =0.38, θ=0, multiple solution case

(f) =0.38, θ=0, single solution case

Figure 7 

Comparisons of displacement time history of transition process between isolated solution and single solution under different initial conditions.

From different initial condition setups, we find the following.

Multiple solutions scenario takes longer to steady state compared with the single solution case. In conjunction with the vector field diagram in Figure 6, it can be found that when the initial displacement is small and initial phase is close to 0 or 2π (Figure 7(a)), the system amplitude maintains slight stable vibration over a period of time. Even when the amplitude decreases at first, the system can still quickly increase to reach the stable node A after a while. Under the same condition for single solution case (Figure 7(b)), the amplitude will reach to minimum near zero and then gradually increase to steady state. We should be aware that the vibration period of transition process in the actual physical domain is much larger than in the dimensionless scale. For the case of Figure 7(a), the riser vibration may deceptively show a nearly stable pattern for a relatively long while at first. And then the riser system may experience a sudden growth of vibration; this could result in structure damage.

When the initial point is located near the saddle point B in the case of multiple solutions (Figure 7(c)), the amplitude of the riser system evolves quite smoothly. The time to reach the stable state is obviously longer than that of the single solution. Under the same initial conditions, the amplitude of the riser system will decrease firstly and then increase dramatically in the single solution (Figure 7(d)) with much shorter time.

When the initial amplitude is large and the initial phase is small, the amplitude will decrease firstly and then increase for both types of solution. In the case of single solution (Figure 7(f)), the amplitude of the riser system decreases rapidly to minimum near zero and then quickly increases to the stable node A. From Figure 7(e), it can be found that the minimum amplitude of the transition process is located near the saddle point B for multiple solutions case. The whole transition process time of multiple solutions case is significantly longer than single solution case. Attention should be paid to the dynamical response shown in Figure 7(f) due to the rapid change of response. Frequent oscillating of vibration amplitude may affect the structure safety of the riser for the long run.

To put all information together, it can be seen that due to the existence of the unstable node C in the case of multiple solutions, the system will maintain small amplitude vibration for a long time at first, and then the vibration increases significantly, which is worth enough vigilance and attention when designing the structure.

4. Conclusions

This paper constructs the nonlinear dynamic equation of VIV of TTR based on the third-order variable lift coefficient hydrodynamic model. The amplitude-frequency response equation under 1:1 primary resonance is developed based on the single mode discrete results. The nonlinear dynamic response characteristics of the TTR vibration under the excitation of vortex-induced resonance are discussed by the analytical solution. The analytical solution obtained in this paper matches the results of the numerical calculation and experimental data with decent accuracy.

Singularity analysis of the TTR system is carried out on the basis of the analytical solution. Influence of different lift coefficient parameters on riser amplitude is studied for key parameters C_L0 and C₁ in the lift coefficient model. Singularity study suggests the boundaries in the C_L0-C₁ plane, which separate the nonlinear behaviors of the riser system into different topological patterns. Both stable and unstable types of responses are discovered within different subregions divided by the transition sets. Numerical implementations agree well with analytical solutions obtained from the averaging method for stable solutions. Comparisons of the time domain responses under different initial conditions are then performed with the discussion of potential risks concealed under typical patterns. This paper studies the nonlinear vibration mechanism of riser system that undergoes vortex shedding excitation. The analysis results can provide guidance for the analytical treatment of VIV problem of top-tensioned risers.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the National Basic Research Program of China (973 Program) Grant No. 2014CB046805 and the National Natural Science Foundation of China (Grant No. 11672349).