Development of a Time-Domain Simulation Code for Cross-Flow Vortex-Induced Vibrations of a Slender Structure with Current Using the Synchronization Model

Abstract

In this paper, the numerical analysis method for the cross-flow vortex-induced vibration (CF VIV) analysis based on the proposed procedure for CF VIV analysis of slender structures is developed. In order to consider the changes in the incoming flow according to the static configuration of the slender structures due to the current, the proposed procedure has three stages. A slender structure is modeled as the lumped-mass line, and the dynamic relaxation method known as the numerical technique for a slender structure with large geometric nonlinearity is applied in the static analysis. The comparison studies with a commercial program are carried out to validate the developed code. The vortex-induced force on slender structures is considered with the synchronization model. To verify the developed CF VIV analysis procedure and numerical method for a slender structure, VIV analysis of the tensioned flexible risers under a uniform and shear current is performed. The simulated results of CF RMS displacement show good agreement with the results of the model test It is found that a tensioned riser vibrates with one dominant frequency in resonance with the nth mode, even though multi-frequencies components of the vortex shedding along the riser due to the shear current occurs.

1. Introduction

Vortex-induced vibration (VIV) is an important phenomenon that occurs around slender structures with the current flows. VIV is occurred due the fluctuating pressure related to the vortex shedding process. Despite the steady incoming flow, the vibrations of slender offshore structures are caused by vortex-induced vibration and lead to a shortened fatigue life due to the accumulation of fatigue damage. Therefore, understanding vortex-induced vibration is essential for the design of slender offshore structures. Extensive research has been carried out for understanding vortex-induced vibration [1,2,3,4,5,6,7,8]. Based on the understanding, various methods for predicting vortex-induced vibration have been developed and applied for design.

The VIV prediction methods are classified into semi-empirical models [9,10,11,12,13,14,15,16,17], wake oscillator models [18,19,20], and computational fluid dynamics (CFD) methods [21,22,23,24]. In the case of wake oscillator models, the wake behind a slender structure is modeled using a Van der Pol oscillator. Simulated results using the wake oscillator models show a qualitative agreement with the model test results. However, it is difficult to find parameters of the wake oscillator models that satisfy both the free and forced vibration phenomenon. CFD methods provide promising results in agreement with model tests due to advanced computer performance, but demand a lot of computational resources for practical calculations. Finally, semi-empirical models are classified into frequency domain and time domain methods. Semi-empirical frequency domain models using databases of hydrodynamic forces from the various model tests are commonly used, but the semi-empirical frequency domain models, which are linear, cannot explain the interaction between different response frequencies. Semi-empirical time-domain models give rise to more possibilities than frequency-domain models because there is no need to linearize components.

The synchronization model, one of the time domain models, was proposed by Thorsen et al. [25]. The synchronization model is simple but accurately mimics the lock-in phenomenon of VIV. The synchronization model predicts the VIV loads using the phase relationship between the VIV loads and the velocity of the oscillating body. The synchronization model, which consists of lift, drag, and added mass term, only predicts VIV loads acting on the cross-flow direction of incoming flow in the cross-section. To predict the hydrodynamic loads other than the cross-flow direction of incoming flow, the Morison equation is used. Since the Morison equation and the synchronization model have physically the same components (drag and added mass) special attention is required when applying both models simultaneously. To apply both models simultaneously, Thorsen et al. [26] propose a modified synchronization model which uses a fixed lift drag coefficient. The advantage of the modified synchronization method [26] can estimate VIV loads when the velocity of incoming flow irregularly changes over time. However, it is difficult to determine the appropriate fixed lift and drag coefficient.

Static configurations of the slender body are deformed under the influence of the current. When the static configuration of the slender body changes, the velocity of the incoming flow at the cross-section of the slender body changes. Therefore, it is important to calculate the exact static configuration of the slender body with the current for predicting accurate VIV loads. Since there are also slender bodies with large geometrical nonlinearities such as SLWR (steel lazy wave riser), a robust prediction method is required to calculate the static configuration of slender bodies with the current. Since it is important to calculate the exact static configuration of SLWR in the basic design, several studies [27,28,29] related to the static configuration of the riser have been carried out. The method of estimating the static configuration with the current is also studied by Wang and Duan [30].

