Two improvements on numerical simulation of 2-DOF vortex-induced vibration with low mass ratio

Till now, there have been lots of researches on numerical simulation of vortex-induced vibration. Acceptable results have been obtained for fixed cylinders with low Reynolds number. However, for responses of 2-DOF vortex-induced vibration with low mass ratio, the accuracy is not satisfactory, especially for the maximum amplitudes. In Jauvtis and Williamson’s work, the maximum amplitude of the cylinder with low mass ratio m*=2.6 can reach as large as 1.5D to be called as the “super-upper branch”, but from current literatures, few simulation results can achieve such value, even fail to capture the upper branch. Besides, it is found that the amplitude decays too fast in the lower branch with the RANS-based turbulence model. The reason is likely to be the defects of the turbulence model itself in the prediction of unsteady separated flows as well as the unreasonable setting of the numerical simulation parameters. Aiming at above issues, a modified turbulence model is proposed in this paper, and the effect of the acceleration of flow field on the response of vortex-induced vibration is studied based on OpenFOAM. By analyzing the responses of amplitude, phase and trajectory, frequency and vortex mode, it is proved that the vortex-induced vibration can be predicted accurately with the modified turbulence model under appropriate flow field acceleration.


Introduction
With the development of marine industry, the phenomenon of vortex-induced vibration has been paid more and more attention (Ren et al., 2014). Numerical simulation is an important method to study vortex-induced vibration. For designers, a practical method is required to provide accurate estimation of the responses of vortex-induced vibration parameters to estimate the dynamic response, nonlinear impact for fatigue analysis and reliability assessment of marine slender bodies (Zhang et al., 2016).
At present, the application of DNS (Direct Numerical Simulation) and LES (Large Eddy Simulations) as tools for engineering predictions still has difficulties to be applied to high Reynolds number condition due to the high cost inherent in need of using very fine grids and small time steps (Sagaut and Deck, 2009). RANS is still considered to be the practical approach for engineering problem whose basic idea is to divide the volume of the system into the average quantity and the fluctuating quantity by the ensemble average of the variables of the turbulent flow field. Among which, SST turbulence model, combined with the advantage of k-ε and k-ω models (Menter, 1994), is considered to have higher precision and reliability of the calculated results compared with other Reynolds-averaged model. However, its performance seems to be not satisfactory in predicting flows where vortex shedding occurs.
RANS-based turbulence models are known to yield less accurate predictions of flows with strong anisotropic turbulence (Ong et al., 2009). When the Reynolds number is large, the numerical response of the lift coefficient of flow around a circular cylinder is often smaller than the experimental values. Stringer et al. (2014) reported poor agreement with experimental data of flow around a circular cylinder when Reynolds number is relatively large, with the SST turbulence model. Four standard RANS-based turbulence models were evaluated by Ünal et al. (2010) on the simulation of the near-wake flow of a circular cylinder, and the results showed that the SST turbulence model is obviously more encouraging to produce more accurate results, but still not satisfactory enough. Meanwhile, in simulations of vortex-induced vibration, the transverse amplitudes of the cylinders are often smaller than the experimental values and the amplitude decays too fast in the lower branch. In Jauvtis and  experiment, when the mass ra-tio is 2.6, the maximum transverse amplitude of the cylinder can reach 1.5D (D is the diameter of the circular cylinder), which has been seldom captured in numerical simulation. Pan et al. (2007) failed to capture the upper branch in simulation of the response of an elastically mounted rigid cylinder at low mass ratio constrained to oscillate transversely to a free stream with the standard SST turbulence model. Similarly, Sanchis et al. (2008) also reported not to capture the upper branch. Similarly, the maximum transverse amplitudes in Gu et al. (2014) and Lin and Wang's (2013) work are only about 0.8, and the amplitude responses decay obviously too fast when the reduction velocity is larger than 9.
Many researchers have investigated the reason for RANS turbulence model to perform poorly in predicting the unsteady separated flows. Medic (1999) pointed out that the standard RANS-based turbulence closure have shortcomings originally in the use of dissipative discretization schemes, and will finally lead to numerical diffusion errors that are adequately large to suppress the vortex shedding process when applied in predicting the unsteady separated flows. Besides, Younis and Przulj (2006) summarized the previous theoretical and experimental results, and proposed that turbulence model should take into account the effects of the interactions between the large-scale organized periodicity of the mean flow and the random, small scale high-frequency motions that characterize turbulence, which has been reflected in none of the standard form of RANS-based turbulence models.
Moreover, according to Jauvtis and Williamson (2004), the upper branch of vortex-induced vibration can only be captured in acceleration process. However, till now, little research has been done in investigating the relationship between acceleration magnitude and transverse amplitude response. In order to search the proper acceleration value, a self-compiled script is adopted in this paper to increase the flow velocity continuously in the process of simulation, the acceleration process of two degrees of freedom (2-DOF) vortex-induced vibration of a cylinder with the reduced velocity U * ( , where U is the inflow velocity, D is the diameter of the cylinder, and f n is the natural frequency) from 0 to 6 under different acceleration magnitudes is simulated, and then the frequencies and transverse amplitudes are analyzed. Generally, in this paper, a modified SST turbulence model is proposed and applied for the numerical simulation based on OpenFOAM. Firstly, the proper acceleration magnitude of flow field is determined by analyzing the responses of the cylinder under different acceleration values. Then the 2-DOF vortex-induced vibration of a cylinder with a mass ratio of 2.6 under different reduced velocities is simulated with a proper acceleration value. The numerical results are compared in detail with the experimental data in-cluding lift and drag coefficients, amplitude, phase angle, frequency, trajectory and vortex model. It is proved that the vortex-induced vibration can be predicted accurately with the modified turbulence model under appropriate acceleration of flow field.

