Marine risers are important equipment for the exploration of oil/gas in deep water, which usually serve as the bridge between an offshore platform and the well head on the seabed (Kaewunruen et al., 2005; Ju et al., 2012). As exploration activities move into deepwater, the long slender risers tend to undergo large-amplitude motions (Xu et al., 2013). Internal waves are believed to be responsible for a great deal of damage. This is because they can create enormous local loads and bending moments on offshore structures. The internal waves have been reported to induce an additional displacement of 200m in the horizontal plane and 10m in the vertical direction (Chakrabarti, 2005). Therefore, the internal waves should be considered in analysis and design of marine risers.

The internal waves occur frequently in South China Sea where water is deep (Du et al., 2001). At the current stage, most investigations focus on the generation, propagation, transformation and other properties of internal waves, few studies are related to their effects on offshore structures. Cai et al. (2003; 2008) firstly introduced Morison empirical formula and modal separation method to estimate the forces and torques exerted by internal solitons on cylindrical piles. They found that the background currents will modify the properties of internal waves and enlarge both force and torque (Cai et al., 2006; 2008). Ye and Shen (2005) calculated and analyzed the force of internal waves on small-scale cylinder in different frequencies and then compared with that of surface waves and currents. Du et al. (2007) found that the maximum total force caused by a soliton of a current with speed of 2.1 m·s^-1 is nearly equal to that by a surface wave with a wave length of 300 m and a wave height of 18 m. A lot of studies have been conducted on the effects of internal solitons on marine risers in time domain by experimental, numerical and even CFD methods (Song et al., 2011; Chen et al., 2011; Liu et al., 2011). Jiang et al. (2012) estimated the dynamic characteristics of a Top Tension Riser (TTR) under internal solitary and non-uniform current by mKDV theory in a two-layer fluid. Guo et al. (2013) computed the dimensionless displacement and stress of TTR under combined excitation of internal solitary waves, surface waves and vessel motion in time domain. However, most of these investigations are concentrated on the nonlinear internal waves such as internal solitons and internal solitary waves by using KDV theory, few are on the dynamic response of marine riser under combined excitation of internal waves of high frequency and background currents.

In this study, the finite element method is applied to calculate the dynamic response of a marine riser conveying fluid subjected to internal waves and background currents. The basic equations of internal waves should be modified if background currents are considered. Thus, it is instructive to demonstrate the effect of background currents on the dynamic response exerted by internal waves on a supposed riser for comparison with no current case. The schematic diagram for a marine riser conveying fluid in internal waves and background currents is shown in Fig. 1. In sections 2 and 3, the governing equation of internal waves, the differential motion equation of a marine riser and the solution methods are presented. The computational results of different cases are discussed in section 4 and the conclusion is given in section 5.

With the Boussinesq and linear approximation, the second-order differential equations govern the vertical structure function W(z) of internal waves in a 2-dimensional Cartesian X, Z coordinate system satisfies the following equation:

where k_v is vertical wave number, z is the vertical coordinate.

In the case of internal waves without background currents, the vertical wave number k_v in Eq. (1) takes the form:

where k_h is the horizontal wave number of internal waves, ω is internal waves frequency, ω_i is inertial frequency and N is the Brunt-Väisälä or buoyancy frequency which is related to the density gradient via:

Consider the dimensional density , which is assumed as Holmboe model, given by (Ye and Shen, 2004):

where is the density of sea water in z₀, z₀ is vertical coordinate of the middle of pycnocline, h is half the thickness of pycnocline, α , β are dimensionless coefficient correlating with the density. The range of β value is [2, 4], here taking 3.5.

However, if the background currents exist, the modified Taylor-Goldstein equation for linear internal waves must be used. Fang et al. (2000) generalize to allow for a background current U₀:

where U₀ is the background current along the progressing direction and vertical curvature is the second derivative of z.

Eqs. (2) and (4) together with Eq. (1) are introduced to calculate the wave function for the case with or without background currents, respectively. Nevertheless, since it is difficult to express N(z) and U(z) analytically, the dispersion relation and wave function are usually decided with numerical method. Herein imposing rigid-lid boundary conditions at the surface, z=0, and at the bottom, z=-H. Then the Thompson-Haskell method (Fliegel and Hunkins, 1975) was used to solve governing equation of internal waves, where the water column is divided into many continuous levels, and if the vertical discrete level is thin enough, the buoyancy frequency N(z) could be regarded as a constant so that the analytic solution of Eq. (1) can be obtained. Finally, the eigenvalues k_v and their functions can be solved numerically by a joint equation transformed with a complete set of matrix equation.

