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

Abstract

This work focuses on flow 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 flow characteristics. In this paper, a high-efficiency very large-eddy simulation method is adopted and verified in order to handle the oscillating boundary. A series of numerical simulations are conducted to investigate the transient flow 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 fluctuations of the lift and drag coefficients. 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 coefficient of the oscillating cylinder shows irregular ”jumping” at a specific 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 coefficient.

1. 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).

Williamson and Rishko studied the flow characteristics of a circular cylinder at various Reynolds numbers (Re), Re = 300 to 1000 [2] and obtained three main vortex shedding modes that existed in their experiments when meeting f_e/f_s ≈ 1, including 2S, 2P, and P + S modes. The 2S mode signifies two single vortices shed per oscillation cycle, which is also the natural Kármán mode of shedding. The 2P mode indicates two count-rotating vortex pairs formed per oscillation cycle, and the P + S asymmetric mode shows a vortex pair on one side and a single vortex on the other side. Moreover, the vortex shedding mode changes from 2P to 2S when f_e/f_s increased through unity. Morse and Willamson experimentally investigated the amplitude-response prediction of a forced oscillation with much higher resolution. In a high-amplitude regime, they found two vortex shedding modes overlap at Re = 4000 and 12,000 [3]. This vortex formation mode includes two vortex pairs per oscillation cycle. Importantly, the primary vortex in each vortex pair is stronger than the secondary, which also is recognized as a 2P overlap mode (2P₀) in the condition of f_e/f_s ≈ 1.

Carberry et al. researched the wake state by using controlled oscillation experiments based on amplitude-to-cylinder ratios (A/D) that increased from 0.25 to 0.6 at Re = 2300, 4410, and 9100 [4], where D is the cylinder diameter. The high- and low-frequency wake states showed the 2S and 2P modes, respectively, which also impacted the near wake structure and the lift. The high- and low-frequency wake states were transited through an intermediate wake state. Pastrana et al. conducted a high-fidelity large-eddy simulation (LES) study to investigate the two-degree-of-freedom (2DOF) VIV phenomenon of a low mass ratio circular cylinder at Re = 3900, 5300, and 11000 [5,6]. Their simulations presented that a 2T vortex formation mode, namely two triplets of vortices per cycle, was an inherent feature of the flow in the super-upper branch. Marble et al. experimentally studied the wake development of a stationary cylinder and the VIV of a pivoted cylinder with elliptical trajectories at a fixed Re = 3027. Their experiments showed that the wake topology of a 2DOF VIV with different amplitudes significantly deviated from the shedding map of an one-degree-of-freedom (1DOF) forced VIV originated from Morse and Williamson [7,8].

Blackburn and Henderson adopted two-dimensional (2D) direct numerical simulations (DNSs) to study the hydrodynamic forces acting on a circular cylinder at Re = 500 [9]. They found that the lift exhibited a variation when the wake vortex changed. However, the state of the wake mode was not reported in their research. Dong and Karniadakis et al. performed a three-dimensional (3D) DNS of a forced oscillating cylinder at Re = 10,000 [10]. Their simulation effectively predicted the lift and drag coefficients through a comparison with the experimental results, but the wake structure was also not illustrated in their study. Although the 2D simulation has some limitations, it has been extensively employed in high Reynolds number VIV studies [10]. Atluri et al. studied the effect of a modulated excitation frequency on a forced oscillating cylinder at Re = 500 to 8000 by using an extensive 2D LES and qualitatively validated their conclusion through a comparison with the experimental data [11]. To accurately simulated the VIV phenomenon, numerical simulation commonly was restricted to low Reynolds numbers. This was caused by the fact that a large number of computational resources and additional modeling works needed to be invested to reasonably describe the flow regime beyond the laminar flow [12,13,14].

Significant improvements have been made to address the difficulties in numerical simulations of the VIV phenomenon. Currently LESs are increasingly applied to calculate the separation flows at subcritical, critical, and supercritical Reynolds numbers. [5,15,16] The LES method can accurately solve the large-eddy structures, because the transient large-scale turbulent motions are explicitly calculated. However, the LES method is usually not computationally feasible, since it is sensitive to the grid resolution requirement near the wall when the global flow dynamics is significantly impacted by the attached boundary layers. Additionally, LES often requires more computing resources [17]. To reduce the computational cost, hybrid turbulence modeling methods have been extensively developed, such as a very LES (VLES) method. The VLES method has also been used in an unsteady flow to handle the transient dynamics of a separated flow [18].

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.

2. Materials and Methods

2.1. 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:

$${\tau }_{ij}^{sub}=Fr{\tau }_{ij}^{RANS}$$