In this paper, the numerical code for CF VIV analysis of a slender structure with the current is developed based on the proposed CF VIV analysis procedure applying the synchronization model for vortex-induced force. Two synchronization models for the CF vortex-induced force on a circular body are observed in the numerical code of one degree of freedom and verified from the comparative study with model tests. After validating the CF vortex-induced force model, a nonlinear numerical analysis model in time history is developed in order to perform CF VIV analysis for long slender structures. The developed numerical code for CF VIV of the slender body reflects the effective incoming flow at each cross-section according to the change in the static configuration with the current. The proposed procedure consists of three steps. The first step is calculating the static equilibrium condition of the slender body due to the current. The dynamic relaxation method, numerically robust and efficient, is used for the estimation of the nonlinear static equilibrium of a slender body with the current. The second step is calculating the effective incoming flow velocity from the static configuration of a slender body with the current. The final step is the dynamic analysis of the CF VIV of the slender body due to the current. CF VIV loads from the synchronization model [25], using effective incoming flow velocity, are considered for the dynamic analysis of CF VIV of the slender body. In this study, unlike the simplified equation of Thorsen et al. [26], a method of eliminating Morrison drag in the CF direction was used to combine with the Morrison equation and the synchronization model [25]. There is no need to simplify the synchronization model [25], and there is no difficulty in selecting a specific lift and drag coefficient. An efficient lumped-mass line model for static and dynamic analysis of the slender body is used in numerically stable and efficient. For the time integration of the dynamic analysis, the stable and accurate fourth-order Runge–Kutta method is applied. To validate the lumped-mass model for the static and dynamic analysis of the slender body, comparisons with commercial code, OrcaFlex 10.1, are carried out. To verify the proposed procedure and developed code for CF VIV analysis of the slender body, several riser model tests with uniform and shear current are simulated and compared with the results of model tests.

2. Synchronization Model for Accurate CF VIV Loads

2.1. Synchronization Model

The flow around the circular cylinder generates forces in the flow direction and also perpendicular to flow. As shown in Figure 1, the force perpendicular to flow is defined as lift, and the force in the flow direction is defined as drag. The lift force of the fixed circular cylinder is caused by the pressure fluctuation on the surface.

Since the lift force on the fixed circular cylinder has oscillated in the frequency of the vortex shedding, it is easy to predict the force. However, since the frequency of the vortex shedding changes in synchronization with the oscillating circular cylinder, it is difficult to predict lift force on an oscillating circular cylinder. Vortex-induced vibration has the characteristic that the vortex shedding is synchronized with the oscillating object. The characteristics of VIV are similar to the synchronization phenomenon of fireflies, in that the firefly’s light is synchronized with the surrounding light cycle. A theory simulating the synchronization of fireflies was developed by Izhikevich and Kuramoto [31] using a weakly coupled oscillator. Based on the theory of Izhikevich and Kuramoto [31], Thorsen et al. [25] derived the synchronization model, which evaluates the load due to vortex-induced vibration. Through comparison with free and forced vibration tests, Thorsen et al. [25] confirm that the derived synchronization model mimics the VIV phenomenon well. The synchronization model [25] (hereafter synchronization I model) provides the force only perpendicular to flow as Equation (1).

$$F_{VIV1}=\frac{1}{2}\rho D\left|\overrightarrow{V}\right|U{C}_L\left(A/D\right)cos\left({\varphi }_L\right)-\frac{1}{2}\rho \left(D{C}_{1}U\dot{y}+A{C}_{2}U\left|\dot{y}\right|\dot{y}\right)-\frac{\rho \pi D^2}{4}{C}_a\ddot{y}$$

(1)

where $\rho$ is the density of the fluid, $D$ is the outer diameter of the circular cylinder, $U$ is the velocity of the incoming flow, and $\overrightarrow{V}$ is the relative velocity of the fluid. $A$, $\dot{y,}$ and $\ddot{y}$ are the amplitude, velocity, and acceleration of the cylinder in a perpendicular direction to the incoming flow. Equation (1) consists of three components. The first term of Equation (1), lift force, consists of a phase function $cos\left({\varphi }_L\right)$ and a lift coefficient $C_L\left(\frac{A}{D}\right)$ determined by the normalized response amplitude $\frac{A}{D}$ perpendicular to the incoming flow. Based on the model test of Gopalkrishnan [32] and analysis data of Larsen et al. [12], Thorsen et al. [25] present $C_L\left(\frac{A}{D}\right)$ as shown in Figure 2. Since the lift constant $C_L\left(\frac{A}{D}\right)$ is a function of the normalized response amplitude $\frac{A}{D}$ perpendicular to the incoming flow, the amplitude $A$ should be calculated and entered in real-time.

Since the circular cylinder oscillates, the mean position of the circular cylinder changes depending on time. For this reason, it is difficult to determine the amplitude $A$ for each time step. To overcome this difficulty, Thorsen et al. [25] calculate the amplitude by integrating the velocity $\dot{y}$ of the cylinder using Equation (2), assuming a narrow-band process with the amplitude $A$.

$$A=\frac{1}{2}{{\displaystyle \int }}_{t_a}^{t_b}\left|\dot{y}\right|dt$$

(2)

where $t_a$ and $t_b$ mean the latest zero-crossing time of velocity $\dot{y}$. The ${\varphi }_L$ of the phase function $cos\left({\varphi }_L\right)$ is calculated by the lift phase model of Thorsen et al. [25] as Equation (3).

