Numerical Simulation of Vortex-Induced Transverse Vibration of a Cylinder with Very Low Mass Ratio

Vortex-induced vibration (VIV) of flow past different types of cylinders is important in many engineering fields and has led to a large number of fundamental studies. Computational fluid dynamics (CFD) has become a powerful tool to investigate VIV problems. In this paper, the VIV of an elastically supported rigid-body circular cylinder was numerically studied by solving the Reynolds-averaged Navier–Stokes (RANS) equations with the shear stress transport (SST) k - ω turbulence model. The dynamic mesh approach was used to tackle the domain mesh change due to the cylinder motion. The characteristics of the amplitude response of the cylinder and the vortex shedding frequency at different reduced velocities were compared with the experimental data, and good agreement was found regarding the vibration amplitude and the frequency lock-in phenomenon. Results show the influence of added mass was perceptible in the free vibration and in the initial branch with small cylinder vibration amplitude, while it is fairly weak in the upper and lower branches considering the fact that the main response frequency of the cylinder is close to the natural frequency in a vacuum. Detailed vorticity distributions clearly indicate the 2P vortex modes at the upper and lower branches, and the 2S mode at the initial branch and their generation mechanisms are analyzed. The effects of the cylinder vibration on the lift and drag forces at different branches are discussed, and the frequencies of peak values in their Fourier spectra are examined. The phase difference between the lift force and the cylinder vibration shows an obvious jump from the initial branch to the upper branch with the vortex mode switching from 2S to 2P.


