Very Large-Eddy Simulations of the Flow Past an Oscillating Cylinder at a Subcritical Reynolds Number

: This work focuses on ﬂow past a circular cylinder at a subcritical Reynolds number. Although this classical study has been a concern for many years, it is still a challenging task due to the complexity of ﬂow characteristics. In this paper, a high-e ﬃ ciency very large-eddy simulation method is adopted and veriﬁed in order to handle the oscillating boundary. A series of numerical simulations are conducted to investigate the transient ﬂow around the oscillating cylinder. The results show that the vortex shedding mode varies with an increase in the excitation amplitude and the excitation frequency. Vortex shedding is a lasting process under the condition of a low excitation amplitude that leads to irregular ﬂuctuations of the lift and drag coe ﬃ cients. For a vortex shedding mode that exhibits a strong vortex pair and a weak vortex pair or a weak single vortex, the temporal evolution of the lift coe ﬃ cient of the oscillating cylinder shows irregular ”jumping” at a speciﬁc time per cycle corresponding to the shedding of the strong vortex pair. The vortex shedding mode and the frequency and time of the vortex shedding co-determine the temporal evolutions of the lift and drag coe ﬃ cient.


Introduction
As the exploitation of offshore oil and gas expands into deep water and regions with strong ocean currents, the significance of vortex-induced vibration (VIV) has attracted the attention of engineers. However, some grave engineering problems are due mainly to VIV, similar to fatigue failure on submarine pipelines. Therefore, the prediction of the hydrodynamic forces on an oscillating cylinder is still paramount to the design of the structures. Numerical and experimental studies on the flow around a cylinder oscillating in the fluid have always been the focus of related studies for many years. However, it is still a significant primary task because of the complexity of flow characteristics. At present, one of the common methods to determine the VIV phenomenon is a circular cylinder that is forced to oscillate with a preset motion including an excitation amplitude (A) and a frequency (f e ) approximating VIV. Given the richness of the research on VIV, only a forced oscillating cylinder problem is included in the following discussion. Bishop et al. conducted an experimental study on the hydrodynamic forces acting on a circular cylinder and found that an abrupt jump existed in the phase angle between the cylinder motion and the lift force when the oscillation frequency (f e ) approximately equaled the natural shedding frequency (f s ) of the stationary cylinder [1]. In addition, the average amplitude of the lift and drag increased when the frequency ratio satisfied f e /f s ≈ 1. This condition is also known as a "lock-in" phenomenon, which signifies that the vortex shedding frequency of the oscillating cylinder (f v ) primarily depends on the oscillation frequency f e , i.e., f v = f e . Additionally, it also means that the vortex shedding frequency f v is no longer determined by the Strouhal number (St).
In this paper, a high-efficiency VLES method is used to simulate the VIV phenomenon of the oscillating cylinder at a subcritical Reynolds number Re = 4000 combined with the dynamic grid technique.

Computational Method and Feasibility Analysis
VLES was proposed in 1998 as a hybrid turbulence method that combined the advantages of diverse turbulence modelings of RANS (Reynolds-averaged Navier-Stokes), LES, and DNS. [18][19][20] It is also known as the flow simulation methodology (FSM) which exhibits significant potential in some application. This method essentially depends on the resolution control function Fr to damp the Reynolds stresses, aimed at modeling the subscale stress tensor, which is expressed as: In this paper, a specific VLES method established on the standard k-ε RANS turbulence model is adopted, whose essence is to modify the formulation of the turbulent viscosity by using a control function Fr [18]: In this VLES method, the transport equations are expressed as: where these transport equations are precisely the same as in the standard k-ε model, k is the modeled turbulent kinetic energy, and ε is the turbulent dissipation rate. In addition, P k is the kinetic energy production item and is defined as: In a computational domain, the distribution of the Fr value depends on the variation of two turbulence length scales and constant terms in the Fr expression which can realize the transformation of the calculation model between LES and RANS. The Fr is defined as: where β is the model constant and equals 2.0 × 10 −3 . In Equation (7), Fr is approximately understood as the ratio of the unsolved turbulent kinetic energy to the total turbulent kinetic energy, and therefore it is between 0 and 1. Here, L i is the integral length scale, Lc is the turbulent cutoff length scale, and L k is the Kolmogorov length scale, which are defined as, respectively: Appl. Sci. 2020, 10, 1870 4 of 18 where (∆x∆y∆z) 1/3 is the equivalent length of a mesh resolution, n = 2, Cx = 0.61, and ν is the kinematic viscosity of the fluid. In the regions near the wall, Fr is close to 1, due to L c > L i , meaning the RANS model is approximately adopted. However, a VLES model is applied in the region away from the wall, due to Fr < 1, or L c < L i . The VLES method is through a resolution control function Fr to determine the computational model, therefore, the user defined function (UDF) in ANSYS Fluent is required to realize an intended function. The VLES method has both the accuracy of the LES model for calculating the large-scale vortex structures and the high efficiency of the RANS model for calculating the flows near the wall. More importantly, it also has good predictability for a coarse mesh. Xiong et al. introduced the VLES method into the dynamic grid technology to calculate the wake structure of an oscillating hydrofoil, and verified the feasibility and high efficiency of the VLES method applied in the dynamic grid technology by comparisons with LES and RANS models [20]. A study by [18] also provided more details about the advantages of the VLES method for the flow separation simulation.