$$\frac{d{\varphi }_L}{dt}=2\pi f_s+H\left({\varphi }_{\dot{y}}-{\varphi }_L\right)$$

(3)

where $f_s$ is the frequency of vortex shedding which is determined by the Strouhal number $S_t$. $H\left({\varphi }_{\dot{y}}-{\varphi }_L\right)$ is an arbitrary function that provides the synchronization relationship between the lift force and velocity $\dot{y}$ of the circular cylinder and is presented by Thorsen et al. [25] as shown in Figure 3. As shown in Figure 3, the lift phase model determines the time derivative of the phase of the lift force by using the difference between the phase of velocity $\dot{y}$ of the cylinder and the phase of lift force, so the phase of the lift force can be synchronized according to the change of the velocity $\dot{y}$ of the cylinder.

The role of $H\left({\varphi }_{\dot{y}}-{\varphi }_L\right)$ is to accelerate the phase of the lift force when the velocity $\dot{y}$ of the cylinder is faster than the phase of the lift force, and slow down the phase of the lift force when the velocity $\dot{y}$ of the cylinder is lower than the phase of the lift force. I this model, the frequency of vortex shedding synchronizes to oscillating objects is different from the classical lift force model which has a fixed frequency of vortex shedding. Therefore, the synchronization model can more realistically mimic the phenomenon of vortex-induced vibration. In the second term of Equation (1), the damping force is defined as a nonlinear function of the velocity. Using the damping coefficients of Venugopal [33] and the least-squares method, the coefficient of the nonlinear damping force can be determined. Thorsen et al. [25] presented the coefficients of the nonlinear damping force model as $C_1=0.485$ and $C_2=0.936$. The last term of Equation (1) is the force due to the added mass, and $C_a=1.0$ derived from the still water potential theory is used. Since the synchronization model only provides the force perpendicular to incoming flow, a different type of model to calculate forces in other directions of the synchronization model is needed. When calculating the hydrodynamic forces of a slender body, the Morison equation is commonly used. For combining the synchronization model with the Morison equation, Thorsen et al. [25] propose the modified synchronization model (hereafter Synchronization II model) as Equation (4).

$$F_{VIV2}=\frac{1}{2}\rho D\left|\overrightarrow{V}\right|U{C}_{L0}cos\left({\varphi }_L\right)-\frac{1}{2}\rho D{C}_D\dot{y}\left|\overrightarrow{V}\right|-\frac{\rho \pi D^2}{4}{C}_a\ddot{y}$$

(4)

where $C_{L0}$ and $C_D$ are the constant lift coefficient and the Morison drag coefficient, respectively.

2.2. Comparison of Synchronization Models

Synchronization models [25,26] described in Section 2.1 are implemented in a numerical code, and simulations are carried out with the free oscillations of the elastically supported circular cylinder and the forced vibration of the circular cylinder using the developed code. Comparative studies between the results of the simulation and model tests are used to verify the developed numerical code and to observe the characteristics of the two types of synchronization models.

2.2.1. Free Oscillations of Elastically Supported Circular Cylinder

In this study, the free oscillations of the elastically supported circular cylinder are simulated as shown in Figure 4.

Table 1. Simulation parameters, free oscillations [34].

	Low Mass Ratio	Heavy Mass Ratio
$m^{\ast }$ [-]	1.19	8.63
$\zeta$ [-]	0.00502	0.00151
D [m]	0.0381
$f_0$ [Hz]	0.33
$S_t$ [-]	0.2
$\rho$ [kg/m3]	1000
$V_r$ [-]	1–20

shows the simulation conditions [34]. In this simulation, the Newmark β method is used for the time integration of the equation of motions. Simulation conditions for mass ratio $m^{\ast }=m/\left(\rho \pi D^2/4\right)$ of 1.19 and 8.63 are selected, referring to model tests of Govardhan and Williamson [34]. The damping ratio $\zeta$ and diameter of the cylinder are set at the same value as in the model tests [34]. The natural frequency in still water $f_0$ and the fluid density $\rho$ are selected in Table 1. The simulations of the free-oscillation of the elastically supported circular cylinder are carried out according to the change of the reduced velocity $V_r\left(U/\left(f_{0}D\right)\right)$. The amplitude and frequency of the simulations are derived from analyzing the time series except for the initial transient. Each simulation is carried out for two types of synchronization models. In the simulation of the synchronization II model, the Morison drag coefficient $C_D$ is fixed at 1.2 and the lift coefficient $C_{L0}$ uses 0.8, 1.0, and 1.2, respectively.