Introduction
Vortex-induced vibration (VIV) of ow past blu bodies is part of a number of disciplines involving uid mechanics, structural mechanics, vibrations, computational uid dynamics (CFD), acoustics, wavelet transforms, and smart materials [1]. Among them, ow past cylinders is often encountered in many elds such as aerospace, nuclear engineering, wind engineering, marine engineering, construction engineering, ground transportation, and sports. VIV should be avoided in most cases of structural safety, as pointed in [2]; the vibrations of small amplitude may lead to fatigue or fretting wear in the long term, while for large amplitude vibration, the damage may occur in the short term. On the other hand, VIV can be used for energy harvesting [3][4][5]. Early experimental results on the characteristics of lift, drag, and dimensionless frequency of vortex shedding (Strouhal number) for the ow around a stationary cylinder, such as Achenbach [6,7], Cantwell and Coles [8], and Schewe [9], have been well summarized in the monographs of Sumer and Fredsoe [10] and Zdravkovich [11,12]. e most interesting phenomena of ow past a cylinder are vortex shedding. Its frequency f vs for a stationary cylinder f vs follows Strouhal's law, i.e., f vs St · U/D, with U being the ow incoming velocity and D being the cylinder diameter. e Reynolds number, defined as Re D � ρUD/μ with ρ and μ being the density and dynamic viscosity of the fluid, respectively, is used to describe the flow regime. In the subcritical regime, with 300 < Re D < 3 × 10 5 , the wake flow is turbulent while the boundary layer over the cylinder surface remains laminar, and St ≈ 0.2 for circular cylinders.
When the cylinder moves in either predetermined forced vibration or flow-induced vibration, the situation will become more complex due to the interactions between vortex motions and cylinder vibrations. e vortex motion may change with the vibration of the cylinder, and vice versa, in the case of VIV. With low incoming velocity, the vortexshedding frequency obeys the stationary-cylinder Strouhal frequency until the reduced velocity U r (its definition is given in the next section) approaches the value of 5. With a further increase in the velocity beyond this point, however, it departs from the Strouhal frequency and begins to follow the natural frequency of the system, remains locked into the natural frequency of the system at U r = 5, and remains locked in until U r reaches the value of about 7 [10].
is phenomenon is known as the "lock-in" or "synchronization" phenomenon, which is characterized by strong fluidstructure interaction.
Although there are different cross-section shapes of cylinders as well as the boundary conditions, the elastically supported circular cylinders are the most fundamental cases. In the study of VIVs of elastically supported cylinders, Feng [13] experimentally investigated the VIV of elastically supported one-degree-of-freedom (1DOF) rigid circular and D-shape cylinders with a large mass ratio m * (ratio of the cylinder mass with displaced fluid mass) in detail in the transverse direction in a wind tunnel and found that when the vortex-induced frequency is close to the natural frequency of the cylinder structure, the cylinder amplitude increases dramatically and the "lock-in" phenomenon will occur, so that the cylinder vibration system can effectively obtain energy from the fluid. Khalak and Williamson [14][15][16][17] measured the VIVs of a cylinder with a low mass and damping ratio. By checking the relationship curve between the maximum vibration amplitude and the reduced velocity, they divided the amplitude response into three regions, namely, the initial branch, the upper branch, and the lower branch. e upper branch is absent in large m * cases, such as in the study by Feng [13]. Khalak and Williamson [16] also found that for the upper and lower branches with large cylinder vibration amplitude, the wake vortex distribution is 2P modes, while the initial branch with small amplitude is 2S modes. In the 2P mode, two pairs of vortices generate in opposite directions in one cycle of cylindrical vibration, while in the 2S mode, there are two independent vortices in one cycle of cylindrical vortex-induced vibration.
In the numerical study of VIVs of elastically supported cylinders, Bahmani and Akbari [18] investigated the dynamic response and vortex shedding of an elastically mounted circular cylinder with Re D ranging from 80 to 160 and found that for m * � 74.5 and 149, lock-in had occurred but not in the case of m * � 298. Only the 2S vortex shedding mode was observed at these low Reynolds numbers, in contrast to the 2P mode of vortex shedding in the turbulent flow regime during lock-in. Li et al. [19], referring to the experimental results of Khalak and Williamson [14], carried out two-dimensional flow field calculations with K-ε and shear stress transport (SST) K-ω turbulence models and found that the SST K-ω model is more appropriate for the simulation of VIV of the elastically mounted rigid cylinder, comparing the cylinder vibration amplitudes with the experimental data, and can successfully capture the 2P wake mode in the upper branch. Guilmineau and Queutey [20] simulated the VIV experiments conducted by Khalak and Williamson [14] with SST K-ω turbulence models with Re D in the range of 900-15 000 and U r between 1.0 and 17.0 and discussed the influence of initial conditions on predicting the maximum amplitude correctly. Results show that under the rest and decreasing velocity conditions, the simulations successfully predict only the lower branch, while with the increasing velocity condition, the maximum vibration amplitude corresponds to the experimental value, although the upper branch does not match experimental data. Pan et al. [21] also employed the two-dimensional numerical simulations of the VIV of an elastically supported cylinder with an SST K-ω turbulence model and discussed the influence of random disturbance on VIV characteristics between experimental and numerical results, and they found different 2P modes reproduced in the simulations were consistent with the experimental results in the upper and lower branches. Khan et al. [22] also made similar simulations of the experimental results of Khalak and Williamson [14] and found that the cylinder vibration amplitude and drag forces are small in the initial branch, and the maximum displacement response of the cylinder is located in the upper branch, with the corresponding dimensionless velocity, which is greater than the experimental results.
Sayyaadi and Motekallem [23] carried out a three-dimensional simulation of the vortex-induced transverse and inline vibrations of three elastic cylinders in a regular triangle arrangement at Re D > 4 × 10 5 and found that the inline amplitude of downstream cylinders could be as large as 0.54D, while the maximum transverse amplitude of three cylinders can reach up to 2.30D. Kinaci [24] argued that 2D simulations provide a fast solution to such a complex flow phenomenon although it is not possible to capture the tip effects of VIV of the elastically supported cylinder and found that the maximum amplitude achieved by the cylinder in the upper branch was not captured with CFD with the K-ω turbulence model as was also mentioned in [20,21]. Recently, Mutlu et al. [25] also numerically investigated the VIV of a smooth cylinder for the transition of shear layer type 2 flow using the SST K-ω turbulence model for the 2D approach and practically discussed the effects of four different mass-damping ratios on the vibration amplitude of the one-degree-of-freedom of the cylinder reveal. ey found that there is an inverse correlation between the massdamping ratios and the transition from the upper to the lower branch.
As pointed out in [25], many past studies dealt with the VIV of flow past a cylinder with high or moderate massdamping ratios, but only a few studies have considered the effects of a low mass-damping ratio so far. erefore, it is still an open topic in this area, especially in the high Reynolds range, such as the influences of the cylinder vibration on the lift and drag forces in each branch and the formation mechanisms of the different vortex modes. e present work numerically studied the vortex induced transverse vibration of an elastically supported cylinder with a very low massdamping ratio, i.e., the mass ratio m * � 2.4, the damping ratio ζ � 0.0054, and the Reynolds number Re D ranged from 2.74 × 10 4 to 1.43 × 10 5 . Experimental data of the VIVs obtained in the cross-flow direction of an elastically supported cylinder with low mass and damping ratios can be found in the study by Khalak and Williamson [14][15][16][17], so that the validation of the present study can be realized by comparing the numerical results with the experimental ones. e unsteady Reynolds-averaged Navier-Stokes equations (URANS) with the SST k-ω turbulence model approach was employed to solving the flow fields, and dynamic mesh technique, specifically the smoothing technique of mesh displacement diffusion method, was used to update the computational grids according to cylinder motion response. e free vibration of the cylinder was simulated first to check the numerical accuracy and the influence of the added mass on the natural frequency. e characteristics of the amplitude response of the cylinder with the influence from the added mass, the vortex shedding frequency, and the fluctuating lift and drag forces at different reduced velocities are explored. Detailed vorticity distributions are examined to show the mechanisms of the 2P and 2S modes.