(1)

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]:

$${\mu }_t=Fr\rho C_{\mu }{k}^2/\epsilon .$$

(2)

In this VLES method, the transport equations are expressed as:

$$\frac{D\rho k}{Dt}=P_k-\rho \epsilon +\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right],$$

(3)

$$\frac{D\rho \epsilon }{Dt}=\frac{\epsilon }{k}\left(C_{\epsilon 1}{P}_k-C_{\epsilon 2}\rho \epsilon \right)+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right].$$

(4)

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:

$$P_k=2{\mu }_t{S}_{ij}{S}_{ij},$$

(5)

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

(6)

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:

$$Fr=\min (1.0,{\left[\frac{1.0-exp(-\beta L_c/L_k)}{1.0-exp(-\beta L_i/L_k)}\right]}^n)$$

(7)

where β is the model constant and equals 2.0 × 10⁻³.

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:

$$L_c=C_x{(\Delta x\Delta y\Delta z)}^{\frac{1}{3}},$$

(8)

$$L_k=v^{\frac{3}{4}}/{\epsilon }^{\frac{1}{4}},$$

(9)

$$L_i=k^{\frac{3}{2}}/\epsilon .$$

(10)

where $\left(\Delta x\Delta y\Delta z\right)$^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.

2.2. 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_et). 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. Elements and nodes distribution of the grid convergence study.

	Nodes (X × Y)	Nodes (Cylinder)	Total Elements
Mesh-I	51 × 34	160	6626
Mesh-II	101 × 68	198	17,894
Mesh-III	201 × 134	240	60,418
Mesh-IV	301 × 201	300	131,426

.

Table 2. Strouhal (St) numbers calculated by the four grids as compared with the results from the literatures.

Mesh-I	Mesh-II	Mesh-III	Mesh-IV	Ref. [21]	Ref. [22]	Experiments [16,23]
0.21995	0.21995	0.22493	0.22495	0.21 (Re = 3900)	0.22 (Re = 3900)	0.2–0.22 (Re = 3900)

shows the comparison between the Strouhal number (St = f_sD/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.

According to the empirical formula provided by Norberg et al., the St number is 0.21 when Re = 3900 [21]. Wornom et al. employed the variational multiscale LES (VMS-LES) to simulate the 3D stationary cylinder at Re = 3900 and obtained a St number of 0.22 [22]. However, Parnaudeau et al. 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.

2.3. 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_esin(2πf_et), where A_e is the oscillating amplitude. The cylinder initially moves towards the negative x coordinate with a velocity U_X = −2πf_eA_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πfA, 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²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]

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_et) 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 Figure 2 and Figure 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.

2.4. Correctness of the VLES Method

Figure 5 exhibits the Fr distribution of the computational domain at A/D = 1 (ratio of the excitation amplitude to the cylinder diameter) and f_e/f_s = 0.75. According to the principle of the VLES method, the Fr value determines the computational model in the computation domain. When Fr approximates 1.0 near the wall, the computational model is recalled into the RANS model to effectively solve the flow near the boundary of the 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.

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.

3. 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. Oscillating amplitudes and frequencies of present numerical simulation.

A/D	fe/fs
0.5	0.5	0.75	0.95	1.05	1.2	1.35
1	0.5	0.75	0.95	1.05	1.2	1.35
1.5	0.5	0.75	0.95	1.05	1.2	1.35

. 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:

$${{\displaystyle \int }}_0^{nT}Cl\left(t\right)e^{i2\pi ft}dt=Cl\left(f\right),$$

(11)

$$\mathrm{PSD}=\underset{T\to \infty }{\lim }\frac{{\left|Cl\left(f\right)\right|}^2}{nT}.$$

(12)

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:

$$Cl=F_l/0.5\rho D{u}^2,$$

(13)

$$Cd=F_d/0.5\rho D{u}^2$$

(14)

and

$$F_l=-{{\displaystyle \int }}_0^{2\pi }\left(p\left(t\right)cos\theta -\frac{2}{R_e}{\omega }_z\left(t\right)sin\theta \right)d\theta ,$$

(15)

$$F_d=-{{\displaystyle \int }}_0^{2\pi }\left(p\left(t\right)sin\theta +\frac{2}{R_e}{\omega }_z\left(t\right)cos\theta \right)d\theta .$$

(16)

where $p\left(t\right)$ and ${\omega }_z\left(t\right)$ are obtained by a dimensionless transformation of the pressure and the axial vorticity on the surface of the oscillating cylinder, respectively.

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.

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₀ 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₀, 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₀ 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₀ 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.

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₀ 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₀ 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₀ 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.

4. 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₀ 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₀ 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.