Figure 5 shows the nondimensionalized oscillating amplitude and oscillating frequency of free oscillation for an elastically supported circular cylinder according to the reduced velocity $V_r$. In simulation results, compared to heavy mass ratio conditions, low mass ratio conditions maintain high amplitude in a longer range of the reduced velocity $V_r$. This phenomenon is also observed in the results of the model test of Govardhan and Williamson [34]. Figure 5b,d show that the lock-in phenomenon is well mimicked where the range of the reduced velocity $V_r$ is 4–6. In the region of the reduced velocity, $V_r$ is high, and a slight difference is observed between the results of simulations and model tests. However, two results show a good agreement in the overall range. These trends are similarly observed in the calculations of Thorsen et al. [25]. In the simulations of the synchronization II model [26], the oscillating amplitude increases as $C_{L0}$ increases. However, there is little effect on oscillating frequency as $C_{L0}$ increases. Through the simulation results, the possibility of predicting the VIV loads of arbitrary sections is expected by tuning $C_{L0}$ and $f_s$ of the synchronization II model [26]. For circular cylinders, simulations are most similar to the test results in the case of $C_{L0}=0.8$.

2.2.2. Forced Oscillations of Circular Cylinder

Following the simulation of free oscillation, simulations for the forced oscillations of a circular cylinder are carried out as shown in Figure 6.

Forced oscillations are also simulated with two types of synchronization models, the synchronization I model [25] and the synchronization II model [26], respectively. For the synchronization II model, $C_{L0}$ and $C_D$ are selected as 0.8 and 1.2 in the simulations of forced oscillations, based on simulation results of free oscillations.

Table 2. Simulation parameters, forced oscillations [32].

Item	Quantity
U [m/s]	0.4
D [m]	0.0254
St	0.2
$\rho$ [kg/m3]	1000
Nondimensional frequency range [-]	0.05–0.35
Nondimensional amplitude range [-]	0.1–1.2

shows simulation conditions for forced oscillations of a circular cylinder. The calculated loads are nondimensionalized using Equation (5).

$$C \left(t\right)=\frac{F \left(t\right)}{0.5\rho D{U}^2}$$

(5)

The nondimensionalized loads are separated into an in-phase and an out-phase component using Equation (6). The separated components are shown as contours in Figure 7.

$$\begin{array}{c}{C}_{y,v}=\underset{T\to \infty }{lim}\frac{2}{T}{{\displaystyle \int }}_0^{T}C\left(t\right)cos\left(\omega t\right)dt\\ C_{y,a}=-\underset{T\to \infty }{lim}\frac{2}{T}{{\displaystyle \int }}_0^{T}C\left(t\right)sin\left(\omega t\right)dt\end{array}$$

(6)

Since the curve, where the value is 0 in the contour of the out-phase component, is a power input or output boundary, this curve is important to analyze the free oscillating cylinder. It is observed that the contours of the out-phase component calculated by the two synchronization models are similar to the model test of Gopalkrishnan [32]. In particular, curves representing power input or output boundary are also well-mimicked. In the simulations of the synchronization II model, characteristics of the power input and output are similarly calculated by simply adjusting $C_{L0}$. The contour of the in-phase component has the characteristic of showing a low slope on the left, while a high slope on the right is bounded by the curve where the value is 0. Compared with the model test of Gopalkrishnan [32] and the two synchronization models, the overall tendencies are similar, although there are slight differences. Through the results of the simulation, it is confirmed that synchronization models qualitatively reflect the physical phenomenon of vortex-induced vibration.

The comparative study shows that, unlike the synchronization I model, it is difficult to select an appropriate lift and drag coefficient to estimate the VIV load, but it is easy to adjust the VIV load of an arbitrary cross-section.

3. Lumped-Mass Line Model for Static and Dynamic Analysis of Slender Structures

3.1. Lumped-Mass Line Model

In this study, the lumped-mass line model is developed for static and dynamic analysis of slender structures [29,35,36]. In the lumped-mass line model, it is discretized into $N$ lines and $N+1$ nodes as shown in Figure 8. Figure 8 shows the right-handed Cartesian coordinate system while the vertical axis z is directed upward for the undisturbed water surface. Each node $X_i$ is defined as a three-dimensional vector such as ${\left[x_i, y_i, z_i\right]}^T$ at $i$ point. The properties of each line, which are outer diameter, inner diameter, dry weight, bending rigidity, axial stiffness, axial damping coefficient, and inner fluid density, are defined at $i+1/2$ point. Through the above-defined line properties, internal forces and external forces are defined at each node of the slender structure as shown in Figure 9.

Internal forces consist of tension, the axial damping force, and shear force due to the bending moment. External forces consist of buoyancy force, weight due to gravity, ground contact force, and hydrodynamic force described Morison equation. Based on the above-defined mass, additional mass, and internal and external forces acting on the slender body, the three-dimensional equation of motion of the slender body is derived in Equation (7) [36].

$$\left[M_i+m_i\right]{\ddot{X}}_i=t_{i+\frac{1}{2}}-t_{i-\frac{1}{2}}+c_{i+\frac{1}{2}}-c_{i-\frac{1}{2}}+s{f}_{i+1/2}-s{f}_{i-1/2}+w_i+b_i+d_{n,i}+d_{t,i}+a{f}_{n,i}$$