Mathematical Model.
e elastically supported 1DOF circular cylinder is schematically shown in Figure 1 with the cylinder mass m, the spring coefficient k, and the damping coefficient c. e cylinder diameter D = 0.1 m. e mass ratio m * � m s /(ρπD 2 /4) � 2.4, which is equal to the experiment parameter in [14]. e structural equation of vibration of the elastically mounted cylinder, which allows it to vibrate in the transverse y direction is where F L (y, _ y, t) is the aerodynamic lift force, which is related to the cylinder transverse displacement y and velocity _ y, C L is the lift force coefficient, U is the incoming flow velocity, and ρ is the density of the fluid. Lift force which causes cylinder motion can be obtained based on pressure distribution integration over cylinder surface. e mathematical method of solution of differential equations (1) to obtain a cylinder displacement value in y direction during each time step is presented in Section 2.2. e natural frequency of the cylinder vibration system without the damping in vacuum is f 0 � ��� � k/m √ /2π. In this work, we tuned k and m so that f 0 � 15Hz. If the added mass m A � ρπD 2 /4 using the potential flow theory is taken into account, then the natural frequency of the vibration system is f N � ���������� k/(m + m A )/2π. It can be seen that with a small value of the mass ratio m * , the added mass m A has a great influence on the natural frequency of the vibration system. Specifically, with m * � 2.4, f N � 12.6 Hz here. e dimensionless damping ratio ζ � 0.0054 is selected according to the experiment parameter value in [26].
With the vortex shedding frequency of flow past a stationary cylinder f vs , the Strouhal number as the nondimensional vortex shedding frequency is defined as and the nondimensional reduced velocity U r is defined as follows: In the study, the characteristic Mach number value corresponding to the incoming flow velocity is much less than 0.3, so that the flow can be considered as incompressible.
e continuity equation and momentum equation of incompressible flow can form a closed problem mathematically without adding and solving the energy equation. e characteristic Reynolds number Re D in the simulations ranges from 2.74 × 10 4 to 1.43 × 10 5 so that flow is in turbulent regime. e unsteady RANS equations with the SST K-ω turbulence model by Menter [27] were solved.
e Reynolds averaged continuity equation and momentum equation (Navier-Stokes equation) are as follows [27]: e overbar is the ensemble average procedure in the URANS approach. In the equations (4) and (5) u is the fluid velocity vector and τ R is the turbulent or Reynolds stress tensor. By the Boussinesq assumption, with S the mean strain rate tensor, μ T the turbulent eddy viscosity, K the turbulent kinetic energy, and I the identity tensor. Shock and Vibration 3 e transport equations for the turbulent kinetic energy K and turbulent specific dissipation rate ω in the SST model are described as where G K represents the generation of turbulence kinetic energy due to mean velocity gradients, c μ , σ K , α, β, σ ω , and σ ω,2 are the model tuning parameters, and F 1 the blending function. e turbulent viscosity μ t is finally defined as follows: where a 1 is another model constant, S � ������� 2(S ij S ij ) is the second invariant of the strain rate tensor, and F 2 is a second blending function. Detailed derivation can be found in [27].