Modified turbulence model
The time-averaged turbulent flow governing equations can be expressed as: where u and p are the time-averaged velocity and pressure, respectively; μ is the molecular viscosity; and ρ is the fluid density.
Turbulence models are used to enclose the governing equations. The k-ε model is one of the most common turbulence models. It is a two equation model, which means that it includes two extra transport equations to represent the turbulent properties of the flow. This allows a two-equation model to account for historic effects like convection and diffusion of turbulent energy. The dimensionless form of standard k-ε model can be expressed as follows: Transport equations for turbulent kinetic energy k: for dissipation ε: Turbulent viscosity: Production of k: where the small-scale fluctuations of the velocity related to the turbulence are reduced as , which is referred to as the Reynolds stresses.
, , , and are the model constants, whose values are shown in Table 1.
As is known, the standard k-ε model does not perform well in the case of reverse pressure gradient (Wilcox, 1998). Younis and Przulj (2006) drew a conclusion that the estab- lishment of a vortex shedding field in a turbulent flow leads to direct energy input from the periodic mean-flow oscillations into the random turbulence motions and this direct energy supply occurs at a frequency that corresponds exactly to the vortex-shedding frequency. The derivation procedure is as follows.
When the vortex shedding occurs, the shape of the distorted spectrum can be expressed as (Reynolds, 1974): where is a constant and is the wavenumber vector. s is a matching index whose precise value is immaterial to the present discussion. must vanish in the steady limit.

E(κ, t)
According to the relationship between kinetic energy k and via (Cant, 2001): (9) An equation for the change of the dissipation rate ε with time is obtained (Younis and Przulj, 2006): where is the coefficient in the standard k-ε turbulence model, and .
According to the above formula, the coefficient is modified as: where Q is the mean-flow kinetic energy per unit mass and C t is obtained by numerical optimization and is assigned to 0.38, and is an additional term, which is aimed to consider the effect of the vortex shedding.
Considering that SST turbulence model has higher accuracy, we transform the additional term R into the term ω based on the relationship , and add it to the specific dissipation rate equation in SST turbulence closures, the expression of R can be obtained as: β ′ where is an empirical coefficient obtained by numerical optimization and is assigned to 0.54 after a large number of test. Add R to the specific dissipation rate equation of standard SST model, the modified SST specific dissipation rate equation is expressed as: Other equations are the same with the ones in standard SST turbulence closures which can be seen in Menter (1994).
The cylinder can be reduced to a mass spring damping system, and the governing equations of its motion can be expressed as: My where M represents the mass of the cylinder, K and C respectively represents the spring stiffness and damping coefficient, and x(t) and y(t) respectively mean the in-line and cross-flow displacement.

