Numerical Simulation of the Vortex-Induced Vibration of A Curved Flexible Riser in Shear Flow

A series of fully three-dimensional (3D) numerical simulations of flow past a free-to-oscillate curved flexible riser in shear flow were conducted at Reynolds number of 185–1015. The numerical results obtained by the two-way fluid–structure interaction (FSI) simulations are in good agreement with the experimental results reported in the earlier study. It is further found that the frequency transition is out of phase not only in the inline (IL) and crossflow (CF) directions but also along the span direction. The mode competition leads to the non-zero nodes of the rootmean- square (RMS) amplitude and the relatively chaotic trajectories. The fluid–structure interaction is to some extent reflected by the transverse velocity of the ambient fluid, which reaches the maximum value when the riser reaches the equilibrium position. Moreover, the local maximum transverse velocities occur at the peak CF amplitudes, and the values are relatively large when the vibration is in the resonance regions. The 3D vortex columns are shed nearly parallel to the axis of the curved flexible riser. As the local Reynolds number increases from 0 at the bottom of the riser to the maximum value at the top, the wake undergoes a transition from a two-dimensional structure to a 3D one. More irregular small-scale vortices appeared at the wake region of the riser, undergoing large amplitude responses.


Introduction
Long flexible risers are employed in the offshore industry to convey fluids from the seabed to sea level. Such slender circular cylinders are vulnerable to undergoing the vortex-induced vibration (VIV) in complex waves and currents. Fatigue damage of risers caused by the VIV may result in huge loss and environmental catastrophe. Therefore, detailed understanding and efficient prediction of the vibration response of risers have been increasingly urgent and hot issues for ocean engineers and researchers.
Considerable efforts have been dedicated to understanding the VIV of rigid cylinders over the past few decades. One may refer to the comprehensive reviews by Sarpkaya (2004), Williamson andGovardhan (2004, 2008) and Bearman (2011). Until the past two decades, the VIV of flexible risers with a large length-to-diameter (L/D) has been investigated experimentally by many researchers (Trim et al., 2005;Song et al., 2011;Xu et al., 2017). Trim et al. (2005) conducted VIV tests of a pipe model with an L/D of 1400 in a towing tank of the Norwegian Marine Technology Research Institute. In the similar experiments performed by Song et al. (2011), the aspect ratio was up to 1750 and multimode response of a bare cylinder was observed. More recently, the response feature of an inclined flexible bare cylinder with an L/D of 350 has been investigated by Xu et al. (2017) through a series of experiments in a towing tank. These experiments provide better insights into some VIV characteristics such as the dominant modes, response amplitudes and orbital trajectories versus the reduced velocity. Since the cylinder models were towed in the experiments, only uniform flow or linearly sheared flow was generated and considered. In our earlier experimental study (Zhu et al., 2016), the VIV of a flexible riser was tested in log-law shear flow generated in a water flume, which better mod-elled the real marine environment in the continental shelf area. Additionally, the riser model was deployed as the catenary type with a curved shape. The response amplitude, dominant frequency and wake mode of scattered markings were observed in the experiments. Though the information obtained by the experiments is limited, it provides some good benchmarks for verifying numerical models.
A number of numerical studies have been carried out on the VIV of flexible cylinders, in order to capture more detailed flow characteristics and vibration features. Multi-strip technique, as a quasi-three-dimensional method, was previously employed to simulate the riser response (Willden and Graham, 2001;Yamamot et al., 2004). The transverse vibration of a flexible cylinder with an aspect ratio of 100 was simulated by Willden and Graham (2001) in a linearly sheared flow at a low Reynolds number. Cellular shedding was observed and the vibration shown significant spanwise correlation over a large length of the cylinder. Yamamot et al. (2004) adopted the same method to simulate the oscillation of a flexible riser subjected to a uniform flow. In their simulations, a 2S (two single vortices are shed per cycle, similar to a von Karman vortex street) mode was found in the regions of low amplitudes, which changed to a 2P (two pairs of counter-rotating vortices are shed per cycle) mode in the regions of larger amplitudes. Since the cylinder is discretized in several panels in such multi-strip method, the three-dimensional (3D) vortex structure and instantaneous shape of the cylinder cannot be treated accurately. Thus, fully 3D numerical simulations have been carried out by several scholars. Kamble and Chen (2016) conducted 3D simulation of the VIV of flexible cylinders using computational fluid dynamics (CFD) approach with an overset grid system, and the results were in good agreement with the experimental data by Trim et al. (2005). Wang and Xiao (2016) carried out a numerical study on the vibration response of a flexible riser using the commercial software ANSYS CFX, and the numerical approach was validated with the experimental data by Lehn (2003). Further understandings of the VIV of vertical flexible cylinders have been obtained through such fully 3D simulations, but which are still quite limited. Moreover, few studies concerned about a curved flexible cylinder in a nonlinear shear flow.
In this paper, the VIV of a curved flexible riser model in log-law shear flow is investigated using a fully 3D fluid-structure interaction (FSI) numerical approach. The numerical results of the flexible riser under such a specific arrangement and flow conditions are validated with the experimental results reported in the earlier study (Zhu et al., 2016). Apart from the vibration responses, the 3D wake structure and instantaneous shape of the flexible riser are discussed in detail.

