The Instability and Response Studies of a Top-Tensioned Riser under Parametric Excitations Using the Differential Quadrature Method

: The differential quadrature method (DQM) is a numerical technique widely applied in structure mechanics problems. In this work, a top-tensioned riser conveying ﬂuid is considered. The governing equation of this riser under parametric excitations is deduced. Through Galerkin’s method, the partial differential governing equation with respect to time t and vertical coordinate z is reduced into a 1D differential equation with respect only to time. Moreover, the DQM is applied to discretize the governing equation to give solution schemes for the risers’ parametric vibration problem. Furthermore, the instability region of Mathieu equation is studied by both the DQM and the Floquet theory to verify the effectiveness of the DQM, and the solutions of both methods show good consistency. After that, the inﬂuences of some factors such as damping coefﬁcient, internal ﬂow velocity, and wet-weight coefﬁcient on the parametric instability of a top-tensioned riser are discussed through investigating the instability regions solved by the DQM solution scheme. Hence, conclusions are obtained that the increase of damping coefﬁcient will save the riser from parametric resonance while increasing internal ﬂow velocity, or the wet-weight coefﬁcient will deteriorate the parametric instability of the riser. Finally, the time-domain responses of several speciﬁc cases in both stable region and unstable region are presented.


Introduction
The risers are pipelines linking the seabed and the platform above the sea. Because of the vertically dynamic excitations induced by the heave motion of the platform, the parametrical resonance of the risers would occur, and fatigue damage will follow. To avoid such risk, it poses an urgent demand to make a parametric instability prediction plan at the designing phase of the risers.
The study of the dynamic properties of pipes dates back to the 1970s. From 1970 to 1972, Chen [1][2][3][4] presented a comprehensive investigation of the free vibration, forced vibration, and stability of pipes conveying fluids. Since then, the vibrations of pipes have been drawing attention from researchers. In 1975, Hsu [5] figured out, for the first time in the engineering application, the parametric instability characters of the straight and vertical configuration of a hanging. He analyzed the parametric instability of hanging string risers with a Ω-Γ instability chart and studied the influences of nonlinear fluid dynamic damping 2 of 23 on the steady-state responses of risers under parametrically unstable conditions. Later, a lot of classic contributions [6][7][8][9][10][11][12] were carried out about the parametric vibration of risers. Recently, Lou et al. [13] studied the horizontal parametric vibration of a compliant vertical access riser, and mode coupling and damping's effects on risers' parametric instability were discussed. Zhang et al. [14] investigates the vibration of a top tensioned riser under both parametrical and vortex-induced excitations by numerical simulation and experiment, and the aggravation of the riser's VIV response affected by parametric excitation has been evaluated. Liang and Lou [15] studied the effects of multiphase internal flows that consider hydrate phase transitions on the parametric stability of marine risers.
In the current study, a top-tensioned riser subjected to parametrical excitation is studied. For the dynamic problems in engineering applications, some numerical methods, such as the finite element method [8,[16][17][18][19][20], Runge-Kutta (R-K) method [7,10], and finite differences method [9], have been verified to be valid and accurate in coping with such problems. Ni et al. [21,22] considered the dynamic problems of a pipe conveying fluid by the differential quadrature method (DQM) and generalized differential quadrature rule (GDQR). They found that DQM and GDQR show good competency in such problems. Illuminated by their works, we decide to consider the parametric instability of the toptensioned riser with DQM. The DQM was first introduced by Bellman and Casti [23] as a new technique to solve the initial-value problems of ordinary and partial differential equations. Based on such work, DQM was extended into various applications such as fluid flow and turbulence, the Hodgkin-Huxley model, structural mechanics and mechanical dynamics, etc. [24][25][26][27][28]. Furthermore, Wu and Liu [29,30] complement the DQM by introducing the generalized differential quadrature rules (GDQR) to equip this method with the ability to treat high-order differential equations and complex boundary conditions more flexibly. In recent years, the study on DQM continues. Essam et al. [31] proposed a new hybrid linearization-differential quadrature method (HLDQM), and some numerical experiments on a Fe3O4 ferrofluid flow over a vertical radiate plated were conducted by this method. Constantin et al. [32] applied the least-squares differential quadrature method (LSDQM) to nonlinear partial differential equations, and solutions with better accuracy were acquired compared to other well-known methods: homotopy perturbation method (HPM) and the Adomian decomposition method (ADM). Additionally, some novel modifications on DQM have been made to accommodate more flexible situations [33][34][35][36].
Based on the aforementioned studies of the risers' parametric vibration and DQM, in the current study, the DQM would be introduced to calculate the vibration properties of risers under parametric excitation. Additionally, the effects of some key factors such as damping, internal flow velocity, and unit wet weight on the risers' parametric instability would be discussed in detail. In Section 2, the governing equation of the parametric vibration of a top-tensioned riser with internal flows and the implementation of solving this governing by DQM was presented, and the deduction of the unstable zone of the Mathieu equation is conducted. In Section 3, the validation of DQM in analyzing risers' parametric instability is verified with the Mathieu equation; also, the instability charts of a top-tensioned riser under different conditions are given to discuss the influences of different parameters on risers' instability. In Section 4, the parametrically excited responses of this riser are given. In Section 5, the conclusions are made.

Governing Equation of Parametric Vibration
In the ocean oil drilling industry, the riser is an indispensable structure transporting oil from the seabed to the platform. Our study is based on a top-tensioned one whose figure is illustrated in Figure 1. As shown in Figure 1, the riser is hinged at the bottom end and the top of the riser is pulled by tension cables. An array of buoyancy tanks is arranged along the riser to counteract the weight of the riser and provide extra tension to keep the riser erect. Since the use of risers is to transport fluid, there is liquid flowing inside risers. When the platform heaves up and down, the riser will suffer axial dynamic tension, which will generate parametric excitation to the riser.

Governing Equation of Parametric Vibration
In the ocean oil drilling industry, the riser is an indispensable structure transporting oil from the seabed to the platform. Our study is based on a top-tensioned one whose figure is illustrated in Figure 1. As shown in Figure 1, the riser is hinged at the bottom end and the top of the riser is pulled by tension cables. An array of buoyancy tanks is arranged along the riser to counteract the weight of the riser and provide extra tension to keep the riser erect. Since the use of risers is to transport fluid, there is liquid flowing inside risers. When the platform heaves up and down, the riser will suffer axial dynamic tension, which will generate parametric excitation to the riser. To transform this engineering model into mathematical expression, some hypotheses are introduced, and a simplification of the riser is given. The material, mechanical properties, and the buoyancy distribution are regarded as uniform and continuous along the length of the riser. Thus, the riser considered in our study can be simplified as a vertical Euler beam hinged at both ends with internal fluid flowing inside it and with uniformly-distributed tension along with it. The simplified structure of this and the force diagram of an infinitesimal unit of this riser are depicted in Figure 2.  To transform this engineering model into mathematical expression, some hypotheses are introduced, and a simplification of the riser is given. The material, mechanical properties, and the buoyancy distribution are regarded as uniform and continuous along the length of the riser. Thus, the riser considered in our study can be simplified as a vertical Euler beam hinged at both ends with internal fluid flowing inside it and with uniformlydistributed tension along with it. The simplified structure of this and the force diagram of an infinitesimal unit of this riser are depicted in Figure 2.
In the ocean oil drilling industry, the riser is an indispensable structure transporting oil from the seabed to the platform. Our study is based on a top-tensioned one whose figure is illustrated in Figure 1. As shown in Figure 1, the riser is hinged at the bottom end and the top of the riser is pulled by tension cables. An array of buoyancy tanks is arranged along the riser to counteract the weight of the riser and provide extra tension to keep the riser erect. Since the use of risers is to transport fluid, there is liquid flowing inside risers. When the platform heaves up and down, the riser will suffer axial dynamic tension, which will generate parametric excitation to the riser. To transform this engineering model into mathematical expression, some hypotheses are introduced, and a simplification of the riser is given. The material, mechanical properties, and the buoyancy distribution are regarded as uniform and continuous along the length of the riser. Thus, the riser considered in our study can be simplified as a vertical Euler beam hinged at both ends with internal fluid flowing inside it and with uniformly-distributed tension along with it. The simplified structure of this and the force diagram of an infinitesimal unit of this riser are depicted in Figure 2. w L (t) The riser is simplified into a 2D beam subjecting to the horizontal forcing displacement excitation u L (t) and the vertical parametric excitation T L (t) at the top end. For this infinitesimal unit, there is tension T and shear force N at both ends. m w g is the unit wet weight, and it equals the difference between the total weight of the unit and the buoyancy of the unit, namely, where f b is the buoyancy provided by the riser and the buoyancy tanks of unit length, and m s and m f are, respectively, pipe weight and internal fluid weight of unit length: where D is the outer diameter of the riser; h is the thickness of the pipe wall; ρ s is the riser wall density, and ρ f is the internal fluid density. To simplify the model, the wet-weight m w g is defined with an adjustable coefficient k mw because the wet weight can be adjusted by the size of buoyancy tanks: U is the flowing velocity of the internal fluid. f nd and f τd are damping forces at normal and tangential directions of the unit. According to the force diagram in Figure 2b, the force equilibrium of this infinitesimal unit could be expressed as [10] in the axial direction: and in the transverse direction: where u and w, respectively, are displacements of the unit in the axial direction and transverse direction; m a represents added mass of unit length, and here, the added mass is hypothesized as the weight of the water in a volume of the unit length of riser [11]: where ρ w is the density of water; D 2 /Dt 2 is the material derivation of the internal fluid: and here, linear damping is utilized [11]: In our study, the riser is assumed as a Euler beam, which means the bending angle ϕ is close to zero, and axial displacement u is quite smaller compared to transverse displacement w. Thus, Equation (4) can be simplified into a static equilibrium: Mathematics 2022, 10, 1331

of 23
Considering small elastic deformation and transverse deflection, the coordinate s could be replaced with vertical coordinate z. As the tension at the top is T L (t), the tension distribution along the riser could be represented as where T 0 is the static pretension at the top, and T d is the dynamic tension at the top. In a Euler beam, the shear force N could be approximated with where EI is the stiffness of the riser. Thus, Equation (6) can be simplified as where M = m s + m f + m a . Equation (14) is a partial differential equation with respect to time t and vertical position z. In the following section, solution schemes for this equation will be discussed.

Galerkin's Method and Derivation of Mathieu Equation
As the governing equation is a partial differential equation with two variables, it is advisable to reduce this govern equation into a differential equation with only one variable. Thus, Galerkin's method and modes' superposition principle could be applied. It can be assumed that the deflection function w(z,t) could be represented with a weighted summation of the orthogonal basis: where ϕ i is the basic shape function of ith mode, and q i is the weight function of ith mode; ϕ 0 is the rigid shape function, which means there is no deformation in this mode, and the displacement of the riser is linear along the vertical coordinate, namely, Substituting Equation (15) into Equation (14), where Here in our study, the mode functions are a series of sine functions: According to the orthogonality of trigonometric functions, π 0 sin(mt) sin(nt) = 0, m = n.
Mathematics 2022, 10, 1331 6 of 23 Thus, applying this orthogonality to the mode functions, With this orthogonality of trigonometric functions, Equation (17) could be reduced into a one-dimensional differential equation by integrating the partial equation along the coordinate z from 0 to L with the weight function ϕ j : In the current research, we assume the parametric excitation T d to be a cosine function with respect to time t: where S is the amplitude of parametric excitation, and Ω is the excitation frequency. Replacing the time t with dimensionless term τ, Substituting Equations (24) and (25) into Equation (23) and evaluating the integrations, the governing equation could be reduced to In Equation (26), ζ is the damping term, µ 1 and µ 2 are wet-weight terms, and µ 3 and µ 4 are internal flow terms. Moreover, P ij and Q ij are coupling coefficient matrices. Especially, when the damping, wet weight, internal flow terms, and forcing excitation are omitted, Equation (26) will be reduced to the classic Mathieu equation:

The Differential Quadrature Solution Scheme for Risers' Parametric Vibration
According to the studies of previous researchers, the DQM is proved to be a powerful technique to solve the initial-value problems. Here, we will give a brief introduction to the DQM.
f (x) is a function in the domain of [a, b] with nth order continuous derivatives with respect to x. According to the Lagrange interpolation formula, f (x) could be approximated to a weighted summation: where Thus, the nth order derivative of f (x) can also be approximated to the derivative of this interpolation function: If {x k } are given nodes, the derivatives of f (x) on these given nodes can be expressed as where With this interpolation and derivation scheme, the value of f (x) and its derivatives can be expressed with weighted summations of the function values at several discrete nodes. Sherbourne and Pandey [25] found that the DQ solution is very sensitive to the node spacing. Runge's phenomenon may easily occur when the nodes are uniformly distributed. Thus, the Chebyshev node spacing scheme is utilized: Then, q i and its derivatives on these node points could be expressed as Mathematics 2022, 10, 1331 8 of 23 where Substituting Equations (41)-(43) to Equation (26), where In Equation (45), index i is the order number, and index k represents the node number of dimensionless time coordinate τ. The discretized differential equation is obtained by using Equation (45), but the equation is underdetermined if the boundary conditions are not given [37]. To avoid the right end of Equation (45) to be constant at zero, the initial value for each mode is given as 0.001, namely The initial velocity of each mode is set to zero: In discretized form, The initial condition can be embedded to the discretized differential equation set Equation (45) by replacing the equations of index l = 1, N with Equations (49) and (50). However, such a replacement would raise the problem that the properties at these replaced indexes will be overlooked. According to [26,27,38], adjacent δ points can be arranged, and the mechanical properties of these replaced points can be approximated by those of adjacent δ points. Additionally, the closeness between the adjacent point and the replaced point is usually set to be of the order δ = 10 −5 . Thus, the Chebyshev node spacing scheme Equation (40) should be reorganized with adjacent δ points: As mentioned in many previous works, for a large number of node points and a large-ordered algebraic equation system, the matrix of DQM will be ill-conditioned [37,39]. However, for the current problem, a large number of node points is needed to acquire a long-time scope numerical solution. Thus, the block-marching technique [39] is introduced to overcome such conflict between the long-time solution and the ill-conditioned matrix. As the level-by-level solution scheme in R-K schemes, the block-marching of DQM divides the temporal problem domain into a sequence of blocks, and the solving procedure is organized block-by-block, which means the initial conditions of each lateral block are determined by the solution of the former block, as shown in Figure 3. For instance, the solution of q i in the  2T], the initial conditions can be defined with the q and q values on the last node of the first block, and the equation set of the second block can be expressed as Mathematics 2022, 10, x FOR PEER REVIEW 10 of 25 Similarly, the sequential blocks can be solved by this block-marching scheme.

Analytical Solution for Mathieu Instability
According to the Floquet theory [40], we have the the Mathieu equation: having a pair of solution basis where p1 and p2 are periodic functions with period π, and σ1 and σ2 are complex coefficients. Additionally, σ1 and σ2 are the roots of the quadratic Equation [40]: namely, When |a| < 1, both σ1 and σ2 are complex roots with modulus 1, and the solution of Equation (56) is stable with respect to τ. On the contrary, when |a| > 1, moduli of σ1,2 are bigger than 1, and that means the solution of Equation (56) is unstable. For the situation |a| = 1, both σ1 and σ2 are equal to 1 or −1, and the solution of Equation (56) is periodic with period π or 2π. Therefore, |a| = 1 is the critical condition between the stable regions and unstable regions, and it means that for the cases at the critical edge between stable and unstable regions, the solutions for Equation (56) are periodic functions with period π or 2π.
To obtain the borderlines of the unstable region, firstly, we assume the solutions to be with period π. Thus, the assumed solution could be expanded with the Fourier series: Substituting Equation (59) into Equation (55), 4 cos 2 4 sin 2 n a n n b n The first block The second block The third block ……. Similarly, the sequential blocks can be solved by this block-marching scheme.

Analytical Solution for Mathieu Instability
According to the Floquet theory [40], we have the the Mathieu equation: having a pair of solution basis where p 1 and p 2 are periodic functions with period π, and σ 1 and σ 2 are complex coefficients. Additionally, σ 1 and σ 2 are the roots of the quadratic Equation [40]: namely, When |a| < 1, both σ 1 and σ 2 are complex roots with modulus 1, and the solution of Equation (56) is stable with respect to τ. On the contrary, when |a| > 1, moduli of σ 1,2 are bigger than 1, and that means the solution of Equation (56) is unstable. For the situation |a| = 1, both σ 1 and σ 2 are equal to 1 or −1, and the solution of Equation (56) is periodic with period π or 2π. Therefore, |a| = 1 is the critical condition between the stable regions and unstable regions, and it means that for the cases at the critical edge between stable and unstable regions, the solutions for Equation (56) are periodic functions with period π or 2π.
To obtain the borderlines of the unstable region, firstly, we assume the solutions to be with period π. Thus, the assumed solution could be expanded with the Fourier series: Substituting Equation (59) into Equation (55), Using the prosthaphaeresis rules, According to the orthogonality the trigonometric functions, and Substituting Equation (61) into Equation (60), and integrating this equation multiplied with weight function cos(2mτ) from 0 to π, a group of linear equation set is formed: Similarly, integrating the equation multiplied with weight function sin(2mτ): To ensure the solution of Equation (59) is a constant 0 value, the determinant of the coefficient matrix of Equations (64) and (65) must be zero. Namely, or For cases where the solution of Equation (56) is of period 2π, there exist the eigen equations

Verification and Validation of DQM
In this section, two simple equations, a first-order linear differential equation with harmonic coefficients and a second-order linear differential equation with constant coefficients, are considered: and The analytical solutions of these equations can be expressed as for Equation (70), and y = sin(x) for Equation (71). As discussed in Section 2.3, the initial-value problems can be studied by the blockmarching DQM. In the current solving scheme for these two equations, the computation domain x ∈ [0, 10 × 2π] is divided into 10 blocks with a length of 2π. Therefore, the long-time-scope initial problem of Equations (70) and (71) is divided into two sequences of short-time-scope subproblems that can be solved by DQM with a small number of nodes. Additionally, the DQM, whose node number N = 51, and analytical solutions for Equations (70) and (71) are given in Figure 4.  It can be acknowledged from Figure 4 that the solutions of DQM reach good agreements with those of the analytical solution. Furthermore, the convergence of DQM is evaluated with the average error where y i is the analytical value and y i is the numerical value; M is the block number and, in the current examples, M is set to 10; N is the DQM node number of each block. Under such conditions, the average errors with different node numbers are computed and presented in Figure 5.   As shown in Figure 5, the average errors of both examples diminish sharply with the increase of node number per block when N remains in relatively small values. Addi tionally, when the node number exceeds a certain value, which is 34 for Equation (70) or 18 for Equation (71), the decreasing trend of the average errors falls to a relatively sligh rate. However, when the node number increases further, the average error begins to al ter in an increasing trend until at N = 737, the DQM solution scheme aborts because o the ill-conditioned matrix. Thus, in application practice, a relatively small node numbe is chosen in DQM. In the following study, 51 nodes are utilized in each block to form the DQM solution scheme.
According to the previous studies, the R-K method is a commonly utilized method in coping with the initial value problems. Here, the R-K solutions for Equations (70) and (71) are carried out to make a comparison with the DQM solutions. The R-K solution schemes in our study for Equations (70) and (71) are: As shown in Figure 5, the average errors of both examples diminish sharply with the increase of node number per block when N remains in relatively small values. Additionally, when the node number exceeds a certain value, which is 34 for Equation (70) or 18 for Equation (71), the decreasing trend of the average errors falls to a relatively slight rate. However, when the node number increases further, the average error begins to alter in an increasing trend until at N = 737, the DQM solution scheme aborts because of the ill-conditioned matrix. Thus, in application practice, a relatively small node number is chosen in DQM. In the following study, 51 nodes are utilized in each block to form the DQM solution scheme.
According to the previous studies, the R-K method is a commonly utilized method in coping with the initial value problems. Here, the R-K solutions for Equations (70) and (71) are carried out to make a comparison with the DQM solutions. The R-K solution schemes in our study for Equations (70) and (71) are: for Equation (71).
The R-K solution scheme utilized in the current examples is an explicit single-step method. To ensure the same node number in each block as that in DQM, the step lengths for the current R-K schemes are: The solution of Equations (70) and (71) are solved by both DQM and R-K, and the average errors per block are evaluated and presented in Figure 6.
for Equation (71). The R-K solution scheme utilized in the current examples is an explicit single-step method. To ensure the same node number in each block as that in DQM, the step lengths for the current R-K schemes are: The solution of Equations (70) and (71) are solved by both DQM and R-K, and the average errors per block are evaluated and presented in Figure 6.  According to Figure 6, it can be informed that for these two examples the average error of the R-K method is much bigger than the DQM. That means with the same number of nodes the DQM could acquire more accurate solutions than the R-K. In addition, it can be told from the figures that the average error of the DQM for Equation (71) is smaller than Equation (70), while the trend of the R-K is the opposite. That means the causes of errors for these two methods are irrelevant.

DQM Solution for Mathieu Equation
In this section, the Mathieu equation (32)  For specific values of α and β as listed in Table 1, the accordingly time-domain solutions for the Mathieu equation are solved by both DQM and R-K, and the solutions are presented in Figure 7.

DQM Solution for Mathieu Equation
In this section, the Mathieu equation (32) is considered. Additionally, the initial values for the current examples are: For specific values of α and β as listed in Table 1, the accordingly time-domain solutions for the Mathieu equation are solved by both DQM and R-K, and the solutions are presented in Figure 7.   Firstly, according to Figure 7, the curves of DQM show good consistency with that of R-K, which is a vastly applied technique for the parametric vibrating problem. Additionally, it can be easily told from the figures that for Examples 1, 2, and 4, the solutions for the Mathieu equation are stable; however, for Example 2, the solution is divergent with respect to time τ. Following such logic, a series of Mathieu equation examples with α ranging from −5 to 20 and β ranging from 0 to 30 are calculated by this DQM solution scheme to form the instability chart of the Mathieu equation. When it comes to how to tell whether an example is divergent, we set up a verification scheme as follows:

1.
Separate the long-time-scope solutions into blocks with lengths of 2π; 2.
Extract the maximum of each block; 3.
Compare the maxima between every neighboring block, and if the maximum is increasing block by block, we define the solution for the current example as divergent;

4.
For the examples that fail the test in step three, calculate the ratio of the maximum of the last block and that of the initial block, and if the ratio is bigger than 100, we define the solution for the current example as divergent; otherwise, it is defined as convergent.
For the solutions defined as divergent, we mark the corresponding (α, β) point in the α-β plane. Thus, after all of the verification schemes of every example are considered, a map of the instability chart for the Mathieu equation is acquired. Additionally, this instability chart is acquired via the aforementioned scheme as well as the borderline solved by Equations (62)-(65) and is depicted in Figure 8. tell whether an example is divergent, we set up a verification scheme as follows: 1. Separate the long-time-scope solutions into blocks with lengths of 2π; 2. Extract the maximum of each block; 3. Compare the maxima between every neighboring block, and if the maximum is in creasing block by block, we define the solution for the current example as diver gent; 4. For the examples that fail the test in step three, calculate the ratio of the maximum of the last block and that of the initial block, and if the ratio is bigger than 100, we define the solution for the current example as divergent; otherwise, it is defined a convergent.
For the solutions defined as divergent, we mark the corresponding (α, β) point in the α-β plane. Thus, after all of the verification schemes of every example are considered a map of the instability chart for the Mathieu equation is acquired. Additionally, this in stability chart is acquired via the aforementioned scheme as well as the borderlin solved by Equations (62)-(65) and is depicted in Figure 8.
. As is shown in Figure 8, the instability region acquired by the DQM fits well with the borderline defined by the Floquet theory. Therefore, it can be easily concluded tha the DQM is a suitable method to evaluate the parametric vibration problems of risers. In the following sections, the parametric instability of risers under different conditions wil be discussed with the solutions of the DQM.

Instabilities of Risers with Different Damping Coefficients
In this section, the parametric instability of a top-tensioned riser with differen damping coefficients is discussed through the solutions of the DQM. Additionally, th riser we study is modified from a mini-TLP in the Gulf of Mexico [11] whose major pa rameters are listed in Table 2. The internal flow velocity U and wet-weight coefficient kmw are set to zero, which means that the fluid inside the pipe is assumed to be still, and the weight of the riser unit and the buoyancy provided by the tanks are assumed to b As is shown in Figure 8, the instability region acquired by the DQM fits well with the borderline defined by the Floquet theory. Therefore, it can be easily concluded that the DQM is a suitable method to evaluate the parametric vibration problems of risers. In the following sections, the parametric instability of risers under different conditions will be discussed with the solutions of the DQM.

Instabilities of Risers with Different Damping Coefficients
In this section, the parametric instability of a top-tensioned riser with different damping coefficients is discussed through the solutions of the DQM. Additionally, the riser we study is modified from a mini-TLP in the Gulf of Mexico [11] whose major parameters are listed in Table 2. The internal flow velocity U and wet-weight coefficients k mw are set to zero, which means that the fluid inside the pipe is assumed to be still, and the weight of the riser unit and the buoyancy provided by the tanks are assumed to be compensating for each other. The parametric excitation is defined through T d (t) = Scos(2πt/P), where the amplitude S ranges from 0 to 5 × 10 6 N and the period P ranges from 0 to 18s. Therefore, the instability charts of six different damping coefficients ζ = (0.00, 0.05, 0.10, 0.15, 0.20, 0.25) are acquired and are presented in Figure 9.
Because the coupling terms µ 1 P ij q i , µ 2 Q ij q i , and µ 3 Q ij q i ' are set to zero, there are no coupling effects between different modes. As a result, the parametric instabilities of different modes are independent of each other. In Figure 9, the instability regions of different damping coefficients of the first ten orders are presented. The instability regions are represented with translucent shades confined to the solid borderlines. It can be told from the pictures that with the increase of damping coefficients, the shaded areas move upwards accordingly. It means that the risers with larger damping coefficients are less prone to vibrate unstably Besides, at a specific heave period, stronger dynamic tension amplitude is needed to induce the parametric instability of the riser. compensating for each other. The parametric excitation is defined through Td(t) = Scos(2πt/P), where the amplitude S ranges from 0 to 5 × 10 6 N and the period P ranges from 0 to 18s. Therefore, the instability charts of six different damping coefficients ζ = (0.00, 0.05, 0.10, 0.15, 0.20, 0.25) are acquired and are presented in Figure 9.  Because the coupling terms μ1Pijqi, μ2Qijqi, and μ3Qijqi' are set to zero, there are no coupling effects between different modes. As a result, the parametric instabilities of different modes are independent of each other. In Figure 9, the instability regions of different damping coefficients of the first ten orders are presented. The instability regions are

Instabilities of Risers with Different Internal Flow Velocity
According to Chen's [1][2][3][4] research, the velocity of the internal flow will alter the vibrating characters of the pipe. Due to the role of the risers being to transport fluid from seabed to the platform, in this section, the instabilities of the riser with different internal flow velocities are studied. The parameters of the riser and the environment are listed in Table 2. The damping coefficient ζ is set to 0.1, and the wet-weight coefficient k mw is set to zero. The instability charts of different internal flow velocities U, which range in array [0-5] m/s, are obtained and depicted in Figure 10.
vibrating characters of the pipe. Due to the role of the risers being to transport fluid from seabed to the platform, in this section, the instabilities of the riser with different internal flow velocities are studied. The parameters of the riser and the environment are listed in Table 2. The damping coefficient ζ is set to 0.1, and the wet-weight coefficient kmw is set to zero. The instability charts of different internal flow velocities U, which range in array [0-5] m/s, are obtained and depicted in Figure 10. The existence of the nonzero coupling term μ3Qijqi' would induce interferences between different modes. Therefore, when one order of the modes vibrates unstably, the divergent response of the current mode will gradually trigger other modes to vibrate unstably. Namely, when the internal flow velocity is not equal to zero, the μ3Qijqi' term would take effect crossing different modes. For instance, in case the first mode of the riser vibrates unstably in the shape of       π sin z L , the divergent vibration of the first order would pump continuous excitations to other modes because of the existence of the coupling term μ3Qijqi'. Thus, the parametrically unstable vibrations of other modes will be triggered. Therefore, the parametrical instabilities of different modes are coexistent, and such an instability chart will no longer be presented individually by mode. In Figure 10, it could be found that with the increase of internal flow velocity, the unstable regions move towards the bottom right. It can be informed that when the internal flow velocity increases, the risers are more easily to be parametrically unstable, and the minimal dynamic tension that could induce the parametric instability would decrease accordingly. Besides, it can also be spotted that the movement of the unstable borderline is accelerating with the increase of internal flow velocity. The acceleration phenomenon is due to the coefficient μ4, which is proportional to U 2 . Thus, the influence of the internal flow velocity on the riser's parametric instability is small at a low-frequency range where U is less than 2 m/s, while its influence would reach a considerable extent otherwise. The existence of the nonzero coupling term µ 3 Q ij q i ' would induce interferences between different modes. Therefore, when one order of the modes vibrates unstably, the divergent response of the current mode will gradually trigger other modes to vibrate unstably. Namely, when the internal flow velocity is not equal to zero, the µ 3 Q ij q i ' term would take effect crossing different modes. For instance, in case the first mode of the riser vibrates unstably in the shape of sin πz L , the divergent vibration of the first order would pump continuous excitations to other modes because of the existence of the coupling term µ 3 Q ij q i '. Thus, the parametrically unstable vibrations of other modes will be triggered. Therefore, the parametrical instabilities of different modes are coexistent, and such an instability chart will no longer be presented individually by mode. In Figure 10, it could be found that with the increase of internal flow velocity, the unstable regions move towards the bottom right. It can be informed that when the internal flow velocity increases, the risers are more easily to be parametrically unstable, and the minimal dynamic tension that could induce the parametric instability would decrease accordingly. Besides, it can also be spotted that the movement of the unstable borderline is accelerating with the increase of internal flow velocity. The acceleration phenomenon is due to the coefficient µ 4 , which is proportional to U 2 . Thus, the influence of the internal flow velocity on the riser's parametric instability is small at a low-frequency range where U is less than 2 m/s, while its influence would reach a considerable extent otherwise.

Instability of Risers with Different Wet-Weight Coefficients
According to the governing equation Equation (29), the terms α j q j , µ 1 P ij q i , and µ 2 Q ij q i , which are relative to the unit wet weight, are considered. Thus, the influences of wet weight on risers' parametric instability are unavoidable. Since the unit wet weight can be adjusted by changing the size of buoyancy tanks distributed along with the risers, a dimensionless coefficient k mw is used in our research to control the unit wet weight. In this section, a series of parametric vibration examples is set up focusing on the riser defined in Table 2 with damping ζ = 0.1 and internal flow velocity U = 2 m/s, and the time-domain vibrating solution is calculated through the DQM. After that, the parametric instability charts of the aforementioned model are obtained and depicted in Figure 11. be adjusted by changing the size of buoyancy tanks distributed along with the risers, a dimensionless coefficient kmw is used in our research to control the unit wet weight. In this section, a series of parametric vibration examples is set up focusing on the riser defined in Table 2 with damping ζ = 0.1 and internal flow velocity U = 2 m/s, and the timedomain vibrating solution is calculated through the DQM. After that, the parametric instability charts of the aforementioned model are obtained and depicted in Figure 11. With the solutions shown in Figure 11, it can be easily inferred that when the wetweight coefficient increases, the instability region of the riser's parametric vibration moves towards the bottom right. Thus, it can be concluded that the increase of wetweight coefficients will decrease the minimum dynamic tensions S, which means the parametric resonance will occur more easily. In addition, it can be inferred that longer excitation periods are needed to induce the risers' parametric resonance for models with larger wet-weight coefficients. Being different from the acceleration phenomenon of the influence of the internal flow velocity, the alteration of the unstable region is uniform with the increase of the wet-weight coefficient. This is because the coefficients αj, μ1 and, μ2 are linearly related to wet weight mw, while the U-dependent coefficient μ4 is quadratically related to the internal flow velocity.

The Dynamic Response of Different Excitation Periods
In this section, the responses of different cases are considered. Firstly, we discuss the responses of a riser with wet-weight coefficient kmw = 0.3 subjecting to parametric excitations of different periods, and the detailed parameters for these cases are listed in Table 3.
According to the instability chart shown in Figure 12, which is excerpted from Figure 11, the Cases I(d)-III(d) are located at the stable region of the riser with kmw = 0.3 (the cyan instability borderline), while the Cases IV(d)-VI(d) are located at the unstable region. Additionally, the response solutions of the first four modes of these cases are presented in Figure 13. With the solutions shown in Figure 11, it can be easily inferred that when the wetweight coefficient increases, the instability region of the riser's parametric vibration moves towards the bottom right. Thus, it can be concluded that the increase of wet-weight coefficients will decrease the minimum dynamic tensions S, which means the parametric resonance will occur more easily. In addition, it can be inferred that longer excitation periods are needed to induce the risers' parametric resonance for models with larger wetweight coefficients. Being different from the acceleration phenomenon of the influence of the internal flow velocity, the alteration of the unstable region is uniform with the increase of the wet-weight coefficient. This is because the coefficients α j , µ 1 and, µ 2 are linearly related to wet weight m w , while the U-dependent coefficient µ 4 is quadratically related to the internal flow velocity.

The Dynamic Response of Different Excitation Periods
In this section, the responses of different cases are considered. Firstly, we discuss the responses of a riser with wet-weight coefficient k mw = 0.3 subjecting to parametric excitations of different periods, and the detailed parameters for these cases are listed in Table 3. According to the instability chart shown in Figure 12, which is excerpted from Figure 11, the Cases I(d)-III(d) are located at the stable region of the riser with k mw = 0.3 (the cyan instability borderline), while the Cases IV(d)-VI(d) are located at the unstable region. Additionally, the response solutions of the first four modes of these cases are presented in Figure 13.       According to the response solutions shown in Figure 13, it can be informed that the solutions of Cases I(d)-III(d) are convergent while the responses of the other cases are divergent with respect to time. For the unstable cases (i.e., Cases IV(d)-VI(d)), the timedomain weight function of the second order q 2 fluctuates with a growing amplitude, while the fluctuating amplitudes of other modes are decreasing at the beginning and start to increase starting at a certain time. It also can be observed that the divergent speed of the second mode is much faster than other modes as it is shown that the vibrating amplitude of Case VI(d) at t = 300 s of the second reaches 0.05 m, while the vibrating amplitude at t = 300 s of other modes is no more than 0.005 m. This is because Cases IV(d)-VI(d) lie in the second order parametric unstable region. Thus, the instability of the second mode is stimulated by the parametric excitation initially. Whereafter, as the existence of the coupling effect between different modes, the vibration of other modes begins to diverge with the growth of the vibrating amplitude of the second mode. In Figure 13, we draw two transversal lines at t = t 1 and t = t 2 . Observing the time-domain responses between t 1 and t 2 , it can be informed that there is only one waveform between these two time-points in the solution of the first order, while there are three waveforms in the third order and four waveforms in the fourth order, regardless of the change of excitation period. This is because for the first, third and fourth modes, the free vibration takes a major part of the responses at the beginning before the coupling effect induced by the parametric resonance of the second mode takes effect. When observing the responses in the amplitude-increasing phase, it can be found that there exist slight phase differences between different cases, and these phase differences can be spotted with their accumulation with the increase of time. Moreover, when observing the responses between t = t 3 and t = t 4 , it can be seen that there is only one waveform for Case VI, regardless of which mode it is. The cause of these phenomena is that the free vibrations die down with the passing of time and the coupling effect induced by the second-order parametric resonance, whose period is related to the excitation period, gradually takes the major part. Thus, the response solutions of different modes vibrate with the same period as the second-order parametric resonance.

The Dynamic Response of Different Wet-Weight Coefficients
To study the influence of wet weight on risers' parametric vibration responses, we set a series of cases with different wet-weight coefficients, as listed in Table 4. The response solutions of these cases are solved with DQM and presented in Figure 14. In Figure 14, the responses of the first four orders of the six cases listed in Table 4 are illustrated. According to the response solution of the second order q 2 , it can be assumed that the vibrating amplitude of the second order responses of Cases III(a)-III(b) are increasing with the time, and those of Case III(c) are keeping stable, while those of Cases III(d)-III(f) are decreasing. The convergences of these six cases are consistent with the prediction of the instability chart Figure 12, that the mark of Case III lies in the unstable region of the instability regions of k mw = 0.0 and 0.1, and on the edge of the instability region of k mw = 0.2, but out of the instability regions of k mw = 0.3, 0.4, and 0.5.
For the first, third, and fourth orders of the modes, the response amplitudes decrease initially. Due to the coupling effect between modes, the divergence of the second-order excites the first and third mode to vibrate divergently from a certain time. At the decreasing phase of the responses where time t is between t 1 and t 2 , the vibrating period is altered in accordance with the wet-weight coefficient. For a bigger wet-weight coefficient, the free vibrating period increases. However, for the time between t 3 and t 4 , the free vibrations die down, and the parametric resonance takes the main part of the responses. At this phase, the vibrating periods of different modes reach an agreement. Additionally, these phenomena can tell that the alteration of the wet-weight coefficients will change the period of free vibration but have little effect on the period of parametric vibration. In Figure 14, the responses of the first four orders of the six cases listed in Table 4 are illustrated. According to the response solution of the second order q2, it can be assumed that the vibrating amplitude of the second order responses of Cases III(a)-III(b) are increasing with the time, and those of Case III(c) are keeping stable, while those of Cases III(d)-III(f) are decreasing. The convergences of these six cases are consistent with the prediction of the instability chart Figure 12, that the mark of Case III lies in the unstable region of the instability regions of kmw = 0.0 and 0.1, and on the edge of the instability region of kmw = 0.2, but out of the instability regions of kmw = 0.3, 0.4, and 0.5.
For the first, third, and fourth orders of the modes, the response amplitudes decrease initially. Due to the coupling effect between modes, the divergence of the secondorder excites the first and third mode to vibrate divergently from a certain time. At the decreasing phase of the responses where time t is between t1 and t2, the vibrating period is altered in accordance with the wet-weight coefficient. For a bigger wet-weight coefficient, the free vibrating period increases. However, for the time between t3 and t4, the

Conclusions
In the current work, the parametrical vibration problem of risers has been considered. The governing equation of parametrically excited risers is deduced. Moreover, the DQM is applied to Galerkin's reduced form of this equation to study the parametric vibration problem of the risers. By testing this solution scheme in two simple equations, the accuracy of the DQM method is evaluated. With the good consistency between the DQM and the R-K for the Mathieu equation, the competency of the DQM in risers' parametric vibration problem has been verified.
According to the numerical experiments on risers with different parameters, some important conclusions can be drawn as follows:

•
The instability charts are given by both DQM fitting well and with the borderlines given by the Floquet theory. • Increasing damping can control the parametric resonance from occurring.

•
The increase of internal flow velocity and the wet-weight coefficients will deteriorate the parametric instability of the riser. Additionally, the influences of internal flow velocity on risers' parametric instability are nonlinear, while the those of wet weight are linear.

•
According to the response solutions, the existence of the coupling term will give a chance to the parametrically unstable modes to excite other modes to vibrate unstably.

•
The alteration of the parametric excitation period will not change the period of free vibration, while the excitation period will have an effect on the parametrically vibrating period.

•
The alteration of the wet-weight coefficients will change the period of free vibration but have little effect on the period of parametric vibration.
Although some results can be obtained as referred above, there are still some limitations in the current study. Firstly, as the node distribution for the current solution scheme is nonuniform, it is to make the frequency domain analysis for the solutions. Secondly, since the block marching technique calculates the solution block by block, the error of this solution scheme would accumulate block by block. Thus, further study could be added to slow down the error accumulation.