CFD Model of the Cylinder VIV System.
In consideration of the high computational cost of 3D simulations, 2D calculations were chosen for the present study as in [19,21,22]. Pan et al. [21] argued that structural vibration produces a high degree of span-wise correlation, giving rise to nearly 2D flow structures for an elastically mounted cylinder subjected to flow, so that the application of 2D numerical simulations are reasonable. Kinaci [24] also suggests that although it is not possible to capture the tip effects with 2D simulations; nevertheless, it provides a fast and acceptable solution to such a complex flow phenomenon. In 2D simulations, the flow is assumed to be uniform along the spanwise direction (independent of z), so that the vortex shedding is uniform along the z direction. e dimensions of the computational domain are shown in Figure 1 with 36D in the cross-flow direction, 18D upstream, and 60D downstream of the cylinder. e domain is larger than that in [19,21,22] in order to reduce the influence of the side boundaries. e calculations were performed in the ANSYS Fluent CFD package, wherein the finite volume method was used to discretize the Reynolds averaged equations (4) and (5) and the turbulent model equations (7) and (8) for K and ω parameters.
e pressure-based solver is used with the pressure-velocity coupling of the SIMPLEC (SIMPLE consistent) scheme, which is an improved SIMPLE (semiimplicit method for the pressure linked equations) algorithm [28]. For the spatial discretization, the gradient terms were discretized with the least square cell-based method, and the pressure equation was discretized with a secondorder scheme. e convective terms in the momentum equation and the turbulent model K and ω equations (7) and (8) were discretized with the second-order upwind technique. e temporal discretization is second-order upwind.
e velocity boundary condition is applied at the inlet boundary, representing the incoming uniform flow as in a wind tunnel. A low turbulent intensity of 1% with a turbulent viscosity ratio of 2% is set at the inlet. An average static pressure of 0 Pa is applied at the outlet boundary. Nonslip wall is set at the cylinder wall, and free slippery wall is applied at the upper and lower sides of the domain.
Quadrilateral mesh cells were used to mesh the computational domain with fine near-wall cell resolution. e mesh independence was checked by increasing the number of cells gradually according to the suggestion in [29], and finally the mesh scheme with a number of 255 102 cells was adopted. Four mesh schemes with element numbers of 64 746, 141 372, 255 102, and 430 192, respectively, were tested. e evolution of cylinder response and the lift force exerted by the fluid with U r � 8 are shown in Figure 2. As it can be seen, the results of the coarse mesh with 64 746 cells have noticeable differences, and the results of 141 372 cells become close to those of the finer cells. e results obtained on grid schemes with 255 102 and 430 192 cells are very close to each other, especially in the steady-state stage after 0.5 s. If the peak-peak value of the cylinder motions in the steadystate stage as the assessment has been taken, and regarding the value from the finest mesh scheme as the precise value, the relative errors of mesh schemes with 64 746, 141 372, 255, and 102 cells are 10.46%, 5.39%, and 1.69%, respectively. So it can be concluded that 255 102 cells are sufficient to obtain the mesh-independent numerical results.
A snapshot of the mesh around the cylinder and in the wake region is shown in Figure 3. e first layer height near the cylinder wall is set to be 10 −5 D and the cell row heights increase gradually with a ratio of 1.02 to ensure that the nondimensional wall distance y+ is around 0.01 and very fine resolution of the near wall flow fields can be obtained. e fine mesh was also created in the wake region along the streamwise direction with a cross-flow distance greater than 4D, considering the vortex motions therein. e time step Δt is set to 10 −5 s, considering the minimum size of the first cell level near the cylinder along the flow direction Δx. For the maximum inlet velocity in the simulation, the Courant-Friedrichs-Lewy condition U max Δt/Δx < 1 is satisfied.
Due to the cylinder motion, the fluid domain changes and the dynamic grid technique is used to tackle the issue. e boundaries of the domain were set at fixed limits. To apply dynamic mesh technology, a cylinder was defined as a rigid body which was allowed to move in a transverse direction.
e fluid domain was specified as the deforming zone. Mesh deformation is determined by unprescribed cylinder motion based on lift force change, which acts on the cylinder. e aerodynamic force on the cylinder was obtained by the integration of the stress over the cylinder wall, and the transverse force (in the y-direction) is the exciting force to the force vibration of the cylinder, as expressed in equation (1). e forced vibration of the cylinder was numerically solved with the Adams-Moulton method [30]. Second-order equation (1) was transformed into a system of two first-order ordinary differential equations: 4 Shock and Vibration For a first-order differential equation _ y � f(t, y), the numerical solution is approximated by where f i is the value of f(t, y) at (t i , y i ). Cylinder transverse displacement value was calculated on each time step and transferred to ANSYS Fluent using the incorporated built-in User Defined Function (UDF) mechanism. Based on the calculated cylinder displacement value, Fluent performed the update for the dynamic mesh for each time step. For the solution of mesh displacement, the diffusion method of smoothing was chosen. It was selected because there is simple geometry, linear motion of the rigid body, and a relatively small amplitude of vibrations of the cylinder. e motion of the mesh nodes in the diffusion method is based on the solution of the diffusion equation.
where u m is mesh displacement velocity and c is the diffusion coefficient, which is determined by the reciprocal of the normalized boundary distance. Based on this equation, a calculation of the mesh cell velocity is performed (velocity at the boundary nodes of the cylinder is used as the Dirichlet boundary condition), and the node positions are then updated according to x new � x old + u m Δt. e simulations started with the cylinder at the equilibrium positions, and once the vortex shedding occurred, the cylinder began to move perceptibly. After the calculation  Shock and Vibration became statistically steady state with the cylinder motion in a periodic manner, the flow data, such as the lift and drag coefficients, and the cylinder displacement were recorded.