Problem description
The curved flexible riser model tested in our earlier study (Zhu et al., 2016) is considered in simulations. The flexible riser is fixed in two ends and is deployed as the catenary type in a vertical plane with its concave facing the incoming flow. The main parameters of the riser model are summarized in Table 1. As shown in Fig. 1, a cuboid computational domain is adopted to simulate the flow past such a curved riser. The height (H) equals the experimental water depth of 0.5 m. The width is 20D (D is the external diameter of the riser) and the two lateral boundaries are 10D away from and parallel to the xoz plane at which the riser model is located. The horizontal distance between the inlet and the lower end of the riser is 10D while the distance is 30D between the outlet and the upper end.
In the experiments, stable log-law shear flows with a mean velocity of 0-0.25 m/s were generated in the flume with the water depth fixed at 0.5 m, and the flow velocity profiles were obtained by an Acoustic Doppler Velocimetry (ADV). Among them, seven representative flow velocity profiles are adopted in simulations, as listed in Table 2. In the table, q is the flow rate per width, U is the free-stream velocity, and ( ) are the mean velocity and the mean reduced velocity, respectively, and Re is the corresponding Reynolds number. The first natural frequencies in the in-line (IL) and cross-flow (CF) directions are f 1x =2.43 Hz and f 1y =1.16 Hz, respectively, which were obtained through the decay tests. The second and third natural frequencies were calculated using the ANSYS modal, which are f 2x =5.13 Hz, f 2y =3.01 Hz, f 3x =9.24 Hz and f 3y =5.96 Hz. The upstream boundary is imposed with a velocity inlet boundary condition defined as the free-stream velocity profile. The downstream boundary is imposed with a pressure outlet boundary condition, which meets a linear static pressure distribution. Apart from the inlet and outlet boundary conditions, the opening boundary condition is applied on the top boundary as p=0 Pa, the symmetry condition is specified on the two lateral boundaries, and the no-slip boundary condition is assigned on the surface of the riser and the bottom wall.
As shown in Fig. 1, 31 black rings were marked along the flexible riser model in the experiments, which represent different locations along the span. The number of markings and the corresponding height (Z/H) are employed in the comparison between the numerical and experimental results in this paper.