Computational Model and Sensitivity Study
A grid sensitivity study is conducted by employing four meshes to analyze the flow past an oscillation cylinder at A/D = 1 and f e /f s = 0.75. The cylinder oscillates as a harmonic mode, expressed as: y(t) = Asin(2πf e t). Computational domain is discretized by a boundary layer and an external unstructured grid, and the boundary layer oscillates as a rigid body along with the boundary of the cylinder. The external unstructured grid is updated by a spring-based smoothing and a remeshing method to avoid a large mesh deformation. The detailed nodes and the total number of elements are shown in Table 1. Table 2 shows the comparison between the Strouhal number (St = f s D/u) obtained by the four meshes and those reported by previous works. The St numbers from the former are all around 0.22 and increase slightly with an increase in the total number of elements. Its growth rate is not more than 2.26%. In the case of Mesh-III and Mesh-IV, the St numbers exhibit a tiny change and tend to a stable state, which indicates that the calculated results are insensitive to a further increase in the number of elements. Additionally, the results also verify the conclusion that the VLES method has good predictability for a coarse grid, and high efficiency for dynamic boundary problems.  [16,23]. Additionally, Dong et al. also studied the flow structures on the onset and development of the shear layer instability at Re = 3900 and 4000 and Re = 10,000. The results showed that the velocity fluctuations and the Reynolds stress from the particle image velocimetry (PIV) results at Re = 4000 agreed well with the simulation studies at Re = 3900 using the DNS method, illustrating that the fluid mechanics between the two had exactly similar characteristics. Thereby, the St number from experimental data at Re = 3900 is reasonably selected to validate the simulation results at Re = 4000 in this section.
Combined with the above analysis, the St number of a 3D stationary cylinder at Re = 4000 should be 0.2 to 0.22 and tends toward 0.22. In this section, the St numbers from the sensitivity study are around 0.22, which also verify the validation of simulation results. The reason for the small difference in the St numbers could be caused by the adoption of a 2D model, thus, deviating from the reported results of the 3D model from Ref. [16,[21][22][23], because the 2D simulation does not capture the 3D effects. However, several 2D studies verified that the computed results are quite consistent with the experimental results, despite this limitation. [10,11,16,23] Moreover, Parnaudeau et al. [16] and Dong et al. [23] studied the flow over a circular cylinder using a statistical estimation, which was shown to need large integration time for most flow characteristics resulting in an uncertainty of about 10% in their experiments. Therefore, numerical results are reliable, and Mesh-III is finally adopted to calculate the flow around the lateral forced oscillation of the cylinder in all simulations. The specific computational domain and the mesh near the cylinder boundary are shown in Figure 1. and Dong et al. found that the St number at Re = 3900 was within the range 0.2-0.22 [16,23]. Additionally, Dong et al. also studied the flow structures on the onset and development of the shear layer instability at Re = 3900 and 4000 and Re = 10,000. The results showed that the velocity fluctuations and the Reynolds stress from the particle image velocimetry (PIV) results at Re = 4000 agreed well with the simulation studies at Re = 3900 using the DNS method, illustrating that the fluid mechanics between the two had exactly similar characteristics. Thereby, the St number from experimental data at Re = 3900 is reasonably selected to validate the simulation results at Re = 4000 in this section.
Combined with the above analysis, the St number of a 3D stationary cylinder at Re = 4000 should be 0.2 to 0.22 and tends toward 0.22. In this section, the St numbers from the sensitivity study are around 0.22, which also verify the validation of simulation results. The reason for the small difference in the St numbers could be caused by the adoption of a 2D model, thus, deviating from the reported results of the 3D model from Ref. [16,[21][22][23], because the 2D simulation does not capture the 3D effects. However, several 2D studies verified that the computed results are quite consistent with the experimental results, despite this limitation. [10,11,16,23] Moreover, Parnaudeau et al. [16] and Dong et al. [23] studied the flow over a circular cylinder using a statistical estimation, which was shown to need large integration time for most flow characteristics resulting in an uncertainty of about 10% in their experiments. Therefore, numerical results are reliable, and Mesh-III is finally adopted to calculate the flow around the lateral forced oscillation of the cylinder in all simulations. The specific computational domain and the mesh near the cylinder boundary are shown in Figure 1.

Verification of the VLES Model
To verify the feasibility of the VLES method to address the oscillating cylinder, we simulated a canonical case of a cylinder oscillating in a fluid initially at rest. We choose the oscillation cylinder adopting the same parameters as reported by Dustch et al., such as Re = 100, KC = 5, U = 0.01 m/s, D = 0.01 m, fe = 0.2 Hz, where KC is a Keulegan-Carpenter number [24]. The cylinder is given by a harmonic motion oscillating in the horizontal direction, and its translational displacement satisfies X(t) = −Aesin(2πfet), where Ae is the oscillating amplitude. The cylinder initially moves towards the negative x coordinate with a velocity UX = −2πfeAe = −U, UY = 0. UX and UY are the velocity components in the x and y direction, respectively. The cylinder surface is set to satisfy the non-slip boundary condition. Similarly, the studied cylinder in the following simulation initially oscillates towards the positive y coordinate with a velocity Uy = 2πfA, Ux = 0, the cylinder surface is also set to meet the nonslip boundary condition. Ux and Uy are the velocity components in the x and y direction, respectively.

Verification of the VLES Model
To verify the feasibility of the VLES method to address the oscillating cylinder, we simulated a canonical case of a cylinder oscillating in a fluid initially at rest. We choose the oscillation cylinder adopting the same parameters as reported by Dustch et al., such as Re = 100, KC = 5, U = 0.01 m/s, D = 0.01 m, f e = 0.2 Hz, where KC is a Keulegan-Carpenter number [24]. The cylinder is given by a harmonic motion oscillating in the horizontal direction, and its translational displacement satisfies X(t) = −A e sin(2πf e t), where A e is the oscillating amplitude. The cylinder initially moves towards the negative x coordinate with a velocity U X = −2πf e A e = −U, U Y = 0. U X and U Y are the velocity components in the x and y direction, respectively. The cylinder surface is set to satisfy the non-slip boundary condition. Similarly, the studied cylinder in the following simulation initially oscillates towards the positive y coordinate with a velocity U y = 2πf A, U x = 0, the cylinder surface is also set to meet the non-slip boundary condition. U x and U y are the velocity components in the x and y direction, respectively. Figure 2 shows the time history of the force coefficient (C F = F/(0.5ρu 2 D)) of the oscillation cylinder in the horizontal direction. The resultant force (solid line), the differential pressure component (dashed line, and the viscous force component (dash dot line) are shown, respectively, which are well predicted as compared with the results calculated by Dustch et al. [24] Dustch et al. (two columns on the left) in Figure 3. Combined with the comparisons in Figures 2 and  3, it can be seen that the VLES method can be used to effectively handle the dynamic boundary problems, such as an oscillating cylinder. To further compare the numerical and experimental results, Figure 4 provides local velocity information for three phase angles. The comparison shows a good agreement between them, and therefore the VLES method is applicable to calculate the flow past an oscillating cylinder.   The instantaneous contours of the pressure and vorticity of the fluid around the oscillating cylinder at four different phase angles are shown in Figure 3. The four phase angles (γ = 2πf e t) adopted are γ = 0 • , γ = 96 • , γ = 192 • , and γ = 288 • , respectively. These vorticity contours show the generation and shedding of the wake vortices of the oscillating cylinder. It is evident that the computed results (two columns on the right) using the VLES method are very consistent with results reported by Dustch et al.
(two columns on the left) in Figure 3. Combined with the comparisons in Figures 2 and 3, it can be seen that the VLES method can be used to effectively handle the dynamic boundary problems, such as an oscillating cylinder. To further compare the numerical and experimental results, Figure 4 provides local velocity information for three phase angles. The comparison shows a good agreement between them, and therefore the VLES method is applicable to calculate the flow past an oscillating cylinder.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 19 component (dashed line, and the viscous force component (dash dot line) are shown, respectively, which are well predicted as compared with the results calculated by Dustch et al. [24] The instantaneous contours of the pressure and vorticity of the fluid around the oscillating cylinder at four different phase angles are shown in Figure 3. The four phase angles (γ = 2πfet) adopted are γ = 0°, γ = 96°, γ = 192°, and γ = 288°, respectively. These vorticity contours show the generation and shedding of the wake vortices of the oscillating cylinder. It is evident that the computed results (two columns on the right) using the VLES method are very consistent with results reported by Dustch et al. (two columns on the left) in Figure 3. Combined with the comparisons in Figures 2 and  3, it can be seen that the VLES method can be used to effectively handle the dynamic boundary problems, such as an oscillating cylinder. To further compare the numerical and experimental results, Figure 4 provides local velocity information for three phase angles. The comparison shows a good agreement between them, and therefore the VLES method is applicable to calculate the flow past an oscillating cylinder.       However, when Fr locates between 0 and 1.0, the turbulence is modeled as the VLES model, considered to be a "coarse LES" model, to capture the large-scale vortices located in the regions away from the wall. The Fr contour and the Fr distribution, in Figure 5, show that the Fr values of most of the regions are within the range 0~0.1, meaning these regions adopt the VLES model to solve the Navier-Stokes equations. In addition, regions behind the cylinder have Fr values from 0.1 to 0.4, due to the larger turbulent kinetic energy and moderate grids, and also employ the VLES model. However, the Fr value around the boundary of the cylinder was basically maintained between 0.75 and 1.0, indicating that the RANS model is recalled. The above distribution law of the Fr values is in good agreement with the principle of the VLES method. It can also be seen in Figure 5 that the distribution of Fr is continuous in the computational domain.

Correctness of the VLES Method
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 19 addition, regions behind the cylinder have Fr values from 0.1 to 0.4, due to the larger turbulent kinetic energy and moderate grids, and also employ the VLES model. However, the Fr value around the boundary of the cylinder was basically maintained between 0.75 and 1.0, indicating that the RANS model is recalled. The above distribution law of the Fr values is in good agreement with the principle of the VLES method. It can also be seen in Figure 5 that the distribution of Fr is continuous in the computational domain. In addition, the VLES method had a clear explanation on the y+ distribution that Fr tended to 1 in the premise of y+ < 3.12, which means that the RANS method has been successfully recalled near the wall [18]. As can be seen in Figure 6, y+ on the cylinder surface is guaranteed to be within 2.0 and its average value is around 1.0, which meets the above requirements.

Results and Discussion
To illustrate the VIV characteristics of the oscillating cylinder at Re = 4000, we conducted a series of simulations for A/D = 0.5 to 1.5 and fe/fs = 0.5 to 1.35 in this paper, and the specific parameter settings are shown in Table 3. Additionally, the temporal evolutions of the lift and drag coefficients, and the frequency and mode of the vortex shedding is discussed in detail, respectively.  In addition, the VLES method had a clear explanation on the y+ distribution that Fr tended to 1 in the premise of y+ < 3.12, which means that the RANS method has been successfully recalled near the wall [18]. As can be seen in Figure 6, y+ on the cylinder surface is guaranteed to be within 2.0 and its average value is around 1.0, which meets the above requirements.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 19 addition, regions behind the cylinder have Fr values from 0.1 to 0.4, due to the larger turbulent kinetic energy and moderate grids, and also employ the VLES model. However, the Fr value around the boundary of the cylinder was basically maintained between 0.75 and 1.0, indicating that the RANS model is recalled. The above distribution law of the Fr values is in good agreement with the principle of the VLES method. It can also be seen in Figure 5 that the distribution of Fr is continuous in the computational domain. In addition, the VLES method had a clear explanation on the y+ distribution that Fr tended to 1 in the premise of y+ < 3.12, which means that the RANS method has been successfully recalled near the wall [18]. As can be seen in Figure 6, y+ on the cylinder surface is guaranteed to be within 2.0 and its average value is around 1.0, which meets the above requirements.

Results and Discussion
To illustrate the VIV characteristics of the oscillating cylinder at Re = 4000, we conducted a series of simulations for A/D = 0.5 to 1.5 and fe/fs = 0.5 to 1.35 in this paper, and the specific parameter settings are shown in Table 3. Additionally, the temporal evolutions of the lift and drag coefficients, and the frequency and mode of the vortex shedding is discussed in detail, respectively.

Results and Discussion
To illustrate the VIV characteristics of the oscillating cylinder at Re = 4000, we conducted a series of simulations for A/D = 0.5 to 1.5 and f e /f s = 0.5 to 1.35 in this paper, and the specific parameter settings are shown in Table 3. Additionally, the temporal evolutions of the lift and drag coefficients, and the frequency and mode of the vortex shedding is discussed in detail, respectively.  Figure 7 shows the vorticities around and behind the oscillating cylinder at A/D = 0.5 and f e /f s = 0.5 to 1.35. Figure 8 shows the temporal evolutions of the lift coefficient (Cl) and drag coefficient (Cd) (left) and the power spectral density (PSD) of the lift coefficient of the oscillating cylinder (right) in the corresponding parameters. The PSD distribution of the lift coefficient of an oscillating cylinder is obtained by a fast Fourier transform (FFT).  As per the Parseval's identity, the lift coefficient can be expressed as follows by FFT: The fluid force acting on the cylinder includes two contributing parts based on pressure and shear stress. Additionally, the lift (coefficient) F l (Cl) and drag (coefficient) F d (Cd) are defined as, respectively: and where p(t) and ω z (t) are obtained by a dimensionless transformation of the pressure and the axial vorticity on the surface of the oscillating cylinder, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 19 In Figure 7, the wake vortices are mainly distributed near the wake center line, showing no significant lateral deviation due to less excitation amplitude at A/D = 0.5. In addition, the number of the wake vortices at fe/fs = 0.5 obviously is greater than that at fe/fs = 0.75. This is caused by the fact that the vortex shedding is influenced by the excitation frequency fe and, according to the PSD at A/D = 0.5 and fe/fs = 0.5 in Figure 8, the Stroulhal frequency fs also impacts the vortex shedding. Both of the two frequencies combine to determine the vortex shedding in the wake, causing an obvious "beating" characteristic in the temporal evolutions of the lift and drag coefficients.
Additionally, the vortex shedding mode also verifies that the frequency ratio fe/fs = 0.5 is not In Figure 7, the wake vortices are mainly distributed near the wake center line, showing no significant lateral deviation due to less excitation amplitude at A/D = 0.5. In addition, the number of the wake vortices at f e /f s = 0.5 obviously is greater than that at f e /f s = 0.75. This is caused by the fact that the vortex shedding is influenced by the excitation frequency f e and, according to the PSD at A/D = 0.5 and f e /f s = 0.5 in Figure 8, the Stroulhal frequency f s also impacts the vortex shedding. Both of the two frequencies combine to determine the vortex shedding in the wake, causing an obvious "beating" characteristic in the temporal evolutions of the lift and drag coefficients.
Additionally, the vortex shedding mode also verifies that the frequency ratio f e /f s = 0.5 is not within the range of the "lock-in" phenomenon at Re = 4000. When f e /f s = 0.75, the scale and strength of the vortices shedding in the wake almost maintains the same level. The vortex shedding frequency is only impacted by the oscillation frequency f e as per the PSD at A/D = 0.5 and f e /f s = 0.75 in Figure 8, thereby being in the "lock-in" range. With an increase in the frequency ratio f e /f s , the strength of the vortex shedding also increases. But the vortex shedding mode keeps a 2S mode in the range of 0.75 ≤ f e /f s ≤ 1.35, and the dominant frequency is always the excitation frequency f e . Additionally, the boundary of the cylinder moves faster as the excitation frequency increases, thereby accelerating the nearby fluid more and causing the wake vortices large lateral displacement in downstream. Along with an increase in the oscillating frequency ratio, the amplitudes of the lift and drag coefficients also increase. Additionally, the temporal evolutions of the lift and drag coefficients show different levels of fluctuation and no significant periodicity at f e /f s ≥ 0.95. It could be that the time of the vortex shedding is not consistent per cycle under the condition of small excitation amplitude or it could be that the vortex shedding is a lasting process, leading the lift and drag coefficients to fluctuate irregularly.
Similarly, Figure 9 shows the vorticities around and behind the oscillating cylinder at A/D = 1 and f e /f s = 0.5 to 1.35. Figure 10 shows the temporal evolutions of the lift and drag coefficients (left) and the PSD of the lift coefficient of the oscillating cylinder (right) in the corresponding parameters. However, the vortex shedding mode and the characteristics of the lift and drag coefficients at A/D = 1 are very different from that at A/D = 0.5.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 19 boundary of the cylinder moves faster as the excitation frequency increases, thereby accelerating the nearby fluid more and causing the wake vortices large lateral displacement in downstream. Along with an increase in the oscillating frequency ratio, the amplitudes of the lift and drag coefficients also increase. Additionally, the temporal evolutions of the lift and drag coefficients show different levels of fluctuation and no significant periodicity at fe/fs ≥ 0.95. It could be that the time of the vortex shedding is not consistent per cycle under the condition of small excitation amplitude or it could be that the vortex shedding is a lasting process, leading the lift and drag coefficients to fluctuate irregularly.
Similarly, Figure 9 shows the vorticities around and behind the oscillating cylinder at A/D = 1 and fe/fs = 0.5 to 1.35. Figure 10 shows the temporal evolutions of the lift and drag coefficients (left) and the PSD of the lift coefficient of the oscillating cylinder (right) in the corresponding parameters. However, the vortex shedding mode and the characteristics of the lift and drag coefficients at A/D = 1 are very different from that at A/D = 0.5.
When fe/fs = 0.5, the vortex shedding mode is a similar "P + S" mode, and the lift and drag coefficients also show the "beat" characteristics. However, its vortex shedding mode actually is a 2S mode placed out of the "lock-in" range. This is due to the fact that the vortex shedding is also influenced by the excitation frequency and the Stroulhal frequency at A/D = 1 and fe/fs = 0.5, which is in line with the case at A/D = 0.5 and fe/fs = 0.5. The detachment of the vortices from the cylinder surface is dominated by a joint action of these two frequencies, leading to the complexity of the vortex shedding in the wake. The temporal evolutions of the lift and drag coefficients exhibit significant fluctuations, which correspond to the vortex shedding mode.  When A/D = 1 and fe/fs = 0.75, the vortex shedding exhibits a 2S mode, and the vortices shed behave regularly with a higher dominant frequency (0.08772) than that of its excitation frequency fe (0.06748). Thus, the vortices shed in the wake closely to the condition of the same excitation frequency. Additionally, the temporal evolutions of the lift and drag coefficients also perform periodicity. The case of fe/fs = 0.75 is an exception to the "lock-in" range, which could be due to the fact that it is in an intermediate state from the 2S mode to the 2P0 mode or it could be that the dominant frequency is a fusion between the excitation frequency fe and the Stroulhal frequency fs.
In the case of fe/fs = 0.95, the vortex shedding mode is transformed into a 2P mode. Moreover, the secondary vortex of each vortex pair has lower strength than that of the primary vortex. This vortex shedding mode corresponds to a 2P overlap mode, 2P0, which is consistent with the vortex shedding mode at fe/fs ≈ 1 obtained by Morse and Willamson [3]. Additionally, the vortex shedding mode is within the range of the "lock-in" phenomenon as a stable state and the vortex shedding frequency of the oscillating cylinder approaches the natural shedding frequency, which can lead to the approximate symmetry of the wake vortices. When f e /f s = 0.5, the vortex shedding mode is a similar "P + S" mode, and the lift and drag coefficients also show the "beat" characteristics. However, its vortex shedding mode actually is a 2S mode placed out of the "lock-in" range. This is due to the fact that the vortex shedding is also influenced by the excitation frequency and the Stroulhal frequency at A/D = 1 and f e /f s = 0.5, which is in line with the case at A/D = 0.5 and f e /f s = 0.5. The detachment of the vortices from the cylinder surface is dominated by a joint action of these two frequencies, leading to the complexity of the vortex shedding in the wake. The temporal evolutions of the lift and drag coefficients exhibit significant fluctuations, which correspond to the vortex shedding mode.
When A/D = 1 and f e /f s = 0.75, the vortex shedding exhibits a 2S mode, and the vortices shed behave regularly with a higher dominant frequency (0.08772) than that of its excitation frequency f e (0.06748). Thus, the vortices shed in the wake closely to the condition of the same excitation frequency. Additionally, the temporal evolutions of the lift and drag coefficients also perform periodicity. The case of f e /f s = 0.75 is an exception to the "lock-in" range, which could be due to the fact that it is in an intermediate state from the 2S mode to the 2P 0 mode or it could be that the dominant frequency is a fusion between the excitation frequency f e and the Stroulhal frequency f s .
In the case of f e /f s = 0.95, the vortex shedding mode is transformed into a 2P mode. Moreover, the secondary vortex of each vortex pair has lower strength than that of the primary vortex. This vortex shedding mode corresponds to a 2P overlap mode, 2P 0 , which is consistent with the vortex shedding mode at f e /f s ≈ 1 obtained by Morse and Willamson [3]. Additionally, the vortex shedding mode is within the range of the "lock-in" phenomenon as a stable state and the vortex shedding frequency of the oscillating cylinder approaches the natural shedding frequency, which can lead to the approximate symmetry of the wake vortices.
As the frequency ratio f e /f s increases to 1.05, the vortex shedding mode steps in an intermediate state from a 2P mode to a P + S mode. Two vortex pairs exist per cycle. The strength of one vortex pair is significantly stronger than that of the other, and the weak vortex pair dissipates faster when it moves downstream. The strong and weak vortex pairs are separated in the process of the downward and upward oscillation, respectively. More importantly, the asymmetry of the strength of the vortex shedding leads in the asymmetry of the temporal evolution of the lift coefficient per cycle. In other ways, an irregular "jumping" phenomenon occurs at the right side of the temporal evolution of the lift coefficient per cycle when the cylinder moves to near a negative extreme displacement. However, the strong vortex pair is separated in the process of the upward oscillation at f e /f s ≈ 1.2, which is different from that at f e /f s ≈ 1.05, although they have the same vortex shedding modes. Accordingly, an irregular "jumping" phenomenon occurs at the left side of the temporal evolution of the lift coefficient per cycle.
From f e /f s = 1.05 to f e /f s = 1.2, the discrepancy in time of separation of the strong vortex pair mainly derives from the difference in the oscillation frequency. Additionally, the irregular "jumping" occurs when the oscillating direction switches. When f e /f s = 1.35, its vortex shedding mode is in line with that at f e /f s = 1.2, but the lateral displacement of the vortex shedding is greater than that at f e /f s = 1.2. From f e /f s = 1.05 to f e /f s = 1.35, the temporal evolutions of the lift coefficient all show a periodic fluctuation, and locations of the irregular "jumping' on them mainly depend on the time of the separation of the vortices.
As above, Figure 11 shows the vorticities around and behind the oscillating cylinder at A/D = 1.5 and f e /f s = 0.5 to 1.35. Figure 12 shows the temporal evolutions of the lift and drag coefficients (left) and the PSD of the lift coefficient of the oscillating cylinder (right) in the corresponding parameters. As A/D further increases to 1.5, the lateral displacement of the vortex shedding is also increased.
Similarly, the vortex shedding mode also shows a similar "P + S" model at A/D = 1.5 and f e /f s = 0.5, and the lift and drag coefficients also show the "beat" characteristics. The temporal evolutions of the lift and drag coefficients also fluctuate drastically. The PSD of the lift coefficient at f e /f s = 0.5 shows that the vortex shedding frequency is affected not only by f e and f s , but also by other nondominant frequencies, as shown in Figure 12a. The vortex shedding mode shows a 2P 0 mode placed out of the "lock-in' range, because the vortex shedding in the wake is either a single vortex or a vortex pair (including a strong vortex and a weak vortex). The phenomenon is also caused by the joint action of the excitation frequency and the Stroulhal frequency.
When f e /f s = 0.75, the vortex shedding mode has transformed into a 2P 0 mode. The PSD of the lift coefficient at f e /f s = 0.75 also verifies that the vortex shedding frequency only rests with the excitation frequency f e . However, the vortex shedding has a larger lateral displacement in downstream as compared with the case at A/D = 1 and f e /f s =0.95, although the vortex shedding mode is the same. Additionally, this phenomenon also indicates that the transition of the vortex shedding mode is more likely to occur at high excitation amplitude.
vortex shedding also increases with the increase of the oscillation frequency ratio. As the excitation amplitude increase to A/D = 1.5, the vortex shedding mode exhibits a P + S mode in high oscillation frequency ratios. An unstable status exists in the mode transition from the 2P0 mode to the P + S mode, which is also considered as an intermediate state. Additionally, the case of a higher excitation amplitude and oscillating frequency ratio has a vortex shedding mode featuring a strong vortex pair and a single weak vortex.  In the case of f e /f s = 0.95 and f e /f s = 1.05, the vortex shedding modes are also in an intermediate state from the 2P 0 mode to the P + S mode. However, the vortex shedding exhibits a small lateral displacement in downstream, converging towards the wake centerline. The temporal evolution of the lift coefficient at f e /f s = 0.95 shows a similar irregular "jumping" in the process of the downward and upward oscillation, respectively, i.e., double "jumpings". The temporal evolution of the lift coefficient at f e /f s = 1.05 also exhibits the double "jumpings". Additionally, the strength of the vortex shedding remains the same in the process of the downward and upward oscillation, although the shedding vortex can be a single vortex or a vortex pair.
When f e /f s = 1.2 and f e /f s = 1.35, the vortex shedding modes keep a P + S mode. The single vortex is still weaker and is dissipated quickly in the downstream. However, the vortex pair enjoys stronger strength, and spreads upward. Therefore, the temporal evolution of the lift coefficient also shows an irregular "jumping" on its right side per cycle. Figure 13 shows a diagram of the vortex formation modes in the above computed range at Re = 4000. The results feature complex vortex formation transitions. A "lock-in" line at f e /f s = 0.75 determines the components of the vortex shedding frequency for a series of simulations of A/D = 0.5 to 1.5 and f e /f s = 0.5 to 1.35. The vortex shedding is dominated by the excitation frequency and the Stroulhal frequency when f e /f s < 0.75, located in the left of the "lock-in" line. However, the vortex shedding is only determined by a dominated vortex shedding frequency when f e /f s ≥ 0.75, which is within the range of the "lock-in" phenomenon. For a lower excitation amplitude, the vortex shedding mode shows a 2S mode. With an increase in the excitation amplitude, the vortex shedding mode turns into a 2P 0 or 2P mode, corresponding to lower oscillation frequency ratios and higher ones, respectively. The mentioned 2P mode characterizes a combined mode including a strong vortex pair and a weak one, which occur in higher oscillation frequency ratios. In addition, the strength of the vortex shedding also increases with the increase of the oscillation frequency ratio. As the excitation amplitude increase to A/D = 1.5, the vortex shedding mode exhibits a P + S mode in high oscillation frequency ratios. An unstable status exists in the mode transition from the 2P 0 mode to the P + S mode, which is also considered as an intermediate state. Additionally, the case of a higher excitation amplitude and oscillating frequency ratio has a vortex shedding mode featuring a strong vortex pair and a single weak vortex.   The symbol × denotes a case that is not within the "lock-in' range and the symbol ■ denotes a case that is within the "lock-in" range. Ps, strong vortex pair; Pw, weak vortex pair; Sw, weak single vortex. The unstable status means that the vortex shedding is a vortex pair and a single vortex without regularity. The red arrow represents that the "lock-in" range is located in the right of the "lock-in" line.

Conclusions
In this paper, a high-efficiency VLES method is adopted to simulate the VIV phenomenon of the oscillating cylinder at a subcritical Reynolds number Re = 4000 combined with the dynamic grid technique. A series of simulations of A/D = 0.5 to 1.5 and fe/fs = 0.5 to 1.35 are conducted, and the temporal evolutions of the lift and drag coefficients, vortex shedding frequency, and wake structures are discussed in detail. The main conclusions are as follows: (1) The VLES method combined with the dynamic mesh technology is adopted to simulate the flow past an oscillating cylinder in this paper, which is also verified to be able to effectively handle the dynamic boundary problems, such as the oscillating cylinder.
(2) The components of vortex shedding frequency of an oscillating cylinder are first determined by a "lock-in" bound, which is labelled as "lock-in" line at fe/fs = 0.75 in this paper. The vortex shedding is dominated by the excitation frequency and the Stroulhal frequency when fe/fs < 0.75, but the vortex shedding is only determined by a dominated vortex shedding frequency when fe/fs ≥ 0.75. The latter is within the range of the "lock-in" phenomenon. The temporal evolutions of the lift and drag coefficients show an obvious "beating" characteristic in the condition of the "lock-in" phenomenon. Even at higher excitation amplitude, multiple nondominant frequencies are also involved.
(3) The vortices located in the "lock-in" line maintain the same shedding mode and level of strength in the wake per cycle, and the temporal evolutions of lift and drag coefficients also perform periodic fluctuation. For a lower excitation frequency, the vortex shedding mode shows a 2S mode. With an increase in the excitation amplitude, the vortex shedding mode turns into a 2P0 or 2P mode, corresponding to lower oscillation frequency ratios and higher ones, respectively. As the excitation amplitude continues to increase, the vortex shedding mode exhibits a P + S mode in high oscillation frequency ratios.
(4) As the oscillating frequency ratio increases, the vortex shedding mode also changes. The 2P mode characterizes a combined mode including a strong vortex pair and a weak one in a higher oscillation frequency ratio and a medium excitation amplitude. The vortex shedding mode varies from a 2P0 or 2P mode to a P + S mode under the condition of a high frequency ratio and Figure 13. The vortex formation modes in the computed range at Re = 4000. The symbol × denotes a case that is not within the "lock-in' range and the symbol denotes a case that is within the "lock-in" range. P s , strong vortex pair; P w , weak vortex pair; S w , weak single vortex. The unstable status means that the vortex shedding is a vortex pair and a single vortex without regularity. The red arrow represents that the "lock-in" range is located in the right of the "lock-in" line.

Conclusions
In this paper, a high-efficiency VLES method is adopted to simulate the VIV phenomenon of the oscillating cylinder at a subcritical Reynolds number Re = 4000 combined with the dynamic grid technique. A series of simulations of A/D = 0.5 to 1.5 and f e /f s = 0.5 to 1.35 are conducted, and the temporal evolutions of the lift and drag coefficients, vortex shedding frequency, and wake structures are discussed in detail. The main conclusions are as follows: (1) The VLES method combined with the dynamic mesh technology is adopted to simulate the flow past an oscillating cylinder in this paper, which is also verified to be able to effectively handle the dynamic boundary problems, such as the oscillating cylinder. (2) The components of vortex shedding frequency of an oscillating cylinder are first determined by a "lock-in" bound, which is labelled as "lock-in" line at f e /f s = 0.75 in this paper. The vortex shedding is dominated by the excitation frequency and the Stroulhal frequency when f e /f s < 0.75, but the vortex shedding is only determined by a dominated vortex shedding frequency when f e /f s ≥ 0.75. The latter is within the range of the "lock-in" phenomenon. The temporal evolutions of the lift and drag coefficients show an obvious "beating" characteristic in the condition of the "lock-in" phenomenon. Even at higher excitation amplitude, multiple nondominant frequencies are also involved. (3) The vortices located in the "lock-in" line maintain the same shedding mode and level of strength in the wake per cycle, and the temporal evolutions of lift and drag coefficients also perform periodic fluctuation. For a lower excitation frequency, the vortex shedding mode shows a 2S mode. With an increase in the excitation amplitude, the vortex shedding mode turns into a 2P 0 or 2P mode, corresponding to lower oscillation frequency ratios and higher ones, respectively. As the excitation amplitude continues to increase, the vortex shedding mode exhibits a P + S mode in high oscillation frequency ratios. (4) As the oscillating frequency ratio increases, the vortex shedding mode also changes. The 2P mode characterizes a combined mode including a strong vortex pair and a weak one in a higher oscillation frequency ratio and a medium excitation amplitude. The vortex shedding mode varies from a 2P 0 or 2P mode to a P + S mode under the condition of a high frequency ratio and a medium or high excitation amplitude. The vortex shedding mode is an unstable status of the P + S mode for a medium frequency ratio and a high excitation amplitude. However, the case of a higher excitation amplitude and oscillation frequency ratio has a P + S mode featuring a strong vortex pair and a single weak vortex. (5) The vortex shedding is a lasting process under the condition of a low excitation amplitude, leading the lift and drag coefficients to fluctuate irregularly. For a vortex shedding mode exhibiting a strong vortex pair and a weak vortex pair or a weak single vortex, the temporal evolution of the lift coefficient of the oscillating cylinder shows an irregular "jumping" at the specific time per cycle corresponding to the shedding of the strong vortex pair. Moreover, the vortex shedding usually occurs when the direction of oscillation switches. The vortex shedding mode, and the frequency and time of the vortex shedding co-determine the temporal evolutions of the lift and drag coefficients.