(7)

where $M_i$ is the mass, $m_i$ is the added mass, $t_{i+1/2}$ is the tension, $c_{i+1/2}$ is the axial damping force, $s{f}_{i+1/2}$ is the shear force due to bending stiffness, $w_i$ is the wet weight, $b_i$ is the ground contact force, $d_{n,i}$ is the normal Morison drag, $d_{t,i}$ is the tangential drag, and $a{f}_{n,i}$ represents the force due to the acceleration of the fluid particle. Each component is defined in each unit segment. Detailed definitions of the internal and external forces in Equation (7) can be found in Oh et al. [36].

In this study, the dynamic relaxation method applies to calculate the static configuration of a slender body with the current. The dynamic relaxation method is a numerical technique used for form-finding of cables and thin-film structures with large geometrical nonlinearity and has been applied and developed by several researchers [29,37,38,39] after Day [40] proposed it. This method is numerically robust to the calculation of slender bodies with large geometric nonlinearities such as SCR and SLWR, as well as tensioned risers. Since there is no need to calculate the stiffness matrix directly, the numerical algorithm is efficient and simple. The dynamic relaxation method satisfies the static equilibrium by using the virtual dynamic equilibrium equation as Equation (8).

$$R_i^t=V{M}_i·a_i^t+V{C}_i·v_i^t$$

(8)

where $V{M}_i$ is virtual mass, $V{C}_i$ is virtual damping coefficient, $a_i^t$ is virtual acceleration, and $v_i^t$ is virtual velocity. $R_i^t$ is the residual force, which means the static internal and external forces acting on each node of the slender body. Equation (9) is defined using the steady component in Equation (7).

$$R_i^t=t_{i+1/2}-t_{i-\frac{1}{2}}+c_{i+\frac{1}{2}}-c_{i-\frac{1}{2}}+s{f}_{i+\frac{1}{2}}-s{f}_{i-\frac{1}{2}}+w_i+b_i+d_{n,i}+d_{t,i}$$

(9)

Assuming that the virtual velocity $v_i^t$ changes linearly with the increment of time $\Delta t$, and defining the virtual acceleration $a_i^t$ as a linear interpolation for $\Delta t$, Equation (8) can be summarized as Equation (10) for incremental velocity $v_i^{t+\Delta t/2}$ [29,38,39].

$$v_i^{t+\Delta t/2}=v_i^{t-\Delta t/2}\left(\frac{V{M}_i/\Delta t-V{C}_i/2}{V{M}_i/\Delta t+V{C}_i/2}\right)+\left(\frac{R_i^t}{V{M}_i/\Delta t+V{C}_i/2}\right)$$

(10)

Using the incremental velocity $v_i^{t+\Delta t/2}$ defined in Equation (10), the incremental displacement is defined as Equation (11).

$$X_i^{t+\Delta t/2}=X_i^t+\frac{\Delta t}{2}{v}_i^{t+\Delta t/2}$$

(11)

Using Equation (11), the incremental displacement is calculated until the incremental velocity becomes less than the reference velocity at the initial position of the slender body. Since the incremental velocity is due to the residual force, when the incremental velocity becomes zero, the residual force has achieved static equilibrium. By appropriately selecting the virtual mass and the virtual damping coefficient, the stability and convergence of the dynamic relaxation method can be improved. In this study, the same virtual mass and virtual damping coefficient are applied as in Oh et al. [29].

The dynamic analysis is performed with static equilibrium as an initial condition using the three-dimensional equation of motion of a slender body as Equation (7). The fourth-order Runge–Kutta method, which is known to have excellent stability and accuracy, is used for time integration.

3.2. Validation of Lumped-Mass Line Model for Statistic and Dynamic Analysis

To validate the lumped-mass line model and dynamic relaxation method, a static and dynamic analysis of SLWR (steel lazy wave riser) is carried out in this study. The results of developed code are compared with those of OrcaFlex 10.1, a commercial software. The main dimensions of SLWR are referenced in the study of Ruan et al. [41] as shown in

Table 3. Parameter of steel lazy wave riser [41].

Item	Hang-Off Riser	Buoyancy Riser	Touch Down Riser
Segment length [m]	1600	370	1200
Outer diameter [m]	0.457	1.137	0.457
Inner diameter [m]	0.406
Dry weight [kg/m]	270	697	270
Bending rigidity [N·m2]	1.66 × 108
Axial stiffness [N]	7.098 × 109
Inner flow density [kg·m3]	881
Seabed stiffness [N/m/m2]	6000

. Unlike the study of Oh et al. [29,36], static and dynamic simulations of SLWR are carried out by considering the current. Since SLWR has strong geometrical nonlinearity, it is a good example to confirm the robustness of the analysis code for the slender body.

The hydrodynamic force coefficients for static and dynamic analysis with the current are defined as shown in