Numerical methods
The vibration response of a flexible riser is a typical fluid-structure interaction issue. It is still a challenging task to conduct fully coupled calculation (He et al., 2012). Thus, the fluid and solid fields are usually solved separately, i.e., there are no iterations between the fluid and solid fields within one time step (Wang and Xiao, 2016). In the present work, a two-way FSI approach is utilized through the multifield solver of the software package ANSYS MFX, in which the CFX is employed to solve the fluid field while the transient structural module is adopted to calculate the structure response. The flow chart of the solution procedure is plotted in Fig. 2.
The unsteady Reynolds-Averaged Navier-Stokes (UR-ANS) equations coupled with a shear stress transport (SST) k-ω turbulence model (Menter, 1994) are used to capture the flow past the curved flexible riser. Additionally, the Arbitrary-Lagrangian-Eulerian (ALE) scheme (Hughes et al., 1981) is adopted to follow the moving boundary of the riser. Thus, the flow governing equations are expressed as: (1) where t is the time; p is the pressure; ρ is the fluid density; x i is the Cartesian coordinate in the i direction; u i and represent the instantaneous Cartesian velocity component and the fluctuation velocity component, respectively; u i or with a top bar represents its time average one; is the velocity of the mesh movement in the x j direction; μ and μ t are the dynamic viscosity and turbulent viscosity, respectively; k t is the turbulent energy and δ ij is the Kronecker delta function (δ ij = 1 for i = j). The SST k-ω turbulence model proposed by Menter (1994) has been demonstrated good performance in capturing turbulent vortices (Zhu and Yao, 2015), and is applied to closure the flow governing equations. One may refer to the literatures by Menter (1994) for more detailed information about this model. The flow governing equations are discretized with the finite volume method. The Rhie-Chow interpolation (Rhie and Chow, 1982) is used to deal with the coupling between the pressure and the flow velocity. A second-order backwind Euler scheme and a high resolution scheme are employed for the temporal discretization and the discretization of the convective terms, respectively.
In the software package ANSYS MFX, a finite element method is used and the riser is discretized into n elements. For each element, the governing equation is described as: where M is the mass matrix with the added mass being considered, C is the damping matrix, K is the equivalent stiffness matrix which is converted from the bending stiffness (EI), X i donates the displacement in the i direction and F i (t) is the hydrodynamic force in the i direction, which is computed by performing an integration involving both the pressure (p) and viscous stress (τ) on the micro element of the riser. The motion equations are solved by the Hilber-Hughes-Taylor method (Hilber et al., 1977).
In one time step, the hydrodynamic forces are firstly obtained by solving the flow governing equations. Then, the motion equations of the riser are solved to find a new position for the structure. Simultaneously, the mesh deformed following the motion of the riser, which is prepared for the flow field calculation at the next time step. The displacement diffusion model is employed for the mesh deformation, which meets , where S i is the displacement of the mesh points in the x i direction and k m is the mesh stiffness. In the simulation, the mesh stiffness is set as 10.
The time step Δt adopted for each case meets the criterion that the Courant number Co is kept below 1.0 in the entire computational domain, where Co is defined as: |u|

U∆t/D
where Δl is the size of a cell in the direction of the flow velocity and is the magnitude of the velocity through the cell. This results in a non-dimensional time-step ( ) of 0.001.
As shown in Fig. 3, a structured hexahedron mesh is constructed for the fluid domain. The regions around the riser model and near the bottom wall are tessellated with a fine mesh. In a xoy plane, the riser perimeter is equally dis-cretized with 116 nodes. The radial size of the first layer of mesh next to the riser is 0.01D (y + <5). The element expansion ratio in the O-ring region around the riser is kept below 1.2 while it is 1.0 for the rest regions. A close-up view of the mesh around the riser in the xoy plane is shown in Fig. 3c. Apart from the fluid domain, the solid domain is tessellated with a two-layer structured mesh, in which the riser perimeter and length are equally discretized with 24 nodes and 60 segments, respectively. The two ends of the riser are fixed, and the outer surface of the riser is set to be a fluid-solid interface for the data transfer.
A mesh independence test is conducted to ensure that the numerical results are independent of the grid size, i.e., the grid size is decreased until the results have no change. Six mesh systems with different sizes are used to simulate the case C12 and the results are compared in Fig. 4. The total number of the cells increases from 263712 (M1) to 1524808 (M6). The cross-flow vibration displacements at five different locations along the span direction are considered, and the locations are distinguished with the marking numbers used in the experiments. As the number of the cells increases, the difference of the numerical results be-  comes smaller, which reduces to only 0.98% between M5 and M6. The displacements predicted by M5 and M6 are in agreement with the experimental results (Zhu et al., 2016) though a small deviation presents at some locations. By taking into account the computational efforts, M5 is selected as the final mesh to carry out the simulations.