e Natural Frequency of the Cylinder Vibration in the Stationary Fluid.
First, the inlet velocity was set small enough so the fluid could be considered stationary. e cylinder is displaced away from the equilibrium position to generate the free vibration. e displacement response is shown in Figure 4, and after examining the period of the free vibration, a natural frequency of 12.2 Hz was found, which corresponds tof N � ���������� k/(m + m A )/2π. By checking the logarithmic decrement δ from Figure 4, and using the equation δ ≈ 2πζ for small damping, we obtained the actual damping ratio ζ act � 0.006, which is a little larger than the predefined experimental value ζ � 0.0054. is is attributed to the fluid viscous damping, i.e., the fluid viscous force over the cylinder behaves as the damping force. In other words, the actual damping coefficient cannot be set a priori since it depends on the fluid viscosity as well as the cylinder vibration.

VIV of the Cylinder.
A Fast Fourier transform was applied to the time series data of the cylinder vibration responses. e peak values in the cylinder response spectra with the corresponding frequencies were picked up to analyze the cylinder response characteristics. Figure 5 shows the amplitude response with a comparison of Khalak and Williamson experimental data [14]. As it can be seen, the numerical results agree well with the experimental data in the initial and lower branches. For the upper branch, the numerical simulation predicted a smaller amplitude with A * max � 0.66 obtained, while the flow pattern is similar to the experimental measurement as explained in Section 3.3. e highest amplitude ratio captured is A * � 0.67 at U r � 5 which was reported by Mutlu et al. [25], claiming the results seem better than those from [21] wherein all values of A * were less than 0.6 with one A * � 0.7 as exception. Kinaci [24] and Mutlu et al. [25] argued that the URANS method cannot capture the randomness of the vortex-induced vibration in the upper branch due to its implemented averaging algorithm. Figure 6 shows the displacement response spectra with different reduced velocities. It can be seen clearly that the vibration amplitude increases with the incoming velocity becoming larger, especially in the region of the cylinder vibration resonance where lock-in phenomena occur, and decreases abruptly after a certain value of the reduced velocity Ur when the lock-in disappears. In the case of small Ur, the cylinder vibration spectra have two obvious peaks: one is related to the vortex shedding frequency (excitation frequency), and the other is the natural frequency with the added mass f N . With increasing Ur, the two peaks merge into a single peak with the frequency around the natural     Shock and Vibration frequency in vacuum f 0 . And with a further increase of Ur, the cylinder vibration amplitude decreases dramatically, and the corresponding frequency is the vortex shedding frequency f vs . A Fast Fourier transform was applied to the time history data of lift and drag coefficients also, and the spectra are shown in Figures 7 and 8, respectively. e abscissa range of Figure 8 is set to 0.1-10, which is different from 0-5 in Figure 7, considering that the main drag force fluctuations' frequency is two times that of the lift force fluctuations. Also, the drag force has a very large nonzero component at 0 Hz. e frequencies relating to the peak values in the lift spectra are the vortex shedding frequencies since each vortex shedding causes a lift pulsation. Figure 9 shows the overall fluctuations of lift and drag forces versus the reduced velocity Ur, with the root mean square (RMS) value as the indicator. e mean values of the drag coefficient are also shown in Figure 9, while the mean values of the lift coefficient are absent since they are sufficiently low and tend to zero due to statically symmetric flow distribution about the horizontal line through the cylinder equilibrium position. e vortex shedding frequencies are presented in Figure 10 with the comparison of Khalak and Williamson's experimental data [14] therein. For Ur � 10, the peaks at the resonance and at the vortex shedding frequency are    Shock and Vibration 7 comparable, so both frequencies are presented in Figure 10. It can be seen clearly that frequency synchronization (lockin) occurs in the upper and lower branches. e phase angle differences Δφ between the cylinder displacement responses and the lift force at the peak value frequencies are obtained by checking the Fourier spectra.
e results are plotted in Figure 11. e conspicuous feature in Figure 11 is the large jump of Δφ when Ur goes through the upper branch, which was found in experimental studies [10,13,26]. It shows the vortex mode transition behind the    Shock and Vibration cylinder accompanied by the phase jump as discussed in Section 3.3.