Table 4. Hydrodynamic coefficients for SLWR analysis.

Item	Quantity
Drag coefficient	1.2
Friction coefficient	0.08
Normal added mass coefficient	1.0
Tangential added mass coefficient	0.08

. The hang-off point and anchor point of SLWR are (0, 0, −10) and (−2340, 0, 1255), respectively. The profile of the current is the shear current which has a maximum velocity at the free surface and decreases linearly to the bottom as shown in Figure 10. A total of 317 elements are used for static and dynamic analysis of SLWR. Static analyses are carried out according to the variation of the maximum velocity of the current to −2.5, 0, and 2.5 m/s. Figure 11 shows the comparison between two results for the static configuration, effective tension, and bending moment. The results of the developed code are in good agreement with those of OrcaFlex 10.1, a commercial code.

In the case of the dynamic analysis, the top of SLWR under shear current is assumed to be excited in motion with an amplitude and period of the oscillating motion. The simulation conditions are 1.54 m of amplitude and 8 s period of the oscillating motion, and 2.5 m/s of the maximum current velocity. The time series of tension and bending moment at four points shown in Figure 12 are observed and compared with the results of OrcaFlex 10.1, commercial code. The tension or bending moment at each point of SLWR oscillates to the excitation period. In particular, the time histories of bending moments show strong nonlinear characteristics. Results of dynamic analysis, calculated by developed code, are also in good agreement with the results of OrcaFlex 10.1. Through comparison studies between the developed code and commercial code, the developed code is validated for static and dynamic analysis of slender bodies with the current.

4. CF VIV Analysis of Slender Body with the Current

4.1. Proposed CF VIV Analysis Model

In this study, CF VIV analysis of the slender body with the current is carried out using the synchronization I model and lumped-mass line model, which are validated in the previous section. To simulate the VIV of the slender body due to the current, the calculation procedure is proposed as shown in Figure 13.

The calculation procedure is summarized in three stages. First, the static configuration of the slender body is calculated due to internal forces and external forces. Internal forces consist of tension and shear forces due to bending moments, and external forces consist of self-weight, buoyancy, ground-contact force, and drag force due to the current. In the second step, based on the static configuration of the slender body caused by the current, the effective incoming flow velocity ${\overrightarrow{U}}_{NV,i}$ at each cross-section is calculated at each node as shown in Figure 14. The unit normal vector ${\overrightarrow{n}}_{cross, i}$ perpendicular to the effective incoming flow direction is also calculated using tangential vector ${\overrightarrow{TV}}_i$ and ${\overrightarrow{U}}_{NV,i}$ as follows Equation (12).

$${\overrightarrow{n}}_{cross, i}=\frac{{\overrightarrow{TV}}_i\times {\overrightarrow{U}}_{NV,i}}{\left|{\overrightarrow{U}}_{NV,i}\right|}$$

(12)

Finally, dynamic analysis of the slender body is carried out using the VIV load, considered by ${\overrightarrow{U}}_{NV,i}$ and ${\overrightarrow{n}}_{cross, i}$. In this study, it is assumed that the static configuration of the slender body with current is the mean position for defining ${\overrightarrow{U}}_{NV,i}$. Therefore, the VIV load always acts perpendicular to ${\overrightarrow{U}}_{NV,i}$ of the static configuration. When the amplitude of the vortex-induced vibration is greater than about 1.2 times the diameter, the lift force decreases dramatically and the damping force becomes dominant. Therefore, it is expected that the maximum dynamic amplitude is very small compared with the length of the slender body. The assumption, that the static configuration is the mean position for defining ${\overrightarrow{U}}_{NV,i}$, is reasonable. Equations of motion (13) are derived by adding the VIV term to Equation (7). Unlike the study of Thorsen et al. [26], applying the simplified synchronization II model and Morison drag in the CF direction to which the VIV load is applied was eliminated to apply the more sophisticated synchronization I model. Since the drag force term of VIV load and Morison equation overlap, the drag force term of Morison equation perpendicular to the effective inflow flow is eliminated.

$$\left[M_i+m_i\right]{\ddot{X}}_i=t_{i+\frac{1}{2}}-t_{i-\frac{1}{2}}+c_{i+\frac{1}{2}}-c_{i-\frac{1}{2}}+s{f}_{i+\frac{1}{2}}-s{f}_{i-\frac{1}{2}}+w_i+b_i+d_{n,i}+d_{t,i}+f_{viv,i}·{\overrightarrow{n}}_{cross, i}-d_{n,i}·{\overrightarrow{n}}_{cross, i}$$

(13)

In Equation (13), $f_{viv,i}$ means the VIV load using the synchronization I model and is defined considering the mean length $dl$ of each node as Equation (14). The added mass term is already defined in Equation (13) and is therefore eliminated in Equation (14).