Computational domain and grid
The grids used in the numerical simulation are shown in Fig. 1. The details of the size of the computational domains are as follows. Definition of coordinate system: the origin is the center point of the cylinder, the directions of the x axis, y axis and z axis are the flow direction, the vertical flow direction, and the thickness direction of the cylinder, respectively.
Range of computational domain: according to the experiment of Jauvtis and Williamson (2004) experiment, set D=0.0381 m, the range of the computational domain is -12D<x<25D, -12D<y<12D and a layer of mesh thickness in the z axial direction.
Grid division: The O type grids are generated in the range of 3D around the cylinder and hexahedral grids in other areas. The number of grids is adjusted according to different Reynolds numbers. The thickness of the first layer on the surface of the cylinder is also adjusted to insure y+<1. The details of the girds are shown in Table 2. The grid independence is analyzed to meet the accuracy requirement of simulation error.
Time step: The time step is determined according to the Courant number (Venkatachari et al., 2008) which can be expressed as: is the time step, and is the grid size. In the numerical simulation, the time step is adjusted according to u to ensure the maximum Courant number is about 0.2. Boundary conditions: The no slip condition is employed on the cylinder surface. On the right exit boundary, the hydrostatic pressure is set to 0. The inflow velocity is controlled by self-script. The rest are symmetric boundaries.
In order to be consistent with the experimental situation, the parameters of the cylinder in the numerical simulation are set as follows: the mass m=1.138 kg, mass ratio m * =2.6, the damping , mass damping ratio and natural frequency both for the in-line and crossflow motions.

Results and discussion
3.1 Response of lift coefficient C l of flow around a fixed circular cylinder According to the mechanism of turbulence model modification, the main differences in the responses of vortex-induced vibration with standard and modified SST models lie in the responses of lift coefficient and transverse amplitude. In order to verify the performance of the turbulence model, Gird G1 is applied in this case, and the cylinder is fixed. Standard and modified SST models are used to simulate the flow around a fixed circular cylinder respectively, with the Reynolds number from 5000 to 100000. The lift coefficient is calculated and compared with the experimental data of Norberg (2003), as shown in Fig. 2.
The change trend of the lift coefficient is similar to the experimental results. It increases rapidly when Reynolds number is about 5000 to 20000 and then slows down, remaining nearly constant. However, it is obvious that when the Reynolds number is larger than 5000, with standard SST turbulence model, the lift coefficient is about 20% smaller than the experimental value, indicating that the shedding process is suppressed as mentioned above. After adding the correction term R to consider the effects of the interactions between the large scale and small scale, the modified SST model significantly promoted the accuracy of lift coefficient. The result of the modified SST model shows better agreement with the experimental data.

Response with different acceleration magnitudes
The gird G1 is used in this case and the detailed steps are as follows. The first stage is the uniformly accelerative stage. In this stage, the flow velocity increases uniformly from zero to the target value. In order to avoid the influence of dimension, the flow rate is expressed as the reduced velocity U * . The magnitude of the acceleration is controlled by the acceleration time of U * from 0 to the target value.
The second stage is the uniform stage. After the flow velocity has reached the target value, it is maintained at a constant speed, and the elastically mounted cylinder keeps oscillating until its movement reaches a steady state. Then the responses of the transverse amplitudes and frequencies are analyzed and compared with the experimental results.