Vortex Patterns.
e vorticity plots at the initial branch are shown in Figure 12, with the red color vortex (positive value) as the counter-clockwise vortex and the blue one as the clockwise vortex. e vortex pattern behind the cylinder is 2S mode, meaning two single vortices are formed per cycle of the cylinder vibration. e agreement between the numerical and experimental results is reasonable in Figure 12, although the mass ratios are different. is may be due to the fact that at low mass ratios, the cylinder vibration responses in these two cases are comparable, and also the flow regimes are similar with the Reynolds number Re D less than the critical value. As can be seen, when the cylinder is in the top position, the upper side counter-clockwise vortex around the cylinder is in the process of cutting the down-side clockwise vortex. When the cylinder passes across the equilibrium position from the upper peak position, the clockwise vortex has nearly a point of separation from its boundary layer. Once the separation is complete, a free vortex is formed, and then it is convected downstream by the flow, creating the 2S vortex street. e vorticity contours of the upper branch are shown in Figure 13, with good agreement with the experimental data. e vortex pattern behind the cylinder is 2P mode, i.e., two pairs of vortices are formed per cycle. As can be seen, in the separation of the down side clockwise vortex, the up side vortex is itself divided into two parts, and one part, together with the free vortex, constitutes the vortex pair. e mechanism of the 2P mode is different from the 2S mode, and this is due to the large cylinder vibration amplitude causing the cutting vortex to be less stable than in the 2S mode case. e vortex plots of the lower branch in Figure 14 also show the 2P mode; however, there is some difference between the vortex modes of the upper and lower branches. e vortices in the vortex pair in the upper branch are much closer to each other, and one vortex is fairly smaller than the other, while in the lower branch, the two vortices are separated by a noticeable distance and the sizes are competitive.