$$f_{VIV,i}=\left[\frac{1}{2}\rho D\left|{\overrightarrow{V}}_i\right|\left|{\overrightarrow{U}}_{NV,i}\right|C_L\left(A/D\right)cos\left({\varphi }_L\right)-\frac{1}{2}\rho \left(D{C}_1\left|{\overrightarrow{U}}_{NV,i}\right|{\dot{y}}_{cross}+A{C}_2\left|{\overrightarrow{U}}_{NV,i}\right|\left|{\dot{y}}_{cross}\right|{\dot{y}}_{cross}\right)\right]dl$$

(14)

In Equation (14), ${\overrightarrow{V}}_i$ means the relative velocity of the slender body. $A$ and ${\dot{y}}_{cross}$ mean the amplitude and velocity perpendicular to the effective incoming flow and are defined using ${\overrightarrow{n}}_{cross, i}$. In this study, since it is limited only to the CF VIV analysis caused by the current, $a{f}_{n,i}$ due to the acceleration of fluid particles is also eliminated in Equation (13).

4.2. Validation of CF VIV Analysis Model

In this study, the CF VIV analysis code for the slender body applying the proposed procedure and synchronization I model is developed. To validate the developed code, tension riser under uniform current and shear current are simulated, where relatively numerous results have been published [5,7].

First, the simulations of tensioned riser under uniform current [7] are carried out. The riser model is divided into 50 elements in these simulations. These tests implement a uniform current through towing equipment moving at a constant speed, as shown in Figure 15. Simulation conditions for the tensioned riser under a uniform current are shown in

Table 5. Parameters of tensioned riser model with uniform current.

	Exp. of Song et al. [7]
Length [m]	7.9
Outer diameter [m]	0.031
Inner diameter [m]	0.027
Bending stiffness [N m2]	1476.763
Axial stiffness [N]	1.45 × 107
Mass per unit length [kg/m]	1.768
Pretension [N]	2943
Current velocity [m/s]	0.4
$S_t$	0.2

. Morison drag and added mass coefficients are assumed as shown in Table 4.

Figure 16a,b show the results of the VIV simulation under the uniform current for the model test case by Song et al. [7]. In the upper panel of Figure 16a, the calculated static configuration in the inline direction is also in exact agreement with that of the commercial code. In the lower panel of Figure 16a, the magnitude of the CF RMS amplitude and the main excitation mode of the RMS of CF displacement in this study are similar to the measurements by Song et al. [7]. The space-time plot of the CF displacement also shows the main excitation mode over time as shown in Figure 16b. The 1st mode in CF VIV is observed in Song et al.’s case [7] with the fairly long model.

In order to validate the VIV simulation of the tension riser under the shear current, simulations for the high-mode VIV model test [5] performed in the Norwegian deep-water program (NDP) are carried out as shown in Figure 17.

Table 6. Parameters of tensioned riser model with shear current.

	Exp. of Kristiansen and Lie [5]
Length [m]	38
Outer diameter [m]	0.027
Bending stiffness [N m2]	599
Axial stiffness [N]	8.118 × 106
Mass per unit length [kg/m]	0.933
Pretension [N]	4300–4600
St [-]	0.19
Maximum current velocity [m/s]	0.6, 1.3, 2.0

shows the parameters of the tensioned riser model of the high-mode VIV test performed in NDP. The riser model is also divided into 50 elements in these simulations. In the riser simulations under shear current conditions, the Morison drag and the added mass coefficient are also assumed as shown in Table 4. The VIV simulations of the riser are carried out under shear current conditions of 0.6, 1.3, and 2.0 m/s, respectively, with maximum velocity $U_{max}$. The simulation results are compared with those of Kristiansen and Lie [5].

For the shear current condition with 0.6m/s of maximum velocity $U_{max}$, the static configuration of the inline displacement and the RMS of CF displacement are shown in Figure 18a.

In Figure 18a, the inline static configuration of the present simulation agrees exactly with the result of OrcaFlex 10.1. The magnitude of the CF RMS amplitude and the main excitation mode of the RMS of CF displacement of the present simulation also agree excellently with the results of the model tested by Kristiansen and Lie [5] and the calculation by Thorsen et al. [25], respectively, in Figure 18a. Figure 18b shows the space–time plot of the CF displacement. Unlike the uniform current condition, the wave of CF displacement under the shear current condition propagates from the section where the velocity of the incoming flow is fast to the section where the velocity of the incoming flow is slow.