∆U *
In order to study whether the vibration response can jump to the upper branch, the target reduced velocity is set to 6, corresponding to the stable stage of the upper branch. The acceleration time is expressed as the normalized time (Ut/D), and the range of normalized acceleration time in this paper is 24-480, corresponding to the acceleration magnitude, which is expressed in a dimensionless form as per normalized time from 0.0125 to 0.25. The constant velocity condition (the initial speed is set to 6 directly) is also included. Fig. 3 presents the response of the transverse vibration frequency ratio response with different acceleration magnitudes, where f y is the transverse oscillating frequency during the induced vibration. When the normalized  (Alam and Sakamoto, 2005), indicating that no lock-in phenomenon occurs in the response. When the normalized acceleration time is smaller than 360, the vibration frequency is locked in the vicinity of the natural frequency (the value of the natural frequency is expressed by the red line), which is one of the characteristics of the transition to the upper branch. The final maximum dimensional transverse amplitude responses with different acceleration magnitudes are shown in Fig. 4. The maximum amplitude value when U * =6 in the experiment of Jauvtis and Williamson is also plotted in the figure. The amplitude is expressed in dimensionless form as A ymax /D, meaning the ratio of the maximum amplitude to the diameter D.
It can be found that when the normalized acceleration time is smaller than 240, corresponding to the acceleration value larger than 0.025 per normalized time, the final maximum transverse response is almost the same with the one under constant velocity condition, which is only about 0.8D, obviously smaller than the experimental value. When the normalized acceleration time is larger than 240, the maximum amplitude increases the increase of the acceleration time until the normalized acceleration time is about 360, corresponding to the acceleration value of 0.017 per normalized time. Then the amplitude reaches the maximum value, equal to about 1D, which is very close to the experimental value and after that it remains constant.
As a kind of nonlinear vibration, the vortex-induced vibration response is not only related to the boundary conditions in the steady state, but also to the process of reaching the steady state. With different flow field accelerations, different vortex-induced vibration responses are presented. From the above analysis, we may draw the conclusion that in order to capture the upper branch in numerical simulation, the flow field acceleration must be smaller than 0.017 per normalized time.