Discussion
In the present study, the vibration response of the elastically supported cylinder in the cross-flow direction due to VIV with a 2D URANS approach was numerically investigated. From the free vibration result, as judged from Figure 4, the added mass of the fluid is precipitable because the natural frequency is consistent with the equation . From the cylinder displacement spectra in Figure 6, at small reduced velocity Ur, the vortex shedding frequency f vs is away from the natural frequency f N of the cylinder vibration system. e cylinder response has both f vs and f N components. With the increase of Ur, f vs becomes close to f N , and the cylinder vibrations become much more vehement. In the case of large cylinder amplitude, the flow around the cylinder certainly cannot be considered as potential flow. In other words, the added mass from the potential flow theory may not be applicable. As can be seen clearly in Figure 6, there is only one peak in the vibration spectra, whose frequency is the natural frequency in vacuum of f 0 . Further increasing Ur, the cylinder vibration amplitude decreases dramatically with a small peak at the vortex shedding frequency f vs .
From the amplitude response Figure 5, the numerical results agree well with the experimental data by Khalak and Williamson experimental data [10], especially the initial and lower branches. For the upper branch, numerical simulation obtained smaller amplitudes. e same result was also found in other numerical studies with 2D approaches, such as [19,21,22], except in Guilmineau and Queutey [20], a sufficiently large amplitude was obtained in a very narrow Ur range using very small steps of Ur. Kinaci [24] and Mutlu et al. [25] explained that the URANS method cannot capture the nature of the randomness of the vortex-induced vibration in the upper branch due to its implemented averaging algorithm. Nevertheless, in the present study, the flow vorticity distribution behind the cylinder is close to the experimental measurements as shown in Figure 13. e 2P vortex modes in the upper and lower branches by URANS and SST K-ω turbulence model are captured successfully, and the vorticity contour plots are consistent with experimental results. e vorticity distribution of the upper branch shows a more fragmental pattern with the much closer vortices of the 2P pair, while the lower branch has a more connected version of the vortex. e lift coefficient spectra in Figure 7 show that for low Ur the vortex shedding component f vs is obviously high. Increasing Ur the conspicuous peak is located at f 0 indicating a "lock-in" phenomenon, and further increase of Ur results in the peak frequency at f vs . For Ur � 4, which is in the initial branch and close to the upper branch, the lift and drag spectra are very complex. As can be seen, the lift has a dominating component at the vortex shedding frequency f vs ., which is different from the cylinder vibration frequency f 0 , showing lock-in has not occurred yet. By comparing Figure 6, it can be seen that the cylinder vibration amplitude is quite large at f 0 . is can induce a certain level of lift force fluctuation at the same frequency. However, the drag force has a prominent peak at the natural frequency of the cylinder vibration system, although some harmonics related to vortex shedding frequency are also noticeable. For the upper and lower branches, lift and drag force spectra have peaks at f 0 and 2f 0, respectively, indicating the frequency lock-in or synchronization happens. For Ur � 10, close to the desynchronization, the vortex shedding frequency component in the lift force spectra is also observable. For Ur � 14, in the desynchronization due to the small cylinder vibration, the lift and drag forces have only the peaks at the f vs and two times of f vs , respectively, showing that the interactions between the vortex shedding and the cylinder motion are negligible.
As it can be seen from Figure 9, the lift force fluctuation increases as Ur approaches the upper branch from the initial branch and decreases gradually when Ur goes down the lower branch to the desynchronization region. While for the drag force, its fluctuation is quite large at the initial branch, and experiences a sharp decline when Ur fully enters the upper branch. Analyzing the lift force coefficient spectra, it is known that in the initial branch, both the vortex shedding and the cylinder vibration have a great influence on the drag force, leading to large fluctuations with different frequencies. However, in the upper and lower branches, the effect of the cylinder vibration on the drag fluctuations is superior. e mean values of the drag coefficient are also shown in Figure 9, and it is interesting to find that its pattern is similar to the RMS lift; i.e., it increases in the upper branch and then decreases in the lower branch and desynchronization region. e phase angle difference Δφ between the lift force and the cylinder vibration clearly shows the jump occurring when Ur enters the upper branch from the initial branch, shown in Figure 11. e power input to the cylinder vibration from the fluid flow is closely related to Δφ. e power input equals the product of the lift force and the cylinder velocity in the cross-flow direction, F L _ y. For simple harmonic vibration, the power input value depends on cos (Δφ − 90 ∘ ) considering the velocity phase angle leads the displacement phase angle by 90°.
It can be seen from Figure 11, in the initial branch region, there is a small phase difference between the lift force and the cylinder displacement, which means cos (Δφ − 90 ∘ ) has small values, so that the energy transition from the fluid into the cylinder vibration system is small enough. When Ur goes into the upper and lower branches, Δφ has certain distance to 180°, so that cos (Δφ − 90 ∘ ) increases substantially, which means fluid energy can be transferred into the vibration system more easily to boost the cylinder vibration. In the desynchronization region, Δφ is very close to 180°, resulting in a sufficiently small value of cos (Δφ − 90 ∘ ), so that the cylinder vibration amplitude decreases dramatically.