Results and discussion
4.1 Vibration response Fig. 5 compares the IL and CF root-mean-square (RMS) amplitudes between the present simulations and the earlier experiments (Zhu et al., 2016). The envelopes of the vibration amplitudes along the riser model at different Reynolds numbers illustrate that the dominant mode varies with the flow velocity. When the Reynolds number is smaller than 370, both the IL and CF oscillations are dominated by the first mode. The dominant mode changes to the second mode when the Reynolds number is in the range of 497-827. This phenomenon of the dominant mode increasing with the growth of the flow velocity is observed in many experimental studies on VIV of flexible cylinders, not only in the cases of uniform flows (Trim et al., 2005;Song et al., 2011) but also in the cases of linearly sheared flows (Huang et al., 2011;Gao et al., 2015). At Re=1015, the IL dominant mode shifts to the third one while the CF dominant mode keeps the second one. It illustrates that the dominant modes are out of phase in the streamwise and transverse directions. Such asynchronous phenomenon was observed by Wu et al. (2016) and Seyed-Aghazadeh and Modarres-Sadeghi (2016) in their studies on the vibration response of a straight flexible cylinder. In terms of the dominant modes and the envelopes of amplitudes, the simulation results are in good agreement with the experimental data. However, there are some deviations in the amplitudes along the riser model. The locations presenting the peak and trough exhibit relatively large deviations. Both the envelopes of the IL and CF amplitudes are asymmetry along the span direction, i.e., the peak point of the first-mode dominated vibration and the trough point of the second-mode dominated vibration do not locate at the half of the water depth. The nonlinear flow velocity profile and the nonlinear deployment of the riser contribute to this result. For the CF amplitude, it is clearly seen that the trough point of the second-mode dominated vibration, as a node exhibiting the minima of the response envelope, locates at Marking #17 (Z/H = 0.37), which is just the midpoint of the riser model. It indicates that the midspan of a flexible cylinder is a characteristic point for vibration response. Additionally, the CF RMS amplitude of the node does not equal zero, illustrating that more than one mode participate in the vibration and the second mode is the dominant one. It is confirmed in Fig. 6 that the time histories of the vibration amplitudes exhibit non-constant values since multiple frequencies take part in the dynamic response. The contributions of different frequencies are variable along the span. This mode competition phenomenon was observed by Gao et al. (2016) for a straight flexible cylinder in uniform and linearly sheared flows, and by Han et al. (2017) for an inclined flexible cylinder in uniform flow. However, the node of the IL RMS amplitude is not located at the midspan of the riser, which is different from the straight flexible cylinders. For this curved flexible riser, the IL motion belongs to an in-plane vibration while the CF motion is an out-ofplane vibration. Therefore, there are some different characteristics in the vibration responses between the two different directions, which are reported by Meng and Chen (2012) and Quéau et al. (2013).
Apart from the RMS amplitudes, the orbital trajectories of the riser model at seven typical positions are also compared in Fig. 5. The trajectories are plotted in the xoy plane with the same scales in the x and y axes. As seen from Fig.  5, some of the trajectories exhibit eight-shape figures indicating the 2:1 ratio on the IL to CF excited frequencies. However, C-shape and line-shape trajectories appear at some positions where the IL amplitudes are much smaller than the CF amplitudes. These trajectories mainly occur at small flow velocities or at the positions near the two ends of the riser model. For the cases dominated by the second mode, some of the trajectories become a bit chaotic, which is attributed to the frequency competition.
The vibration frequencies of the riser model at seven typical positions are shown in Fig. 7. In the figure, the frequencies are normalized with the first natural frequencies (f 1x and f 1y for the IL and CF directions, respectively), and f s is the Strouhal vortex shedding frequency defined as: where St is the Strouhal number taken as 0.18. The oblique dashed lines in Fig. 7 describe the normalized Strouhal vortex shedding frequency which increases linearly with the reduced velocity. U r It is clearly seen that only one vibration frequency appears at low reduced velocities ( ≤7.921) while more than one frequency take part in the vibration when the reduced velocity is larger than 8.853. There are two frequencies participating in the CF oscillation for different positions along the riser in the range of Re=370-1015. It is interesting to find the occurrence of the frequency competition at Re=370 (C6), which is dominated by the first mode. The reason is that it undergoes the mode transition from the first mode to the second one. When the vibration is domin- ated by the second mode (the dominant frequency becomes higher), the low frequency corresponding to the first mode still appears although its contribution is lower. Additionally, the frequencies at different positions do not keep the same value in the mode-transition cases such as C6, C8 and C16. It indicates that the transition of the vibration frequencies is out of phase along the span direction. Though some part of the riser keeps the dominant frequency of the first mode (such as the midspan in C8), the dominant mode of the riser shifts to the second one when most parts vibrate in a dominant frequency close to the second natural frequency.
The most obvious difference between the IL and CF vibration frequencies is that there are three frequencies participating in the IL vibration at Re=1015 (C16), which is dominated by the third mode. It further illustrates that the lower modes still take part in the vibration though it is dominated by a higher mode. Besides, the difference in the IL frequencies at different locations is more obvious than that in the CF direction, indicating that the frequency competition is more intense in the IL direction. Fig. 8 depicts the instantaneous shapes of the flexible riser and its response envelopes at different flow velocities. It is clearly seen that the instantaneous shape of the riser model varies over time with different amplitudes along the span. Since the IL displacements are much smaller compared with the horizontal span, the instantaneous shape of the riser presented on the left column of the figure is seen from the upstream along the flow direction, which mainly shows the CF displacements. In the figure, t 0 is the time when the riser reaches its equilibrium position and t 0 -T/4 and t 0 +T/4 are two other typical moments during a vibration cycle, where T is the CF vibration cycle. Besides, the surface of the riser is colored by the transverse velocity (v) of the fluid around the riser, which may reflect the fluid-structure interaction. It is seen from Fig. 8 that the  Fig. 8. Instantaneous shapes of the flexible riser seen from the upstream (along the flow direction) and the response envelopes of the riser (the IL displacements of the envelopes on the right column are magnified by a factor of 120). transverse velocity of the ambient fluid reaches its maximum value when the riser reaches the equilibrium position. For the first-mode dominated cases (C5), the red surface around the middle part of the riser indicates that the ambient fluid has the maximum velocity along the positive y axis. For the second-mode dominated cases, both the red surface and the dark blue surface of the riser illustrate that the ambient fluid reach local maximum velocities, though one is along the positive axis and the other one is along the negative axis. The local maximum transverse velocities of the ambient fluid occur at the peak CF vibration amplitudes. Additionally, these local maximum transverse velocities have relatively large values when the vibration is in the first and second resonance regions (C5 and C12, respectively). It implies that the transverse velocity of the ambient fluid is mainly induced by the CF vibration of the riser. When the riser reaches the place where the CF displacement reaches the largest value (t 0 ±T/4), the transverse velocity of the ambient fluid becomes smaller.