Response of 2-DOF vortex-induced vibration
According to the analysis above, the acceleration is set to 0.015 per normalized time in the simulation of 2-DOF vortex-induced vibration of the cylinder. The acceleration steps are similar to those introduced in Section 3.2.
Only the increasing condition is considered since we mainly focus on the maximum transverse amplitude of the upper branch, which the occurs only in the acceleration process. The responses of amplitudes versus the reduced velocity U * from 2 to 14 under constant and accelerating conditions as well as the experimental data of Jauvtis and Williamson (2004) are shown in Fig. 5.
It can be seen from Fig. 5 that the amplitude under the constant condition can only reach about 0.8D, combined with the analysis result in Section 3.2, lock-in phenomenon is also failed to be captured under this condition, which causes the small amplitude.
On the other hand, an obvious upper branch is captured under the increasing condition with the modified SST model with the flow field acceleration of 0.015 per dimensionless time. In many existing numerical results such as Pan et al.'s (2007) and Sanchis et al.'s (2008), the upper branches are failed to be observed and the maximum amplitudes are obviously smaller than the experimental value. Moreover, the amplitudes are always attenuated too fast when U * is larger than 9. While in this paper, the maximum amplitude is about 1.4D, much closer to the experimental value and the amplitude response curve in this paper agrees well with the experimental data, except for a slight difference in the reduced velocity in which the upper branch shifts to the lower branch, which may be caused by the three-dimensional effect and other interference factors in the actual experiment.
The responses of the in-line amplitude versus the reduced velocity with standard SST and modified SST turbulence models under the increasing condition are shown in Fig. 6. Similar to the amplitude responses in the cross flow direction, three branches are captured under the increasing condition, while in the upper branch, the maximum vertical amplitude with the standard SST turbulence model is only  about 0.22D, siginificantly smaller than the experimental value, which is about 0.32D, and the one with the modified SST turbulence model reaches about 0.29D, much closer to the experimental value. In the initial and lower branches, the responses of the vertical amplitude with the two turbulence models are similar and all agree well with the experimental data.
Seen from above, the amplitude responses of the vortexinduced vibration with the standard SST model in the existing literature often decay too fast when the reduced velocity is large, and always fail to capture the upper branch. While the results with the modified SST model and appropriate flow field acceleration in this paper are in good accordance with the experimental results.
The frequency analyses accompany with three groups of the trajectories and vortex modes of the cylinder under the critical reduced velocities which equal to 3, 6.5 and 10, corresponding respectively to the stable stage in the initial branch, the upper branch, and the lower branch are presented in Fig. 7, as well as the trajectory of Jauvtis and Williamson (2004) and DPIV results. Notice that the results of Jauvtis and Williamson are shown above and the results in this paper are shown below.
Three distinct branches are captured in the frequency response obviously, corresponding to the initial branch, upper branch, and lower branch. As in the initial branch, the transverse vibration frequency ratio increases gradually from about 0.5 to 0.6 with the reduced velocity U * . When U * reaches about 4, a jump phenomenon appears in the frequency response, corresponding to the beginning of the transition period. During this period, the frequency ratio increases rapidly, and is close to f st with the increase of the velocity, until U * reaches about 5. Then the frequency ratio remains stable at about 0.9, indicating that the transverse frequency is locked in the natural frequency. When U * reaches about 7, the vibration progress shifts to the lower branch, and the frequency ratio suddenly shifts to about 1.3 and remains almost stable. The frequency response is in good agreement with the results of Jauvtis and Williamson (2004) on the whole, except for a slight deviation in the transition period from the upper branch to the lower branch, which may be caused by three dimensional effect and other interfering factors in real experiment.
In the initial branch, the 2S mode wake is observed. In the 2S mode, the vortices are shed from the upper surface of the cylinder, and two vortices are generated per vibration cycle, and fall off from the cylinder alternately. With the increase of U * , in the upper branch, the 2T mode  is captured in the numerical simulation. In this vortex mode, two pairs of vortexes are formed in each period behind the cylinder, each vortex pair consists of three vortices which have been marked in the figure and the rotation direction of the vortices is not exactly the same. When it comes to the lower branch, the 2P mode wake is captured. It is seen from Fig. 7, two vortex pairs are generated per cycle. In each vortex pair, there are two vortices rotating in opposite directions and the vortices expand more laterally than for the 2S mode. The key features of the vortex mode are captured in the numerical simulation perfectly and agree well with the DPIV results of Jauvtis and Williamson (2004).
Then we concentrate on the trajectories. In the initial branch, the trajectory of the cylinder is in traditional "figure 8" shape. With the increase of U * , after the vibration shifts to the super-upper branch, the phase angle between the transverse and longitudinal vibration changes, and the trajectory of the "figure 8" shape gradually inclines to the inflow direction and finally becomes a "crescent shape". Finally, when it shifts to the lower branch, the trajectory turns to be "thin figure 8" shape, indicating that the longitudinal vibration is suppressed more strongly, so that the transverse amplitude is significantly larger than the longitudinal one. Similar shape and change trend of the trajectory were also observed in Kang et al.'s (2016) The responses of the lift coefficient , drag coefficient and displacement y/D at some critical reduced velocities are ploted in Fig. 8, in which U * =3, 6 and 7, respectively correspond to the initial branch, the upper branch, and the lower branch. Notice that stages of the acceleration, and uniform speed are included in the numerical simulation process. For U * =3 and 6, the response shown in the figure is in  uniform stage, while for U * =7, the curve shown in Fig. 8c is the response near U * = 7, from the acceleration to uniform stage.
As can be seen from Fig. 8, when U * increases from 3 to 6, the lift coefficient, drag coefficient and displacement of the cylinder significantly increase with the movement shifting from the initial branch to the super upper branch. The maximum value of increases from about 2 to 2.8; the mean value of increases from about 1.7 to 2.5, and the maximum transverse amplitude increases from about 0.2D to 1D. Combined with the frequency analysis above, at this time, the locking phenomenon occurs, the transverse vibration of the cylinder is locked in the vicinity of the natural frequency, and the resonance is generated, so that the amplitude increases rapidly. It can be seen from Fig. 8c, when the reduced velocity reaches about 7, the maximum value of decreases from about 2.8 to 0.5; the mean value of decreases from about 2.5 to 2, and the maximum transverse amplitude decreases from about 1D to 0.6D, which indicates that the vibration response shifts from the upper branch to the lower branch.
The phase angle between the lift force and displacement during the vortex-induced vibration process is investigated through duration curves of the lift coefficient and displacement under critical reduced velocities U * which are drawn in Fig. 8. When the vibration is in the initial branch (U * =3) and upper branch (U * =6), the phase angle between the lift force and displacement is tend to be 0°, and then the lift force promotes the oscillation of the cylinder, so that the amplitude continues to increase in this stage. As shown in Fig. 8c, when U * reaches about 7, an obvious sudden phase change from 0° to 180° after a very short transition period can be observed, corresponding to the shift from the upper branch to the lower branch. Since then, the lift force and displacement are nearly in the same phase, and the vibration is suppressed, which agrees well with the basic theory, as well as the experimental phenomena of vortex-induced vibration .  Fig. 9. It can be seen that three distinct branches are captured in numerical results and the results of numerical simulation are in agreement with the experimental results overall, except for the hopping point from the upper branch to the lower branch, and the reason has been analyzed in the amplitude response analysis shown in Fig. 5.