Figure 19a and Figure 20a show the inline static configuration and the RMS of CF displacement in the conditions of the maximum velocity $U_{max}$ of shear current 1.3 m/s and 2.0 m/s, respectively. The magnitude of the CF RMS amplitude and the main excitation mode of the RMS of CF displacement in this study also agree well with the results of the model test by Kristiansen and Lie [5] and the calculation by Thorsen et al. [25]. Figure 19b and Figure 20b also show the space-time plot of the CF displacement. The wave of CF displacement also propagates from the section of fast incoming flow to the section of slow incoming flow, such as the conditions of the maximum velocity $U_{max}$ of shear current 0.6 m/s. However, in the case of $U_{max}=1.3, 2.0 \mathrm{m}/\mathrm{s}$, it is observed that standing waves are generated due to the reflection of the fixed boundary condition at the section of the slow incoming flow. These standing waves have also been reported by the calculation of Thorsen et al. [25]. These propagated waves due to VIV mean the propagation of energy from fast velocity region to a slow velocity region and this phenomenon is not observed from the uniform current condition.

In order to analyze the VIV characteristics of a tensioned riser under shear current, VIV responses at three points (near the top, middle, and bottom) of the riser under the shear current of $U_{max}=1.3 \mathrm{m}/\mathrm{s}$ are presented in the time and frequency domains in Figure 21. Figure 21a shows the VIV responses at the three points in the time history. Unlike regular wave response in the CF on the uniform current condition, an irregular response with multi-frequencies is observed at all three points. Multi-frequencies components of VIV responses from other points due to the shear current are presumed to propagate and be superposed along the riser. The phenomenon is further analyzed with the frequency results in Figure 21b. The peak frequency of 7.5 (Hz), resonance frequency with the 9th natural frequency (7.335 Hz) of the riser, is observed at different three points. Considering the natural frequencies of the riser model, VIV components with 7.5 (Hz) of frequency induced by 1.04 m/s of flow speed at −7.6 m of the riser excites the riser, resulting in resonant with the 9th natural frequency of the riser. It is found that a tensioned riser under the shear current vibrates in crossflow with multi-frequencies riding one dominant frequency.

Figure 22 shows the CF VIV snapshots of the riser in the time history. It is confirmed that the tensioned riser is dominantly excited with the 9th natural frequency by vortex shedding.

In this study, CF VIV analysis of the tensioned riser for the uniform current and shear current condition is carried out and compared with the model test. Through the comparative study, it is observed that the magnitude of the CF RMS amplitude and the main excitation mode of the RMS of CF displacement, which is calculated by the developed code in this study, are similar to the results of the model tests. Therefore, it is confirmed that the developed code provides a reasonable prediction of VIV of a slender structure.

5. Conclusions

In this paper, the numerical analysis method to predict CF VIV for slender structures under the current is developed based on the proposed three stage procedure. A slender structure is modeled as the lumped-mass line, and the dynamic relaxation method known as the numerical technique for slender structures with large geometric nonlinearity is applied in the static analysis. The fourth-order Runge–Kutta method is used for the time integration. The vortex-induced force on slender structures is considered with the synchronization model. In order to apply the more sophisticated synchronization I model [25] compared to synchronization II model [26], Morison drag in the CF direction to which the VIV load is applied was eliminated. The following conclusions from the case studies with the developed numerical method are found.

(1)

The developed numerical code with the synchronization models exactly predicts the CF VIV of the elastically supported circular cylinder. It is shown that the amplitudes of the low mass ratio are maintained higher than the amplitudes of the heavy mass ratio in the long range of the reduced velocity. In addition, it is observed that the numerical code predicts the lock-in phenomenon well in the four–six regions of reduced velocity. From the simulations for forced oscillations of the circular cylinder, it is observed that the power input or output boundary is qualitatively well mimicked using the two synchronization models. It is found that the developed numerical code with two synchronization models provides an accurate VIV load. It is difficult to select appropriate lift and drag coefficient with the synchronization II model without the experimental data. However, since the expression of the synchronization II model is simple, modeling the VIV load of the arbitrary section from the data of experiments is promising when $C_{L0}$, $C_D$, and $f_s$ are adjusted.

(2)

The lumped-mass line model for the slender structures is used for the static and dynamic analysis. The results for the static and dynamic analysis of the SLWR with highly geometric nonlinear are in good agreement with commercial code OrcaFlex 10.1. The accuracy and robustness of the developed code are confirmed.

(3)

In the validation of the developed numerical method to predict CF VIV of a tensioned riser under uniform and shear current, it is shown that CF RMS amplitudes and main excitation modes of the developed numerical code show fairly good agreement with the measurements of the model tests. In the shear current condition, it is observed that the wave of CF displacement propagates from the section of fast incoming flow to the section of slow incoming flow, unlike the uniform flow condition. It is found that a tensioned riser vibrates with one dominant frequency in resonance with the n^th mode, even though multi-frequencies components of the vortex shedding along the riser due to the shear current occur. By comparison between the results of developed numerical code and model tests, the proposed VIV analysis procedure and developed numerical method are validated.

(4)

It is confirmed that the developed numerical method is able to predict the VIV of a slender structure, and will contribute to estimating the VIV fatigue of a slender structure.

There is a plan to perform CF VIV analysis of SCR or SLWR, which has highly geometric nonlinearity, using the proposed procedure and developed code in the future.