In the two-dimensional Cartesian x, z coordinate system, the water pulsation can be assumed as wave solution and displayed in the form (Ye and Shen, 2005):

where W(z) is the normalized dimensionless amplitude of vertical velocity distributed along the depth in ocean and A is the maximum amplitude of vertical velocity.

According to the relation between the vertical displacement of water particles η and vertical velocity w, one obtains:

Note that A=η_maxω, it means that once the maximum wave amplitude η_max is determined, the vertical velocity profile is obtained.

Based on the continuity equation for incompressible flow, the horizontal relative velocity u_h and acceleration of internal waves are obtained with some manipulation:

where the small dot ( ⋅ ) represents the derivative with respect to time t.

The following assumptions are stipulated in the present mathematical model: the material property of the marine riser is homogeneous and linearly elastic and the cross-section are uniform and remain plane perpendicular to the axis at all state; longitudinal strain is large, but the shear strain, torsional deformation and rotational inertia are neglected; the internal fluid is inviscid, incompressible, irrotational, flowing upward and the density is uniform along the riser.

The nonlinear free vibration is determined through the virtual oscillations along the marine riser. The riser will move nonlinearly from the static equilibrium configuration due to small perturbation, as shown in Fig. 1. In the Cartesian coordinates, the normal displacement x=x(z,t) is considered to define the deformation of any point on the neutral axis of the static equilibrium configuration that was obtained from the static analysis.

According to mathematical model and strain-displacement relationship, when the effect of shear deformation is considered negligible, the curvature of the riser and the axial strain can be expressed as follows:

where s is the arc-length of riser element, the superscript ( ' ) represents the derivative of the parameter with respect to z, note that ; ε₀ is the initial axial strain due to the top tension.

Total strain energy

The total strain energy of a marine riser, which consists of the strain energy due to axial deformation and bending, is given by (Chucheepsakul et al., 2003):

where E is the Young’s modulus, A_r is the area of the cross-section, the initial stress σ_ov = 2υ (P_eA_re − P_iA_ri) / Ar , and υ is the Poisson’s ratio, P is the hydrostatic pressure, the additional subscripts e and i refer to the exterior and interior of the riser. I is the moment of inertia of the cross-section.

By substituting Eqs. (10) and (11) into Eq. (12), the variational expression of the total strain energy is written as:

where the effective tension T_e = EA_rε + 2υ[P_eA_re sgn(H -z) − P_iA_ri] and the sign function is defined as .

Virtual work done by external forces

Needed for dynamic response analysis, the virtual work done due to external forces included inertial force, structural resistance and the hydrodynamic force are taken into account for accurate reckons, i.e., δW =δW_I +δW_c +δW_f . Based on the newton’s second law, the virtual work done due to inertial force is given by (Kaewunruen et al., 2005):

where m represent the mass per unit length, the subscripts r, i, a stand for the riser, the internal flow and the added mass, respectively; V is the internal flow velocity,

The virtual work done by the structural resistance is written as:

here c is the structural resistance coefficient.

Following is the virtual work done due to hydrodynamic forces that are decomposed into two directions of the water flow velocities. The horizontal direction is considered only and the forces on risers are calculated by Morison formula. Then the work done by hydrodynamic loading is given by:

where

In Eq. (17), F_l is hydrodynamic force due to the internal waves, is the mean density of the static seawater, C_D is the drag force coefficient, C_M is the inertia coefficient, D is the outer diameter of the marine riser.

Based on the principle of virtual work, the total virtual work of the apparent system is zero:

By substituting Eqs. (13)~(16) into (18) and integrating by part, one obtains the total virtual work energy of the riser system expressed as below:

where m = m_r + m_i + m_a , and the superscript (4) represents the fourth derivative of displacement u .

Since δu is arbitrary, the differential motion equation of the marine riser in the horizontal direction can be expressed as below:

While both ends of the marine riser are assumed to be hinged, the boundary conditions are given by:

Discretize Eq. (21) by Hermite interpolation function, obtaining the general form of the finite element method:

where [M], [C], [K] and [F] represent the total mass matrix, the total damping matrix, the total stiffness matrix, and the external force matrix, respectively:

where [N] is the Hermite interpolation function, l is length of finite element of the riser.

Newmark-β method

To determine the dynamic responses of the riser, the numerical direct integration method is employed in this study. Newmark integration method is found suitable for this problem. The integration constants used are

The initial conditions for determining the dynamic responses of the marine riser are:

Programs Fre.m and resp.m in MATLAB are developed to solve the natural vibration frequency and the dynamic response of riser in consideration of background currents only. The results of ANSYS analysis which are taking the pipe59 element to simulate and solved by the prestress method are chosen to verify the accuracy of the programs. The parameters of riser and fluid, which are applied in this work, are given in Tables 1 and 2 respectively. The background currents are chosen from the one-year return currents profile in the Liuhua oilfield of South China Sea, as shown in Fig. 2.

The results in the first six natural frequencies of marine riser are illustrated in Table 3. In addition, the observational background currents and its mean velocity are employed to discuss the effect on the marine riser, as shown in Figs. 3 and 4, respectively. At first, it is found that the riser has a static offset in the uniform background currents. However, the maximum response displacement changes in some vary or degree when the background currents profile are different. The reversal point moves upwards as the background currents concentrate on upper layer, while locates on the center of riser in the mean currents. The Table 3 and Figs. 3 and 4 show a good agreement between the program results and ANSYS software analysis.

The parameters of Holmboe density distribution in ocean are given in Table 4 (Ye and Shen, 2005). Then, the density distribution profile and the computed buoyancy frequency N are obtained and as shown in Fig. 5. According to vast field observations and investigations of the internal waves in the South China Sea in recent decades (Bole et al., 1994; Qiu et al., 1996; Fang et al., 2000; Li et al., 2003), the internal waves of high frequency (about 3~5 cph) and huge amplitude (about 15~100 m) are commonly observed. In this study, a typical internal wave with frequency of 5.2358×10^-3 rad/s and maximum amplitude of 30 m is selected to study the dynamic response.

By solving the governing differential equation of internal wave using the Thompson-Haskell method with a vertical resolution of 1m, the distribution of waves eigenfunction for the first three modes are obtained. The horizontal wave number k_h are 0.00378 m^-1, 0.01553 m^-1 and 0.02628 m^-1, respectively. Meanwhile, the wave amplitudes of 2^nd and 3^rd modes are assumed as 1/2 and 1/3 that of 1^st mode, thus, the horizontal velocity distribution of internal waves versus depth are computed, as shown in Fig. 6. This figure also displays the transfiguration of a riser at different time in different modes. It can be seen from Fig. 6 that the configurations of riser act out in different shape. The reversal point locates on upper layer in the 1^st mode while shifts downward to the mixing layer in higher mode.

Moreover, Fig. 7 illustrates the time history of response displacement at the reversal point under internal wave alone corresponding to the first three modes. It shows that the oscillating cycles of riser are same as the wave period, and the maximum response displacements are in the same magnitude at the first two modes while rapidly diminish at the third mode. Hence, it suggests that the first two modes of internal waves play dominated role in the dynamic response of marine riser.

To study the dynamic response of riser under combined excitation of internal waves and background currents, the following five cases are considered as the follows. It should be noted that the main data of internal wave are kept the same as the previous section and the main parameters of uniform currents are taken form Fig. 2

(1) Mean currents, i.e., the velocity U₀(z) = 0.3757 m/s is constant.

(2) Observational currents, current speed U₀(z) is kept the same as case 1, but the vertical curvature(z) are considered in calculation.

(3) Observational currents, mean currents speed U₀(z) are used to solve the governing equation of internal waves and (z) are in consideration, while the observational data is applied in computing the hydrodynamic force.

(4) Observational currents, i.e., all U₀(z) and (z) are taken from the observational data in Fig. 2.

(5) Same as case 1 but the currents are in the opposite direction of internal waves.

For all the shear current cases, the local Richardson number Ri( z )= (N( z )/(z))² for wave breaking is also computed, and the results Ri( z ) is always >1/4, which suggests that the flow is always stable.

The horizontal wave numbers k_h are calculated by Eqs. (1) and (5), and Table 5 are filled with the results corresponding to the first modes. It is found that the wave numbers are smaller in case 1-4 compared with the case that no background currents exist, while greatly increase in case 5. This phenomenon means that the background currents make waves going along them stretched and wave against them shrunk.

Fig. 8(1) shows the solution in case 1. It is found that the variations in the horizontal relative velocities of internal wave can be ignored. Compared with Figs. 3 and 6, the transfigurations of riser vary greatly. The maximum displacement of riser x_max up to 0.86 m in the first mode and have a certain increase in higher modes. This figure demonstrates that when background currents exist, the displacements of riser exhibit a significant growth, especially in the first mode of internal waves. Furthermore, one can see that the maximum displacement of the riser varies versus time.

When the vertical curvature ( z )is incorporated, the curves of horizontal relative velocities profile of internal waves are not smooth in the first three modes, as shown in Fig. 8(2). Compared with case 1, the only difference is the slight variation in the first mode, where the maximum displacement of the riser is increased by 0.27 m. No differences are observed in higher modes. In fact, the same consequence can be obtained by the wave numbers of case 1 and 2 in Table 5. It means that the vertical curvature (z) has certain influence on dynamic response of a riser under combined excitation of internal waves and background currents and should not be ignored in the first mode.

In case 3, mean currents are used to solve the governing equation of internal waves, however, the observational data are served to evaluate hydrodynamic forces on marine riser. As shown in Fig. 8(3), it shows that the maximum displacements x_max = 1.62 m. Compared with case 1 and 2, the maximum response displacement has a great increase in the first mode but change little in higher mode. In addition, it is found that the transfigurations of the riser are different from that in mean currents and internal waves. The reversal point move upwards for the shear flow is mainly concentrating on the upper layer.

When the observational data is applied to solve the governing equation of internal waves, the eigenfunction (i.e., the vertical wave structure) will break in high mode and only be obtained in the first mode. Fig. 8(4) marks the horizontal relative velocities of internal waves versus depth in the ocean and its effect on the dynamic response of riser. It can be seen that the curve of velocity is rougher than case 2. The maximum displacement reduces by 0.08 m, but the riser has a similar shape. It can be concluded that the method used in case 3 is best reasonable and applicable.

Fig. 8(5) displays the horizontal relative velocity profile of the internal waves and the displacements along the riser at different times. In this case, the waves go against the uniform currents. It is found that the relative speeds of internal waves are lower than the other cases. Meanwhile, the maximum response displacement x_max = 0.76 m is slightly less than other cases but still lager than that either internal waves or background currents did alone.

As mentioned above in 5 cases, it is found that the internal waves will exhibit significant effects on riser motions in the 1^st mode, when background currents exist. In order to observe and compare the changes in different cases more conveniently, herein, the velocity profile of internal waves and maximum transfigurations of riser corresponding to 1^st mode are plotted on a graph, respectively, as shown in Figs. 9 and 10. It can be seen from Fig. 9 that the maximum horizontal relative velocity profile of internal waves have fewer changes on the whole in different cases except the case that internal waves propagate in opposite direction to currents, while the extreme values, which are 0.50 m/s, 0.64 m/s and 0.63 m/s in case 1, case 2&3 and case 4 respectively, have significant changes partially. Meanwhile, the horizontal velocities tend to be unsteady at some positions when the vertical variations of currents are incorporated. It means that both the distributions of background currents U₀(z) and the vertical curvatures (z) have significant influence on the dynamic response of a marine riser under internal wave excitations, as described in Fig. 10. The displacement increases greatly compared with case 1, while slightly lower than case 3.

In this paper, the dynamic performance of a marine riser under combined excitation of internal waves and background currents are presented. The dispersion relations and relative velocities of internal waves in horizontal plane are obtained, by solving the geometrically nonlinear equation for riser motion. The numerical solutions about displacement and transfiguration of riser with elapsing time are analyzed. The following conclusions have been drawn:

(1) Without background currents, the internal waves in the first two modes play a dominated role in dynamic response of a marine riser. The displacement of the riser diminishes rapidly in higher modes.

(2) The background currents will make waves that go along them get stretched and waves against them shrunk, which is in agreement with other investigations (Fang et al., 2000). Moreover, the displacements of the riser caused by internal waves in the first mode have an enormous increase when background currents exist. In this case, the distribution of both the currents speed U₀(z) and the vertical curvatures (z) make great contribution.

(3) When background currents exist, the method that the mean value of observational currents speed are used to calculate the internal waves velocity field, while actual value is applied to evaluate the hydrodynamic force is reasonable and applicable in dynamic response analysis of riser, and it is simple and convenient for application in engineering. Meanwhile, the vertical curvature should not be neglect.

(4) The displacements of the riser enlarge greatly, no matter whether the internal waves propagate along or against the direction of background.