Instantaneous shape of the riser
As shown in Fig. 8, the response envelopes more clearly present the dominant modes both in the IL and CF directions. When the vibration is dominated by a high mode, it is seen that the nodes of the envelopes deviate from the original position of the riser model, which is more obvious in the IL direction. It further confirms that multiple frequencies take part in the response and the frequency competition in the IL direction is more intense.
4.3 Flow field and wake structure Fig. 9 shows the flow velocity and vorticity around the flexible riser in one vibration cycle at Re=827. In the figure, the bold black line represents the riser, whose location varies over time. It is seen that the vortices are shed alternatively from the two sides of the riser, and the streamwise flow velocities in the core of the vortices are negative. The minimum flow velocity is -0.03 m/s, appearing in the core of the vortex just shed from the riser surface. In addition, the velocity contours at different levels reflect the shear flow profile. The maximum flow velocity occurs at two sides of the riser near the water surface. However, the velocity difference is not obvious at the lowest plane close to the bottom, which is attributed to the relatively small incoming flow velocity. Therefore, the corresponding vortices in the wake region are also relatively small.
The flow wakes at six different Reynolds numbers are compared in Fig.10. The 2S, 2P and P+S (a pair of the vortices are shed on one side and a single vortex is shed on the other side) vortex shedding mode are captured in the simulations. As seen from the contours of ω z at different xoy planes, different vortex shedding patterns are appeared at different positions along the span, where ω z is defined as: |ω z | where u and v are the IL and CF flow velocities, respectively. In the first-mode dominated cases, the majority of the vortex shedding exhibits a clear 2S pattern, whereas a P+S pattern is observed in the wake of the midspan where the maximum amplitude appears. In the second-mode dominated cases, a P+S mode appears at the positions corresponding to the peak amplitudes while a 2S mode mainly presents at the position corresponding to the trough. Such observations are in agreement with the conclusions of , Meneghini et al. (2004) and Wang and Xiao (2016) for the vibration response of flexible risers that a 2S mode appeared behind the parts of small amplitudes. Moreover, the asymmetric P+S mode at most positions along the riser is attributed to the multi-frequency responses. Apart from the wake structure, the numerical results provide some quantitative information about the vortex intensity and size. As shown in Fig. 10, becomes larger as the flow velocity increases, indicating the riser extracts more energy from the fluid flowing in a higher velocity. Additionally, the vortex size is relatively large behind the parts of the riser undergoing large amplitudes. It implies that the motion of the riser, in turn, affects the vortex shedding. Therefore, the numerical simulations give some complements for further understanding the experimental results.  In addition, the 3D wake structure that is hard to capture in experiments is obtained by the simulations, as shown in Fig. 11. The left column of the figure exhibits the iso-surfaces of , while the right column presents the iso-surfaces of (streamwise vorticity) which is defined as: where w is the vertical flow velocity. As seen from Fig. 11, the vortex columns are shed nearly parallel to the axis of the riser. It indicates that the vortex shedding is related to the yaw angle of the riser. This phenomenon was observed by Miliou et al. (2007) in their simulations of flow past a curved cylinder and by Assi et al. (2014) and Seyed-Aghazadeh et al. (2015) in their experiments of flow over a fixed concave cylinder. Besides, the wake of the lower part of the riser near the bottom wall is well organized with no obvious shedding, while the wake flow of the remainder part is entirely 3D. The nonlinear shear flow profile contributes to this variation, since the depth-dependent incoming |ωz| |ω x | Fig. 11. Iso-surfaces of and for various Re values (the flow is from the left to the right), in which the red and blue colors denote positive and negative values of ω z , respectively, and the green and yellow colors denote positive and negative values of ω x respectively. flow velocity decreases quickly as Z/H is reduced. From the bottom of the riser to its top, the local Reynolds number increases from 0 to the maximum value. Therefore, the wake structure undergoes a transition from 2D to 3D. It is similar to the wake transition of a straight circular cylinder as the incoming flow velocity increases (Williamson, 1992;Jiang et al., 2016). The wake at Re=827 (C12) is dominated by more irregular small-scale vortices, which agrees with the observation of Zhang et al. (2017) for a straight circular cylinder undergoing the VIV in the upper branch. It indicates that the riser experiencing a strong vibration obviously affects its 3D vortex structure.