Response of 1-DOF vortex-induced vibration
The performance of the modified turbulence model in predicting the amplitude response of 1-DOF vortex-induced vibration is also tested. The same grids are used in this process, and then change the degree of freedom from two-dimensional to one-dimensional and set . With the similar methods, the response of the transverse amplitude versus the reduced velocity U * from 2 to 14 is shown in Fig.  10. Moreover, the experimental results of Khalak and Williamson (1996) as well as the numerical simulation results of Lin and Wang (2013) with the standard SST turbulence model are also plotted in Fig. 9.
It is obvious that in the SST results of Lin and Wang (2013), the upper branch is not captured. When U * > 4, the transverse amplitude is suppressed, and is much smaller than the experimental value. While the modified SST result captures the upper branch obviously and agrees well with the experimental data.

Conclusions
Two improvements for the vortex-induced vibration nu-  merical simulation are carried out in this paper, including the modification of turbulence model and investigation of the effect of flow field acceleration on the response of the vortex-induced vibration based on OpenFOAM. The specific work and conclusions of this paper can be summarized as follows.
Firstly, in view of the defects of standard RANS-based turbulence model, which will generate and suppress the separation of the vortex in predicting the unsteady separation flow, a modified turbulence model is proposed.
Secondly, the relationship between the acceleration of the flow field in numerical simulation with responses of the transverse vibration frequency ratio and the maximum transverse amplitude is investigated. It is found that when the acceleration value of flow field is larger than 0.017 per normalized time, the transverse vibration frequency is close to the Strouhal frequency of the fixed cylinder, and the maximum lateral displacement is obviously smaller than the experimental value. Only when the acceleration value is smaller than 0.017 per normalized time, the upper branch will be captured and then the amplitudes be locked in the vicinity of the natural frequency.
Finally, the numerical responses of 2-DOF vortex-induced vibration of the cylinder with the modified SST turbulence model are compared in detail with the experimental data including the lift and drag coefficients, amplitude, frequency, phase, trajectory and vortex model. The features such as the maximum transverse amplitude, mutations of frequency and phase angle between the lift force and displacement, the characteristics and change trends of the trajectory as well as the 2T vortex model are captured. The performance of the modified SST turbulence model in predicting the response of 1-DOF vortex-induced vibration is also tested.
The modified SST turbulence model performs better to numerically simulate the flow around a fixed circular cylinder, both for 2-DOF vortex-induced vibration and 1-DOF vortex-induced vibration, than the standard one. The acceleration magnitude of the flow field is also an important factor.