Conclusions
e vortex-induced vibration (VIV) of an elastically supported cylinder was numerically studied by computational fluid dynamics (CFD). e unsteady Reynolds averaged Navier-Stokes equations (URANS) with the turbulence modelling of shear stress transport (SST) K-ω was applied to calculate the flow fields, and the dynamic mesh approach was used to tackle the domain change due to the cylinder motion. e characteristics of the amplitude response of the cylinder and the vortex shedding frequency at different reduced velocity were compared with the experimental data, and good agreement was found. e added mass was perceptible in the free vibration and in the initial branch with small cylinder vibration amplitude, while in the upper and lower branches the influence of added mass is weak regarding the fact that the main response frequency of the cylinder is close to the natural frequency in vacuum. Detailed vorticity distributions clearly show the 2P vortex modes at the upper and lower branches and 2S mode in the initial branch.
In the initial branch and close to the upper branch, the lift and drag spectra are very complex with peaks related to both vortex shedding frequency f vs and cylinder dominating vibration frequency f 0 . For the upper and lower branches, lift and drag forces have peaks at f 0 and 2 f 0 indicating the frequency lock-in or synchronization happens. In the desynchronization region, which is characterized by the small cylinder vibration, the lift and drag forces have only their peak values at the f vs and 2f vs , respectively, showing that the interactions between the vortex shedding and the cylinder motion are negligible. e root mean square (RMS) value of lift force increases as the reduced velocity Ur approaches the upper branch from the initial branch and decreases gradually when Ur goes down the lower branch to the desynchronization region. However, for the drag force, its RMS value is quite large at the initial branch, and experiences a sharp decline when Ur fully enters the upper branch. e pattern of the changes in the mean drag coefficient is similar to that of the RMS lift, i.e., it increases in the upper branch and then decreases in the lower branch and desynchronization region. e phase difference between the lift force and the cylinder vibration, which influences the power transfer from the flow to the cylinder vibration, was analyzed, and the phase jump was found from the initial branch to the upper branch with the vortex mode switching from 2S to 2P. e difference between the generation mechanisms of the 2S and 2P vortex modes lies in whether the cutting vortex on one side of the cylinder separates into two parts or not. With the large cylinder vibration, the cutting vortex becomes two parts, one part together with the free vortex to form the vortex pair. e 2P modes in the upper and lower branches show some differences. For the upper branch, the vortices in the vortex pair are much closer to each other with one very small vortex, while for the lower branch, the two vortices are separated by a noticeable distance with comparable vortex sizes.
Data Availability e data that supports the findings of this study are available within the article.
Shock and Vibration 11