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

Authors: Yuancen Wang; Zhiqiang Wu; Xiangyun Zhang
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 Figures1(a) and 1(b), respectively [8–10].Figure 1
Schematic diagram of vortex-induced vibration of the TTR system.
(a)
VIV of a riser subject to current (b)
Simplified force diagram of riserBased 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:(1)EI∂4yz,t∂z4-∂∂zTz∂yz,t∂z+mr+mf+ma∂2yz,t∂t2+C∂yz,t∂t=FyFy=12ρwDCLU2-12ρwDCDy˙y˙whereEI is the bending stiffness; y(z,t) is the transverse displacement of the riser.T is the top tension, which can be expressed as T(z)=Asg(ρs-ρw)(ftopL-z), in which As is the riser cross-sectional area, g is the acceleration of gravity, ρs and ρw are the riser material density and seawater density, respectively, ftop is the top tension coefficient, and L is the length of the top tension riser [11]; mr is the riser mass per unit length, mf is the riser fluid mass per unit length, ma is the additional mass per unit length (ma=CaρwπD2/4, in which Ca 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,CL is the lift coefficient, andCD is the drag coefficient. In this paper, the third-order variable lift coefficient model is used, which models the lift coefficient as CL(y˙(z,t))=CL0sin⁡(ωst)+C1y˙(z,t)+C2y˙2(z,t)+C3y˙3(z,t) [12], in which ωs is vortex shedding frequency (as ωs=2πStU/D, and St is the Strouhal number). The values of coefficients in the CL(y˙(z,t)) expression can be found from reference [7, 12], such asCL0=0.5,C1=1.82,C2=-1.29,C3=-0.71 orCL0= 0.22,C1= 1.62,C2=-2.31,C3= 0.75. Note that these coefficientsCL0,C1,C2, andC3 are dimensionless.Riser boundary conditions are as follows:(2)y0,t=yL,t=0,d2ydz2z=0,z=L=0


### 2.1. Dimensionless Equation

We introduce the following dimensionless quantities:y-=y/D, z-=z/L, and τ=ω·t (ω is time scale and has ω=EI/ML4). Let M=mr+mf+ma, by substituting the dimensionless transformation into (1), we get the dimensionless equation as follows:(3)∂4y-∂z-4+α1∂y-∂z--α1ftop-z-∂2y-∂z-2+∂2y-∂τ2+α2∂y-∂τ=α3CL0sin⁡ντω2D+α3C1ω∂y-∂τ+α3C2D∂y-∂τ2+α3C3ωD2∂y-∂τ3-α4∂y-∂τ∂y-∂τwhere α1=Asg(ρs-ρw)/Mω2L,α2=C/Mω,α3=ρwDU2/2M,α4=ρwD2CD/2M,ν=ωs/ω (dimensionless frequency).


### 2.2. Galerkin Discretization

Lety-=∑i=1Nφi(z-)qi(τ), qi(τ) is generalized coordinates, and let the ith order mode function φi(z-)=sin⁡(iπz-) 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, asN = 1. Multiply both side of (3) by the mode function φj(z-) and integrate both sides over the interval [0,1]. Using the orthogonality of the mode function, the ordinary differential equation with the single mode vibration of the TTR can be expressed as [13](4)q¨+ω12q+β2q˙-β3q˙2-β4q˙3+β5q˙q˙=β1sin⁡ντin which ω12=π4+(ftop-1/2)π2α1, β1=4CL0/πω2Dα3, β2=α2-C1/ωα3, β3=8C2D/3πα3, β4=3C3ωD2/4α3, β5=8/3πα4.


### 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 ν≈ω1 [17–19]. Assume the solution of (4) as follows:(5)q=a-·cos⁡ντ+θ=a-·cos⁡φq˙=-a-ν·sin⁡ντ+θ=-a-ν·sin⁡φwhere ν=ω1-εσ with ε being the small quantity that characterizes the perturbation and φ=ντ+θ.Note the fact thatq¨=-da-/dτνsin⁡ντ+θ-a-νcos⁡(ντ+θ)·(ν+dθ/dτ); by substituting this expression into (4) and (5), we have(6)da-dτ=-β2a-2+38β4a-3ν2-43πβ5a-2ν-β12νcos⁡θdθdτ=ω12-ν22ν+β12a-ω1·sin⁡θIn order to determine the amplitude-frequency response equations in steady state, let us setda-/dτ=0,dθ/dτ=0. 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 sin2θ+cos2θ=1.(7)ν2-ω122a-2β12+4ν2β12-β2a-2+38β4a-3ν2-43πβ5a-2ν2=1With given dimensionless frequencyν, (7) enables the calculation of steady state amplitude a- via solving nonlinear algebraic equation, when ν is close to ω1, which implies the occurrence of near resonance vibration. As a result, one can obtain the fixed point (a-,θ) 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 J=∂Y1/∂a-∂Y1/∂θ∂Z1/∂a-∂Z1/∂θ can be obtained from (6) withY1 andZ1 representing the right-hand side of the differential equations. For the first-order approximation solution of the averaging method, da-/dτ=Y1dθ/dτ=Z1, the characteristic equation can be written as follows:(8)∂Y1∂a-∂Y1∂θ∂Z1∂a-∂Z1∂θ-λΕ=∂Y1∂a--λ∂Y1∂θ∂Z1∂a-∂Z1∂θ-λ=λ2-∂Y1∂a-+∂Z1∂θλ+∂Y1∂a-·∂Z1∂θ-∂Z1∂a-·∂Y1∂θ=0Stability criterion of the system at each of its fixed points can be determined via its eigenvaluesλ1 and λ2 computed from (8). The eigenvalues of stable solution should satisfy (9)λ1+λ2=∂Y1∂a-+∂Z1∂θ<0λ1·λ2=∂Y1∂a-·∂Z1∂θ-∂Y1∂θ·∂Z1∂a->0


## 2.1. Dimensionless Equation

We introduce the following dimensionless quantities:y-=y/D, z-=z/L, and τ=ω·t (ω is time scale and has ω=EI/ML4). Let M=mr+mf+ma, by substituting the dimensionless transformation into (1), we get the dimensionless equation as follows:(3)∂4y-∂z-4+α1∂y-∂z--α1ftop-z-∂2y-∂z-2+∂2y-∂τ2+α2∂y-∂τ=α3CL0sin⁡ντω2D+α3C1ω∂y-∂τ+α3C2D∂y-∂τ2+α3C3ωD2∂y-∂τ3-α4∂y-∂τ∂y-∂τwhere α1=Asg(ρs-ρw)/Mω2L,α2=C/Mω,α3=ρwDU2/2M,α4=ρwD2CD/2M,ν=ωs/ω (dimensionless frequency).


## 2.2. Galerkin Discretization

Lety-=∑i=1Nφi(z-)qi(τ), qi(τ) is generalized coordinates, and let the ith order mode function φi(z-)=sin⁡(iπz-) 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, asN = 1. Multiply both side of (3) by the mode function φj(z-) and integrate both sides over the interval [0,1]. Using the orthogonality of the mode function, the ordinary differential equation with the single mode vibration of the TTR can be expressed as [13](4)q¨+ω12q+β2q˙-β3q˙2-β4q˙3+β5q˙q˙=β1sin⁡ντin which ω12=π4+(ftop-1/2)π2α1, β1=4CL0/πω2Dα3, β2=α2-C1/ωα3, β3=8C2D/3πα3, β4=3C3ωD2/4α3, β5=8/3πα4.


## 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 ν≈ω1 [17–19]. Assume the solution of (4) as follows:(5)q=a-·cos⁡ντ+θ=a-·cos⁡φq˙=-a-ν·sin⁡ντ+θ=-a-ν·sin⁡φwhere ν=ω1-εσ with ε being the small quantity that characterizes the perturbation and φ=ντ+θ.Note the fact thatq¨=-da-/dτνsin⁡ντ+θ-a-νcos⁡(ντ+θ)·(ν+dθ/dτ); by substituting this expression into (4) and (5), we have(6)da-dτ=-β2a-2+38β4a-3ν2-43πβ5a-2ν-β12νcos⁡θdθdτ=ω12-ν22ν+β12a-ω1·sin⁡θIn order to determine the amplitude-frequency response equations in steady state, let us setda-/dτ=0,dθ/dτ=0. 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 sin2θ+cos2θ=1.(7)ν2-ω122a-2β12+4ν2β12-β2a-2+38β4a-3ν2-43πβ5a-2ν2=1With given dimensionless frequencyν, (7) enables the calculation of steady state amplitude a- via solving nonlinear algebraic equation, when ν is close to ω1, which implies the occurrence of near resonance vibration. As a result, one can obtain the fixed point (a-,θ) 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 J=∂Y1/∂a-∂Y1/∂θ∂Z1/∂a-∂Z1/∂θ can be obtained from (6) withY1 andZ1 representing the right-hand side of the differential equations. For the first-order approximation solution of the averaging method, da-/dτ=Y1dθ/dτ=Z1, the characteristic equation can be written as follows:(8)∂Y1∂a-∂Y1∂θ∂Z1∂a-∂Z1∂θ-λΕ=∂Y1∂a--λ∂Y1∂θ∂Z1∂a-∂Z1∂θ-λ=λ2-∂Y1∂a-+∂Z1∂θλ+∂Y1∂a-·∂Z1∂θ-∂Z1∂a-·∂Y1∂θ=0Stability criterion of the system at each of its fixed points can be determined via its eigenvaluesλ1 and λ2 computed from (8). The eigenvalues of stable solution should satisfy (9)λ1+λ2=∂Y1∂a-+∂Z1∂θ<0λ1·λ2=∂Y1∂a-·∂Z1∂θ-∂Y1∂θ·∂Z1∂a->0


## 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 U=0.85m/s for comparison.The black dashed line in Figure2 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 atU=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.
parameter symbol values Elastic modulus E 210Gpa Seawater density ρ w 1025kg/m3 Material density ρ s 7850kg/m3 Internal fluid density ρ f 800kg/m3 Outside diameter D 0.325m Inside diameter D i 0.305m Riser length L 1500m Top tension coefficient f t o p 1.3 Additional mass coefficient C a 1.0 Strouhal number S t 0.2


### 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. Figure3 shows the bifurcation diagram of the TTR vortex-induced vibration response. Amplitude A varies with the external excitation frequency ωs 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ωs and the reduced velocity Ur.The reduced velocityUr is also an important reference for studying the VIV of the cylinder. In this section the explicit expression of the reduced velocity is Ur=U/f1·D, in which f1 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 Figure3 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 asCL=CL0sin⁡(ωst)+C1y˙(z,t)+C2y˙2(z,t)+C3y˙3(z,t); for the calculation of rigid riser and flexible riser, the choice of parameters inCL 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 parametersC2 andC3 have very little effects onCL. The influences brought byC2 andC3 are orders of magnitude smaller than that ofCL0 andC1. Therefore, in this paper, we only study the impact of varyingCL0 andC1 on the riser nonlinear dynamics via singularity theory.C2 andC3 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 B, the hysteresis setH, and the double limit setDL [24]. The bifurcation equation based on (7) can be expressed as (10)Ga-,ν,CL0,C1=ν2-ω122a-2β12+4ν2β12-β2a-2+38β4a-3ν2-43πβ5a-2ν2-1where a- and ν are treated as the state variable and the bifurcation parameter, respectively; CL0 and C1 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 (CL0, C1) plane that cateogize different solution tolopologies.(11)B:Ga-,ν,CL0,C1=Ga-a-,ν,CL0,C1=Gνa-,ν,CL0,C1v=0,a-≥0;H:Ga-,ν,CL0,C1=Ga-a-,ν,CL0,C1=Ga-a-a-,ν,CL0,C1a-=0,a-≥0;DL:Ga-i,ν,CL0,C1=Ga-a-i,ν,CL0,C1a-=0,a-1≠a-2,a-i≥0;by replacing the dimensionless frequencyν and the dimensionless amplitude a- with vortex shedding frequency ωs and amplitude a 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.(12)B:CL0≈7.56·10-7·251-2000ωs·251-1000ωs1.5ωs5C1≈0.157·251-2000ωs0.5ωs2.5(13)H:CL0≈6.21·10-5·2000ωs-2512ωs2C1≈0.0560·2000ωs-251ωs2.5From the above formula we can draw the boundaries of transition set as shown in Figure4. 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 divideCL0 -C1 plane into three regions (①, ②, ③).Figure 4
Transition set of parametersCL0 andC1 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 ofCL0 andC1 (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.Figure 5
Amplitude-frequency response of the riser VIV calculated with parametersCL0 andC1 in different regions.
(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 ③)From Figure5 we find that whenC1 stays unchanged, the amplitude will reduce withCL0 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 largeC1/CL0 ratio, which corresponds to region ③, the amplitude-frequency response may have multiple solutions. It can be found from the lift coefficient expression thatCL0 represents forced excitation due to vortex shedding andC1 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 Figure5, 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 withCL0=0.01,C1=8 andCL0=0.06,C1=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.Figure 6
Comparisons of amplitude-frequency response and vector field diagram between the isolated solution and single solution at dimensionless frequencyν=825.
(a)
Amplitude-frequency response of isolated solution (CL0=0.01,C1=8, and ν=825) (b)
Amplitude-frequency response of single solution (CL0=0.06, C1=2, and ν=825) (c)
Vector field diagram of isolated solution (CL0=0.01,C1=8, and ν=825) (d)
Vector field diagram of general single solution (CL0=0.06, C1=2, and ν=825)We now investigate the case of isolated solutions in detail. Substituteν=825 into (7); the amplitude a- corresponding to the fixed points A, B, and C in Figure 6(c) can be calculated. The coordinates of fixed point A (a-A,θA)=(0.26,3.735), B (a-B,θB)=(0.183,5.879), and C (a-C,θC)=(0.052,6.171) 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 Section2.4. Eigenvalues at each fixed point under the isolated solution case are listed as follows: A(λA1=-25.189,λA2=-2.729), B(λB1=-11.736,λB2=3.886), and C(λC1=11.797,λC2=14.683). 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 Figure6(d), it is found that the fixed point A (a-A,θA)=(0.265,3.236) is with eigenvalue (λA1=-41.38,λA2=-18.477), which is stable. By comparing the vector fields of Figures 6(c) and 6(d), the following phenomena are found:( 1 ) 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.( 2 ) 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.( 3 ) 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.( 4 )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 withCL0=0.01,C1=8 andCL0=0.06,C1=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.Figure 7
Comparisons of displacement time history of transition process between isolated solution and single solution under different initial conditions.
(a)
a -=0.046, θ=0, multiple solution case (b)
a -=0.046, θ=0, single solution case (c)
a -=0.135, θ=0, multiple solution case (d)
a -=0.135, θ=0, single solution case (e)
a -=0.38, θ=0, multiple solution case (f)
a -=0.38, θ=0, single solution caseFrom different initial condition setups, we find the following.( 1 ) 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.( 2 ) 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.( 3 ) 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.


## 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. Figure3 shows the bifurcation diagram of the TTR vortex-induced vibration response. Amplitude A varies with the external excitation frequency ωs 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ωs and the reduced velocity Ur.The reduced velocityUr is also an important reference for studying the VIV of the cylinder. In this section the explicit expression of the reduced velocity is Ur=U/f1·D, in which f1 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 Figure3 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 asCL=CL0sin⁡(ωst)+C1y˙(z,t)+C2y˙2(z,t)+C3y˙3(z,t); for the calculation of rigid riser and flexible riser, the choice of parameters inCL 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 parametersC2 andC3 have very little effects onCL. The influences brought byC2 andC3 are orders of magnitude smaller than that ofCL0 andC1. Therefore, in this paper, we only study the impact of varyingCL0 andC1 on the riser nonlinear dynamics via singularity theory.C2 andC3 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 B, the hysteresis setH, and the double limit setDL [24]. The bifurcation equation based on (7) can be expressed as (10)Ga-,ν,CL0,C1=ν2-ω122a-2β12+4ν2β12-β2a-2+38β4a-3ν2-43πβ5a-2ν2-1where a- and ν are treated as the state variable and the bifurcation parameter, respectively; CL0 and C1 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 (CL0, C1) plane that cateogize different solution tolopologies.(11)B:Ga-,ν,CL0,C1=Ga-a-,ν,CL0,C1=Gνa-,ν,CL0,C1v=0,a-≥0;H:Ga-,ν,CL0,C1=Ga-a-,ν,CL0,C1=Ga-a-a-,ν,CL0,C1a-=0,a-≥0;DL:Ga-i,ν,CL0,C1=Ga-a-i,ν,CL0,C1a-=0,a-1≠a-2,a-i≥0;by replacing the dimensionless frequencyν and the dimensionless amplitude a- with vortex shedding frequency ωs and amplitude a 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.(12)B:CL0≈7.56·10-7·251-2000ωs·251-1000ωs1.5ωs5C1≈0.157·251-2000ωs0.5ωs2.5(13)H:CL0≈6.21·10-5·2000ωs-2512ωs2C1≈0.0560·2000ωs-251ωs2.5From the above formula we can draw the boundaries of transition set as shown in Figure4. 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 divideCL0 -C1 plane into three regions (①, ②, ③).Figure 4
Transition set of parametersCL0 andC1 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 ofCL0 andC1 (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.Figure 5
Amplitude-frequency response of the riser VIV calculated with parametersCL0 andC1 in different regions.
(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 ③)From Figure5 we find that whenC1 stays unchanged, the amplitude will reduce withCL0 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 largeC1/CL0 ratio, which corresponds to region ③, the amplitude-frequency response may have multiple solutions. It can be found from the lift coefficient expression thatCL0 represents forced excitation due to vortex shedding andC1 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 Figure5, 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 withCL0=0.01,C1=8 andCL0=0.06,C1=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.Figure 6
Comparisons of amplitude-frequency response and vector field diagram between the isolated solution and single solution at dimensionless frequencyν=825.
(a)
Amplitude-frequency response of isolated solution (CL0=0.01,C1=8, and ν=825) (b)
Amplitude-frequency response of single solution (CL0=0.06, C1=2, and ν=825) (c)
Vector field diagram of isolated solution (CL0=0.01,C1=8, and ν=825) (d)
Vector field diagram of general single solution (CL0=0.06, C1=2, and ν=825)We now investigate the case of isolated solutions in detail. Substituteν=825 into (7); the amplitude a- corresponding to the fixed points A, B, and C in Figure 6(c) can be calculated. The coordinates of fixed point A (a-A,θA)=(0.26,3.735), B (a-B,θB)=(0.183,5.879), and C (a-C,θC)=(0.052,6.171) 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 Section2.4. Eigenvalues at each fixed point under the isolated solution case are listed as follows: A(λA1=-25.189,λA2=-2.729), B(λB1=-11.736,λB2=3.886), and C(λC1=11.797,λC2=14.683). 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 Figure6(d), it is found that the fixed point A (a-A,θA)=(0.265,3.236) is with eigenvalue (λA1=-41.38,λA2=-18.477), which is stable. By comparing the vector fields of Figures 6(c) and 6(d), the following phenomena are found:( 1 ) 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.( 2 ) 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.( 3 ) 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.( 4 )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 withCL0=0.01,C1=8 andCL0=0.06,C1=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.Figure 7
Comparisons of displacement time history of transition process between isolated solution and single solution under different initial conditions.
(a)
a -=0.046, θ=0, multiple solution case (b)
a -=0.046, θ=0, single solution case (c)
a -=0.135, θ=0, multiple solution case (d)
a -=0.135, θ=0, single solution case (e)
a -=0.38, θ=0, multiple solution case (f)
a -=0.38, θ=0, single solution caseFrom different initial condition setups, we find the following.( 1 ) 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.( 2 ) 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.( 3 ) 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.


## 3.3.1. Evolution of Vector Field

Under dimensionless condition two parameters set withCL0=0.01,C1=8 andCL0=0.06,C1=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.Figure 6
Comparisons of amplitude-frequency response and vector field diagram between the isolated solution and single solution at dimensionless frequencyν=825.
(a)
Amplitude-frequency response of isolated solution (CL0=0.01,C1=8, and ν=825) (b)
Amplitude-frequency response of single solution (CL0=0.06, C1=2, and ν=825) (c)
Vector field diagram of isolated solution (CL0=0.01,C1=8, and ν=825) (d)
Vector field diagram of general single solution (CL0=0.06, C1=2, and ν=825)We now investigate the case of isolated solutions in detail. Substituteν=825 into (7); the amplitude a- corresponding to the fixed points A, B, and C in Figure 6(c) can be calculated. The coordinates of fixed point A (a-A,θA)=(0.26,3.735), B (a-B,θB)=(0.183,5.879), and C (a-C,θC)=(0.052,6.171) 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 Section2.4. Eigenvalues at each fixed point under the isolated solution case are listed as follows: A(λA1=-25.189,λA2=-2.729), B(λB1=-11.736,λB2=3.886), and C(λC1=11.797,λC2=14.683). 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 Figure6(d), it is found that the fixed point A (a-A,θA)=(0.265,3.236) is with eigenvalue (λA1=-41.38,λA2=-18.477), which is stable. By comparing the vector fields of Figures 6(c) and 6(d), the following phenomena are found:( 1 ) 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.( 2 ) 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.( 3 ) 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.( 4 )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 withCL0=0.01,C1=8 andCL0=0.06,C1=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.Figure 7
Comparisons of displacement time history of transition process between isolated solution and single solution under different initial conditions.
(a)
a -=0.046, θ=0, multiple solution case (b)
a -=0.046, θ=0, single solution case (c)
a -=0.135, θ=0, multiple solution case (d)
a -=0.135, θ=0, single solution case (e)
a -=0.38, θ=0, multiple solution case (f)
a -=0.38, θ=0, single solution caseFrom different initial condition setups, we find the following.( 1 ) 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.( 2 ) 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.( 3 ) 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 parametersCL0 andC1 in the lift coefficient model. Singularity study suggests the boundaries in the CL0-C1 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.



---

*Source: 6491517-2018-12-18.xml*