Conclusions
A series of fully 3D numerical simulations have been carried out to investigate the VIV of a curved flexible riser subject to shear flow. The two-way FSI simulations are accomplished by using the multi-field solver of the software package ANSYS MFX. The numerical results are in good agreement with the experimental results, indicating that the present numerical method is reliable and capable of predicting vibration responses of nonlinear-deployed flexible cylinders in nonlinear shear flow. The main findings of the present work are summarized as follows.
(1) The dominant mode increases with the growth of the flow velocity, but it is not synchronous in the streamwise and transverse directions. Moreover, the transition of the vibration frequencies is out of phase along the span direction. The lower modes still participate in the vibration though it is dominated by a higher mode. Such mode or frequency competition results in the non-zero nodes of the RMS amplitude and the relatively chaotic motion trajectories of the riser.
(2) The transverse velocity of the ambient fluid reaches the maximum value when the riser reaches the equilibrium position, while it becomes smaller when the riser leaves its equilibrium position. Moreover, the local maximum transverse velocities occur at the peak CF amplitudes, and the values are relatively large when the vibration is in the resonance regions. It reveals that the transverse velocity of the ambient fluid is mainly induced by the CF oscillation. (3) Apart from the wake structure, the numerical results provide quantitative information about the vortex size and the vorticity. As the flow velocity increases, becomes larger, illustrating that the riser extracts more energy from the flow with a higher velocity. In addition, the vortex size is relatively large for the parts undergoing large amplitudes, indicating that the motion of the riser, in turn, also affects the vortex shedding.
(4) The vortex columns are shed nearly parallel to the axis of the curved flexible riser. The wake structure undergoes a transition from 2D to 3D as the local Reynolds number increases from 0 at the bottom of the riser to the maximum value at the top. For the riser undergoing VIV in the lock-in region, the wake is dominated by more irregular small-scale vortices.