3.1. Critical ${\textit {Re}}$ for the onset of vortex shedding

Before investigating the loss of stability of the Kármán wake to 3-D perturbations, it is useful to determine the underlying periodic base-flow states. Although this dependence of the transition on reduced velocity has been examined previous by Mittal & Singh (Reference Mittal and Singh2005) for the two-degree-of-freedom case, they only considered mass ratios $m^* \gtrsim 5$. For the current simulations, this range is extended to lower mass ratios where the effect on transition is amplified.

For a stationary circular cylinder, the wake undergoes a transition from a stable to periodic Kármán wake at ${\textit {Re}} \sim 47$ (e.g. Dusek, Gal & Fraunie Reference Dusek, Gal and Fraunie1994; Williamson Reference Williamson1996; Le Gal, Nadim & Thompson Reference Le Gal, Nadim and Thompson2001). This instability leads to unsteady lift and drag forces that cause an elastically mounted cylinder to undergo VIV. The onset of this periodic shedding occurs at different Reynolds numbers for a cylinder undergoing VIV. The value of this critical ${\textit {Re}}$ depends on the reduced velocity.

Based on the stability analysis of § 2.5, figure 7 provides growth rate contours for mass ratios $m^*=1$, $2.546$ and $10$. The zero growth rate curves are marked, which indicate the transition Reynolds number as a function of reduced velocity.

Figure 7. Growth rate contour maps for the steady to unsteady transition of an elastically mounted cylinder showing the effect of mass ratio on transition as a function of reduced velocity. The rows correspond to (a) $m^* = 1$, (b) $m^*=2.546$ and (c) $m^*=10$. The zero growth contour is marked indicating the variation of the transition Reynolds number as $U^*$ is varied.

As $U_r$ is increased from very small values, the critical Reynolds number, $Re_{co}$, initially rises from a value consistent with the stationary cylinder transition Reynolds number of $Re_{co} \simeq 46$ to reach a peak, prior to rapidly dropping after a different instability mode becomes dominant (corresponding to the system frequency rather than the unsaturated shedding frequency). This shift corresponds to the kinks in the contour lines of figure 7. For the lightest mass ratio considered, $m^*=1$, the critical Reynolds number is increased to approximately $Re_{co}=51$ at $U_r=4.6$. For larger mass ratios, the effect is less; for $m^*=10$, the increase is only to $Re_{co} \simeq 47$ at $U_r=5.7$ before again dropping rapidly. In all cases, beyond the initial peak, the critical Reynolds number reaches a value of slightly less than 20, noting that the minimum occurs at different reduced velocities for the different mass ratios. This result is consistent with the findings of Mittal & Singh (Reference Mittal and Singh2005), who also found a minimum critical Reynolds number of ${\sim }20$.

3.2. Variation of the mode growth close to 3-D transition

In this section, for the system undergoing VIV, loss of stability of the Kármán wake to 3-D perturbations is analysed for a mass ratio $m^* = 2.546$, as a function of Reynolds number. The Reynolds number range chosen covers the range where analogues of both modes A and B become unstable. Figure 8 shows the magnitude of the predicted Floquet multiplier, $|\mu |$, vs wavelength, $\lambda /D$, for undamped free vibration. The natural frequency is fixed at $f_n^* = 0.2$, close to shedding frequency of the fixed cylinder; this corresponds to $U_r = 5.0$. Figure 8 shows that for ${\textit {Re}}= 300$ and 320, the flow is unstable to 3-D perturbations, whereas at ${\textit {Re}} = 280$ and $260$ the 2-D wake flow remains stable ($|\mu | < 1$), with no Floquet modes showing positive growth. The amplitude of oscillation increases slightly with decreasing ${\textit {Re}}$, with $A^* = 0.49, 50$, $0.50$ and $0.51$, for ${\textit {Re}}= 320$, $300$, 280 and 260, respectively. The figure shows that for ${\textit {Re}} \geq 280$, two clear peaks can be observed although both remain stable at $Re=280$, whereas for ${\textit {Re}} = 260$ only the mode B peak has become visible. The modes corresponding to the two branches are analogous to modes A and B for the stationary cylinder. We observe that mode B becomes critical first then mode A; a similar result was found for $A^* > 0.3$ under forced cylinder oscillation by Leontini et al. (Reference Leontini, Thompson and Hourigan2007).

Figure 8. Floquet multiplier vs wavelength at ${\textit {Re}} = 320$, $300$, $280$ and $260$; $\zeta = 0$ (undamped system), $f_n^* = 0.20$ and $m^* = 2.546$.

In the case of ${\textit {Re}}= 300$, only mode B is unstable, with the maximum growth rate occurring at $\lambda \approx 0.9 D$. The multiplier branch corresponding to Mode QP, which has a complex Floquet multiplier, remains stable but is dominant for $1.4D < \lambda < 1.9 D$. At ${\textit {Re}} = 300$, mode A is also stable with a maximum $|\mu | \approx 0.9$ occurring at $\lambda \approx 2.5D$. In fact, for $\lambda > 3.7 D$, the dominant mode has a complex Floquet multiplier corresponding to a QP mode. Since it is observed at a longer wavelength than mode A, to make a clear distinction from the previously observed QP mode, it is designated as mode QPL, named for long-wavelength QP. For ${\textit {Re}} = 280$, the flow is near to unstable to mode B, again with a wavelength centred around $\lambda \approx 0.9D$, whereas the branch corresponding to mode QP remains stable and exists for $\lambda > 1.4 D$ although it is not shown explicitly here. Thus, modes A and B become critical at lower Reynolds numbers than those of modes QP and QPL. This latter behaviour is similar as that reported by Blackburn, Marques & Lopez (Reference Blackburn, Marques and Lopez2005) for a stationary cylinder, who estimated the critical ${\textit {Re}}$ for mode QP to be about 377 with $\lambda _{cr} \approx 1.8 D$. Mode QP has no local maximum over this Reynolds number range, unlike modes A and B; however, we have not examined the stability for ${\textit {Re}} > 300$.

Figure 9 shows the spanwise and streamwise perturbation vorticity fields at ${\textit {Re}} = 300$ in order to indicate the similarities and differences in the structures of modes A and B. The solid lines in the figure outline the base-flow vorticity with two single vortices shed per oscillation cycle (referred to as a $2S$ wake mode). This mode occurs at low amplitudes of oscillation, in particular for $A^* \lesssim 0.5$ at ${\textit {Re}} =300$ (Leontini et al. Reference Leontini, Stewart, Thompson and Hourigan2006a). Mode B in figure 9(a,b) is the first occurring unstable mode showing a maximum Floquet multiplier at $\lambda = 0.9D$. Mode B obeys the spatiotemporal symmetry of the base flow about the wake centre line, with spanwise perturbation vorticity repeated (on reflection through the centreline) every half base period. The successive spanwise vortex structures from each side of the wake have the same sign. The still stable mode A correspondingly shows a maximum $|\mu |$ at $\lambda = 2.5D$. This mode is shown in figure 9(c,d). The primary vortex at the rear of the cylinder has both positive and negative regions of perturbation vorticity; further downstream, the perturbation becomes concentrated in the primary vortex cores. Mode A has spatiotemporal symmetry opposite to that of mode B. The perturbation field structures for modes A and B are found to be similar to those of modes A and B for a fixed circular cylinder (Barkley & Henderson Reference Barkley and Henderson1996), although the oscillation certainly has a strong effect on the underlying base flow. In particular, the near-wake base-flow vortex structures are quite different from those for a stationary cylinder in that they are stretched and elongated as the cylinder moves up and down rather than more elliptical in character for the stationary cylinder. Since mode A has been at least partially attributed to elliptic instability (Leweke & Williamson Reference Leweke and Williamson1998; Thompson et al. Reference Thompson, Leweke and Williamson2001), this may help to explain why mode A remains stable until higher Reynolds numbers under oscillation. Since mode B is likely to depend more on the shear of the braid regions between Kármán wake vortices (Leweke & Williamson Reference Leweke and Williamson1998), which of course persist in the oscillatory case, this is consistent with the smaller effect on the onset of mode B. However, this question is further considered in the following through an examination of the perturbation kinetic energy generation.

Figure 9. Spanwise perturbation vorticity depicting (a,b) mode B and (c,d) mode A at ${\textit {Re}} = 300$, $\zeta = 0$, $f_n^* = 0.20$ and $m^* = 2.546$. Red (blue) represents positive (negative) spanwise vorticity. The solid lines outline base-flow vorticity.

The effect of mass ratio on 3-D stability under VIV at ${\textit {Re}}= 300$ is examined for mass ratios $m^* = 1$, 2.546 and 10 in figure 10. The reduced velocity is adjusted to result in $f^* = 0.20$ in each case. The response amplitude is close to 0.49 for each mass ratio. Indeed, the dominant Floquet multiplier curves are virtually identical as the mass ratio is varied by a factor of 10, again showing a collapse against the actual system response frequency (rather than $U_r$ or $U^*$), and consistent with figure 6.

Figure 10. Floquet multiplier vs wavelength for (a) $m^* = 10$, 2.546 and 1, at $Re=300$ for $f^* = 0.2$.

3.3. Relationship of the 3-D instability modes to stationary cylinder modes A, B and QP

Figure 11 plots the spatial distributions of the perturbation kinetic energy ($k'$) and generation rate term ($\dot {K}'$) given by (2.10) for the different instability modes. This latter term ($\dot {K}'$) provides the rate of transfer of perturbation kinetic energy to or from the base flow, thus providing insight into which parts of the wake contribute to the development of the instability modes. This analysis was undertaken for modes A, B and QP for the stationary circular cylinder and the analogous A, B and QP instability modes for a cylinder undergoing VIV.

Figure 11. Perturbation kinetic energy, $k'$, (a,c,e,g,i,k), and kinetic energy generation rate, $\dot {K}'$, (b,d, f,h, j,l), from (2.10), which indicates the generation of perturbation kinetic energy from the base flow: (a–d) mode A comparison for the stationary and VIV cases; (e–h) same for mode B, bottom two rows; (i–l) same for mode QP. See the text for further details.

Each pair of rows of this figure shows side-by-side comparisons of the dominant stationary-cylinder and VIV oscillating-cylinder instability modes close to onset. The left column shows perturbation kinetic energy and the right column the perturbation energy generation rate. The top, middle and bottom pair of rows show comparisons for modes A, B and QP, respectively. Clearly, for mode A of the stationary cylinder, energy input into the mode occurs primarily in the wake vortex cores, consistent with the interpretation of elliptic instability assisting the maintenance of the instability (Leweke & Williamson Reference Leweke and Williamson1998; Thompson et al. Reference Thompson, Leweke and Williamson2001). On the other hand, for mode B perturbation energy transport from the base flow is primarily in the sheared region between the near-wake vortices and extending into the braids between the row of vortices further downstream. This is consistent with the interpretation that the instability is supported mostly by an instability of the braids (Leweke & Williamson Reference Leweke and Williamson1998). Mode QP appears to be supported through a more equal combination of elliptic and hyperbolic instabilities (Sheard et al. Reference Sheard, Thompson and Hourigan2003).

In terms of the VIV instability modes, mode A shows distinct differences from the $k'$ and $\dot {K}'$ distributions of the stationary cylinder modes. Under oscillation, there appears to be little generation associated with the vortex cores and more generation in the braids. In addition, further downstream, the instability mode appears supported by energy transfer in the vortex cores, while for the stationary cylinder mode A the generation remains small. The study by Leontini et al. (Reference Leontini, Thompson and Hourigan2007) indicates that the mode A seen for a stationary cylinder wake is continuously transformed into the mode A seen under high-amplitude oscillations as the forced oscillation amplitude is steadily increased. Thus, an interpretation of the perturbation kinetic energy analysis is that the stretching of vortex structures in the near-wake caused by the body oscillation suppresses the presence of near-wake elliptic instability, which causes the mode A instability to occur at considerably higher Reynolds numbers ($Re_c\simeq 305$ at $U_r=5$) than for the stationary case ($Re_c \simeq 190$).

For mode B similar spatial distributions of $k'$ and $\dot {K}'$ are seen for both the stationary and oscillating cylinder modes. Given this, it is not surprising that the instability mode onset is similar in both cases: $Re_c = 259$ for the stationary cylinder and $Re_c=277$ at $U_r=5$ for the oscillating cylinder. Thus, mode B is relatively minimally affected by large amplitude oscillations and the non-trivial change to the wake structure, noting a significantly increased wake width.

As shown later in § 3.6, a QP mode is the first occurring VIV instability mode at $U_r=7$. The comparison of the stationary and VIV QP modes in figure 11 indicates that the VIV QP mode more closely resembles mode B than the stationary cylinder QP mode. Hence, it is unsurprising that the critical Reynolds number for these modes is so different ($Re_c\simeq 377$ for the stationary cylinder mode and $Re_c\simeq 220$ for the VIV mode).

Further discussion on how these modes change with reduced velocity is given in § 3.6.

3.4. Response curves for a circular cylinder under VIV at ${\textit {Re}} = 150$, 200 and 250

The VIV response of a circular cylinder for an undamped system is shown in figure 12 at $m^* = 2.546$ for ${\textit {Re}} = 150$, 200 and 250. The following variables, the maximum oscillation amplitude ($A^*$), peak lift coefficient ($C_{L,max}$), normalised frequency ( $f^*$), total phase difference ($\phi _{tot}$) and vortex phase difference ($\phi _{vor}$), are plotted as a function of reduced velocity ($U_r)$ in figure 12(a–d), respectively. Normalised frequency is defined as $f^* = f_y/f_n$, where $f_y$ is the frequency corresponding to the highest power in the power spectral density of displacement, and $f_n$ is the natural frequency of the structural system. The total phase $(\phi _{tot})$ and the vortex phase ($\phi _{vor}$) are the phase differences between fluid force and displacement signals, obtained using a Hilbert transform (Khalak & Williamson Reference Khalak and Williamson1999). In particular, $\phi _{tot}$ and $\phi _{vor}$ are the mean phase difference between lift force and cylinder displacement, and the vortex force (total $-$ added mass force) and cylinder displacement, respectively.

Figure 12. VIV responses of the cylinder at $m^* = 2.546$ for ${\textit {Re}} = 150$, 200 and 250: (a) maximum oscillation amplitude, $A^*$; (b) peak lift coefficient, $C_{L,max}$; (c) normalised frequency, $f^*\ (= f/f_n)$ (green dotted line represents the vortex shedding frequency for a stationary cylinder); (d) mean phase difference between lift force and cylinder displacement, $\phi _{tot}$; (e) mean phase difference between vortex force and cylinder displacement, $\phi _{vor}$.

We observe a typical three-branch response, the initial (red colour), upper (blue colour) and lower (black colour) branches, followed by a desynchronisation (magenta colour) regime. These branches are not necessarily one-to-one related to those observed for high-Reynolds number VIV at a low mass-damping ratio. Based on the observation by Govardhan & Williamson (Reference Govardhan and Williamson2000), Leontini, Thompson & Hourigan (Reference Leontini, Thompson and Hourigan2006b), Soti et al. (Reference Soti, Zhao, Thompson, Sheridan and Bhardwaj2018), Mishra et al. (Reference Mishra, Soti, Bhardwaj, Kulkarni and Thompson2020) and Sharma, Garg & Bhardwaj (Reference Sharma, Garg and Bhardwaj2022), the criteria to classify different VIV branches are as follows. The initial branch is demarcated with monotonically increasing $A^*$ and $C_{L,max}$, multiple frequencies in $f_y$, the vibration frequency overlaps the vortex shedding frequency of a stationary cylinder ( $f^* \approx f_v/f_s$), $\phi _{tot} \approx 0$ and $\phi _{vor} \approx 0$. The range of the upper branch is identified with large $A^*$, constant $f^*$ slightly lower than 1 ( $f^*<1$), large $C_{L, max}$, $\phi _{tot} \approx 0$ and transition of $\phi _{vor}$ from 0 to ${\rm \pi}$. Beyond this, conditions for the existence of the lower branch is large $A^*$, $f^* =1$, monotonic decrease in $C_{L, max}$, transition of $\phi _{vor}$ from 0 to ${\rm \pi}$ and $\phi _{vor} \approx {\rm \pi}$. Finally, desynchronisation is characterised by negligible $A^*$, $f^*$ parallel to $f_v/f_s$, $\phi _{tot} \approx {\rm \pi}$ and $\phi _{vor} \approx {\rm \pi}$.

As evident from figure 12(a) $A^*$ is a continuous function of $U_r$, except there are sudden jumps at $U_r = 3.3$, 3.1 and 3.0 for ${\textit {Re}} = 150$, 200 and 250, respectively, thus allowing the identification of the transition between the initial and the upper branches. Although displacement amplitude is continuous, figure 12(b) shows a sudden jump in the lift-force amplitude, thus demarcating the upper and lower branches at $U_r = 4.0$ (${\textit {Re}} = 150$), 3.9 (${\textit {Re}} = 200$) and 3.8 (${\textit {Re}} = 250$). Figure 12(c–e) shows the branch-to-branch jumps in the oscillation frequency and phase difference. We further quantify the branch ranges for ${\textit {Re}} = 150$, 200 and 250 in table 1, together with the criteria used to distinguish the branches.

Table 1. Range of $U_r$ and criteria distinguishing the different VIV branches for ${\textit {Re}} = 150$, 200 and 250.

Figure 12(a) shows that the peak oscillation amplitude at all three Reynolds numbers is ${\simeq }0.58D$, and the overall response curves are similar for ${\textit {Re}} = 150$, 200 and 250. These findings agree with Govardhan & Williamson (Reference Govardhan and Williamson2006), who showed that peak oscillation amplitude is almost constant for low-Reynolds-number flows. However, Williamson & Govardhan (Reference Williamson and Govardhan2004) has noted that such low-${\textit {Re}}$ simulations, in general, predict peak oscillation amplitudes around $0.6D$, whereas experiments for higher ${\textit {Re}}$ produce amplitudes up to 1$D$. Figure 13 shows wake vorticity contours at a reduced velocity, $U_r = 4.0$, corresponding to the maximum oscillation amplitude. At ${\textit {Re}} = 150$, the wake vortices assemble into a double-row, $2S$ (two single vortices per shedding period), configuration consistent with the high amplitude oscillation (figure 13a). At ${\textit {Re}} = 200$ and 250 (figure 13b,c), the newly shed vortices show signs of a $2P$ (two vortex pairs per shedding period) structure close to the rear of the cylinder, although the second vortex of each pair is very weak. However, we certainly do not claim that the wakes have a $2P$ structure for any of these Reynolds numbers, since the experimentally observed $2P$ structures show approximately similar strength opposite-signed vortices in each pair.

Figure 13. Vorticity contours at $U_r = 4.0$ for Reynolds number (a) ${\textit {Re}} = 150$, (b) ${\textit {Re}} = 200$ and (c) ${\textit {Re}} = 250$. The red (blue) colours represents clockwise (anticlockwise) spanwise vorticity.

On this point, Blackburn, Govardhan & Williamson (Reference Blackburn, Govardhan and Williamson2001) compared low Reynolds number experimental results ($Re \sim O(1000)$) with 2-D and 3-D numerical predictions. The 3-D simulations matched the experiments well. Furthermore, demarcation of the upper and the lower branches based on amplitude was not observed for the 2-D simulations, noting that they were clearly evident by the gaps between branches in the 3-D simulation. Therefore, and as has been recognised for a long time, 2-D simulations are inadequate for accurately modelling higher-${\textit {Re}}$ VIV, when the flow is 3-D and turbulent.

3.5. Effect of reduced velocity on different modes

In § 3.2, stability analysis for VIV is presented for an oscillation frequency of $f_n^* = 0.2$ at ${\textit {Re}} = 300$, corresponding to the strongly resonant case. In the current section, results for fixed ${\textit {Re}} = 300$ for a standard spring–damper VIV system are presented as a function of reduced velocity. Figure 14 shows the maximum amplitude of vibration vs reduced velocity at ${\textit {Re}} = 150$ and for $m^*=2.546$ depicting the periodic and QP nature of the oscillations. The situation is similar at other Reynolds numbers. Since FSA is only applicable to a periodic base flow, only the reduced velocity range $U_r \in [4.0, 8.0]$ is considered, noting periodic flow together with a strong resonant response.

Figure 14. Variation of amplitude ($A^*$) vs reduced velocity ($U_r$) for undamped VIV in the range $U_r \in [2.0, 8.0 ]$, depicting the periodic and QP range.

Figure 15 presents $|\mu |$ as a function of spanwise wavelength for $U_r \in [4.0- 8.0]$. The $U_r = 4.0$ case has a complex Floquet multiplier corresponding to a stable $2S$ wake mode, as shown in figure 16. At $U_r = 4.0$ the oscillation amplitude $A^* = 0.56$, corresponds to the highest amplitude during VIV. Interestingly, Leontini et al. (Reference Leontini, Thompson and Hourigan2007) observed a $P+S$ (base-flow) wake for forced oscillations at $A^* = 0.55$ and $f_s^* = 0.2$. That case corresponds to a forcing frequency of $St = 0.2$ or $U_r \simeq 5.0$. They observed the transition from a $2S$ to $P+S$ wake configuration at ${\textit {Re}} = 280$. Thus, the system frequency, which controls the shedding frequency over the resonant range, does have a non-negligible effect on the base-flow structure, with increasing the forcing/shedding frequency leading to a higher-Reynolds-number transition from the $2S$ to the $P+S$ state.

Figure 15. Floquet multiplier vs wavelength, at ${\textit {Re}} = 300, m^*=2.546$, $\zeta = 0$ for: (a) $U_r = 4.0$, (b) $U_r = 4.5$, (c) $U_r = 5.0$, (d) $U_r = 5.5$, (e) $U_r = 6.0$, ( f) $U_r = 6.5$, (g) $U_r = 7.0$, (h) $U_r = 7.5$ and (i) $U_r = 8.0$.

Figure 16. Base-flow vorticity contours at ${\textit {Re}} = 300$ for reduced velocity $U_r = 4.0$. The red (blue) represents clockwise (anticlockwise) spanwise vorticity.

For $U_r \geq 4.5$, modes A, B, QP, QPL (a long-wavelength QP mode) and SH are observed. Of interest, $U_r = 4.5$ shows that a subharmonic mode (SH) is the first to become unstable. This occurs for a wavelength range typically dominated by Mode A. Modes B, QP and QPL are not unstable, whereas mode SH is unstable and grows fastest at $\lambda \approx 1.4D$. Therefore, at $U_r = 4.5$, mode SH is responsible for 3-D transition. It is important to note that the Floquet multiplier for the SH mode is always real and negative consistent with its subharmonic nature. The spanwise perturbation vorticity field is shown in figure 17. The subharmonic mode nature is highlighted by comparing images at each base-flow time period ($T$), verifying that the perturbation field is repeated over two periods ($2T$) and the structure of the mode is the same at the end of a single period except for a sign change (Sheard et al. Reference Sheard, Thompson and Hourigan2005).

Figure 17. Spanwise perturbation vorticity indicating mode SH for time (a) $T$, (b) $2T$, (c) $3T$ and (d) $4T$ at $\lambda = 1.4 D$, $U_r = 4.5$, $Re = 300$ and $m^* = 2.546$. Here, $T$ is the base-flow field period.

For the $\lambda < 5 D$, mode QPL is only observed at reduced velocity $U_r = 4.5$ and 5.0, whereas mode QPL is absent for other reduced velocities $U_r \in [5.5, 8.0]$. It is also observed that as $U_r$ increases, the stability of the modes are affected. For $U_r = 4.0$ the wake is stable, whereas for $U_r = 4.5$ mode SH is just critical and modes B, QP and QPL can be observed. For $U_r=5.0$, mode A is observed for the first time but remains stable, however mode B is unstable. Furthermore, for $U_r \geq 5.0$, mode B is unstable ($|\mu | > 1$) and $|\mu |$ for this mode increases up to $U_r = 6.5$. The maximum $|\mu | \approx 1.7$ occurs at $\lambda \approx 0.9D$ for $U_r = 5.0$, increasing to a maximum $|\mu | \approx 4.0$ at $\lambda \approx 1.0D$ for $U_r = 6.5$. The growth rate for mode B further decreases at $U_r= 7.0$, with a maximum of $|\mu | \approx 3.5$ at $\lambda \approx 1.0D$, which remains constant on further increasing $U_r$. The maximum $|\mu | \approx 3.5$ occurs at $\lambda \approx 1.0D$ for $U_r = 7.0$, 7.5 and $8.0$.

Similar behaviour can be observed for mode A. At $U_r = 5.5$, the maximum $|\mu | \approx 1.03$ occurs at $\lambda \approx 2.7D$, indicating that the mode is just supercritical. The growth rate increases on further increasing $U_r$ up to $U_r = 7.5$ with a maximum $|\mu | \approx 1.8$ at $\lambda \approx 3.6D$. Beyond $U_r \geq 7.5$, the growth rate remains constant. Mode QP is seen to be stable for $U_r \in [4.0, 6.0]$ and unstable for $U_r \geq 6.0$. The maximum $|\mu |$ increases up to $U_r = 7.0$, with a highest growth rate of $|\mu | \approx 1.9$ at $\lambda \approx 1.9 D$. On increasing $U_r$ beyond 7.5, the maximum $|\mu |$ slightly decreases with a maximum of $|\mu | \approx 1.5$ and $|\mu | \approx 1.3$ at $U_r = 7.5$ and 8.0, respectively, centred around $\lambda = 1.9 D$. The flow is observed to be more stable at $U_r = 4.0$, 4.5 and 5.0, the highly resonant region. The synchronisation between cylinder vibration and vortex shedding stabilises the flow and hence it is more stable over this range.

3.6. Critical Reynolds number and wavelength variation with reduced velocity

This section presents the maximum ${\textit {Re}}$ for which the wake remains 2-D for the undamped VIV (standard spring–damper) system as a function of reduced velocity, $U_r$. We have focused our analysis on the lower branch of VIV, where the base flow exhibits periodic wake oscillations. In the desynchronisation branch, the amplitude is relatively small and periodic with the oscillation frequency close to but lower than the natural shedding frequency. On the other hand, for the initial and (pseudo-)upper branch, the oscillation is not periodic, so FSA, which requires a purely periodic signal, cannot be applied. Thus, we have limited our consideration to the resonant reduced-velocity range, $U_r \in [4.5, 8]$, where the base-flow amplitude variation is periodic, noting that at the upper end of this range corresponds to the desynchonisation branch.

Figure 18 shows the critical Reynolds number, ${\textit {Re}}_{cr}$, and corresponding critical wavelength, $\lambda _{cr}$, as a function of reduced velocity, $U_r$, for an undamped VIV system. Here, ${\textit {Re}}_{cr}$ is the value of the Reynolds number at which a given mode first becomes unstable, i.e. $\mu$ first exceeds $|\mu | = 1$. In addition, $\lambda _{cr}$ is the wavelength at which ${\textit {Re}}_{cr}$ occurs. To find an accurate ${\textit {Re}}_{cr}$ at each $U_r$, simulations were run for each ${\textit {Re}}$ for a series of $\lambda$ in steps of 0.1. For each ${\textit {Re}}$, $\mu$ locally varies quadratically about each local maximum (as shown in previous $\mu$ vs $\lambda$ plots). We determine the coordinate ($|\mu |, \lambda$) for the local maximum. The value of ${\textit {Re}}$ is then incremented in steps of five, noting down the maximum $|\mu |$ and the corresponding value of $\lambda$ for each ${\textit {Re}}$. This process is repeated until resolving values that give a maximum $|\mu | = 1$. The values of ${\textit {Re}}_{cr}$ are then found by interpolating ${\textit {Re}}$ and corresponding maximum $|\mu |$ values using quadratic interpolation. The same process is repeated for each reduced velocity.

Figure 18. (a) Critical Reynolds number, ${\textit {Re}}_{cr}$, and (b) critical wavelength, $\lambda _{cr}$, against reduced velocity for the undamped system at $U_r \in [4.5, 8.0]$, and $m^* = 2.546$.

Figure 18 shows the transition boundary for modes A, B, SH and QP, corresponding to reduced velocity $[4.5, 8.0]$. The critical Reynolds number, ${\textit {Re}}_{cr}$, decreases with reduced velocity up to $U_r=7.25$. At $U_r = 4.5$, mode SH is the only mode that becomes critical, which occurs at ${\textit {Re}} \approx 299$. For $U_r \in [5.0, 6.0]$, modes A and B become critical, with mode B being the leading 3-D mode. At $U_r = 6.25$, mode QP is critical first, with the order of inception of modes being mode QP, mode B and mode A. It is observed that for $U_r = 7.0$ and 7.5, there are two distinct QP modes (referred to as QP and QP$_0$) together with mode A that become critical, with mode QP$_0$ becoming critical first. Interestingly, mode QP$_0$ does not become critical at $U_r=6.25$. The perturbation kinetic energy analysis (see figure 11 and the discussion in § 3.3) suggests that the physical triggering mechanism for this mode has more in common with mode B rather than mode QP. It also occurs for the same wavelength band as for mode B, however, the transition Reynolds number is much lower than observed for mode B. Also for this $U_r$ range, at $U_r=7.25$ and $Re>270$, the wake develops asymmetry about the centreline and the wake period doubles. This part of parameter space is where the transition to mode B would be expected to occur. Transition to the desynchronisation branch occurs between $U_r = 7.5$ and 8.0, resulting in the relatively sudden shifting of shedding and cylinder oscillation to the be close to the natural shedding frequency, and a significant decrease in oscillation amplitude. At $U_r = 8.0$, mode A precedes the occurrence of modes B and QP, with modes B and QP transitioning at approximately the same ${\textit {Re}}$. Interestingly, the transition Reynolds numbers and wavelengths of modes A and QP at $U_r=8.0$ are consistent with the trends at lower $U_r$ values, despite the significant change in oscillation frequency and amplitude. Mode B, on the other hand, shows a significantly lower transition Reynolds number for this case.

The transition boundary, corresponding wavelength and amplitude for modes A, B, QP and SH, and for $U_r \in [4.5, 8.0]$, are summarised in table 2. These findings are in accordance with Leontini et al. (Reference Leontini, Thompson and Hourigan2007) for the forced oscillation case: for $A^* < 0.3$ mode inception order (with increasing ${\textit {Re}}$) is A, B and QP; whereas for $0.3< A^*<0.55$, the order of mode inception is B then A; and for $A^*>0.55$ mode SH is the only unstable mode. Table 2 indicates that for $U_r \in [4.5, 7.5]$, $A^*$ decreases approximately from 0.52 to 0.42 with $U_r$, whereas for $U_r = 8.0$, $A^* \approx 0.11$. Mode SH occurs at $A^* \approx 0.52$ for $U_r =4.5$, slightly varying from the prediction of Leontini et al. (Reference Leontini, Thompson and Hourigan2007), who observed mode SH at $A^* >0.55$. For $U_r \in [5.0, 6.0]$ with $0.42 \leq A^* \leq 0.50$ results are in agreement with Leontini et al. (Reference Leontini, Thompson and Hourigan2007). It is also observed that for $U_r = 6.5$ and 7.5, $0.40 \leq A^* \leq 0.44$, mode QP first becomes critical, which is different from Leontini et al. (Reference Leontini, Thompson and Hourigan2007); however, of course, the system frequency is then significantly different from the forcing frequency of $St = 0.2$ used in that study. On the other hand, at $U_r = 8.0$ with $A^* \simeq 0.11$ and the shedding frequency is close to 0.2, modes A and B undergo transition at close to the values seen for the forced oscillation case: $Re_{cr}=243$ at $St=0.18$ (VIV) compared with $Re_{cr}=255$ at $St=0.20$ (forced) for mode B; and $Re_{cr}=212$ at $St=0.18$ (VIV) compared with $Re_{cr}=230$ at $St=0.20$ (forced) for mode A. Overall, there appear to be non-negligible effects of the oscillation frequency and amplitude on transition.

Table 2. Transition boundary characteristics: corresponding critical Reynolds numbers, wavelengths and amplitudes for modes A, B, SH and QP. The post-simulation oscillation frequency is shown in the last column corresponding to the Reynolds number of the first transition.

The envelope of the transition lines in figure 18(a) shows the variation of critical Reynolds number with reduced velocity for which the flows remains 2-D. The critical Reynolds number decreases with reduced velocity over the entire range of $U_r$ considered. Specifically, ${\textit {Re}}_{cr} \approx 299$ at $U_r = 4.5$, decreasing monotonically to ${\textit {Re}}_{cr} \approx 208$ at $U_r = 7.5$.

Figure 18(b) shows the variation of $\lambda _{cr}$ with reduced velocity. In terms of the first modes to become unstable, the figure shows that $\lambda _{cr}\approx 1.4 D$ for mode SH at $U_r = 4.5$. For $U_r = 5.0$ to 6.5, $\lambda _{cr}\approx 0.9D$ for mode B, however, for mode A there is an increase in $\lambda _{cr}$ from $2.4D$ to $4.3D$. On further increasing $U_r$, $\lambda _{cr}$ for mode A decreases towards $3.9D$ at $U_r = 8.0$. Mode QP, which appears for $U_r \geq 6.5$, decreases from $\lambda _{cr}\approx 2.0D$ to $\lambda _{cr}\approx 1.0D$ over the interval $6.5 \leq U_r \leq 7.5$. It then increases to $\lambda _{cr}\approx 2.0D$ at $U_r = 8.0$. The critical values for $U_r = 8.0$, where the oscillation amplitude is only 0.11, are close to those for the non-vibrating case (Barkley & Henderson Reference Barkley and Henderson1996) as should be expected.

Overall, the synchronisation between cylinder vibration and vortex shedding has a significant role in maintaining the two-dimensionality of the flow. Clearly, high-amplitude oscillations delay onset of 3-D transition but there is also a distinct dependence on reduced velocity. In summary, for the standard elastically mounted VIV system, the maximum ${\textit {Re}}$ for which the wake remains 2-D over the resonant and periodic reduced-velocity range $U_r \in [4, 8]$ is 208 and occurs at $U_r = 7.5$, corresponding to a critical wavelength of $\lambda _{cr} \approx 1.2D$. The maximum delay in the transition from 2-D to 3-D flow occurs for $U_r = 4.5$, with the wake remaining 2-D up to ${\textit {Re}}_{cr} \approx 299$.

3.7. Influence of mass ratio on the amplitude response and transition

Figure 19 shows the influence of mass ratio on the amplitude response under VIV as a function of Reynolds number. The subfigures in the first column correspond to (a) $Re=150$ (where the wake is certainly 2-D), (b) $Re=200$ (with the wake on the verge of transition) and (c) $Re=300$ (where the wake is mostly 3-D). The response curves have a similar shape but with the onset reduced velocity of resonance shifted to the right as the mass ratio is increased. In addition, desynchronisation occurs at smaller $U_r$ values for higher mass ratios. The right-hand column provides the equivalent amplitude curves but plotting against the ratio of actual system oscillation frequency (and vortex shedding frequency under VIV) to vortex shedding frequency for a stationary cylinder. The displacement response is only plotted for the resonant $U_r$ range where the oscillation is periodic. Clearly, for each Reynolds number, there is an excellent collapse of the amplitude response for the wide range of mass ratios considered. As already commented, higher mass ratios result in the resonance dropping out at higher oscillation amplitudes as $U_r$ is increased; however, the shape of the curves remains the same. The only slight variation where this does not occur is in the neighbourhood of the maximum amplitude response for $Re=300$, where there are differences in the maximum response amplitude for the different mass ratios. The upshot of this collapse suggests that the transition scenario discussed in the previous section should hold for other mass ratios by using the post-simulation measured oscillation frequency as the scaling parameter rather than $U_r$. To this end, table 2 also provides the oscillation frequency for the transitions for $m^*=2.546$. Thus, to find the equivalent transition for, say, $m^*=1.0$, that occurs at $U_r=6$ for $m^{*}=2.546$ (initially to mode B at $Re=261.3$), simulations can be undertaken at $m^*=1.0$ for $Re=261.3$ varying $U_r$ to find when $f(D/U)= 0.1729$. For $m^*=1$, this occurs for $U_r=6.37$, and indeed determining the Floquet multiplier for this $Re$–$U_r$ combination for $\lambda /D = 1$ gives $\mu = 1.0$ in agreement with the $m^*=2.546$ result.

Figure 19. (a–c) Amplitude response for different Reynolds numbers and mass ratios. (d–f) Corresponding amplitude response against $f/f_v$ for the same cases in the left column.

3.8. Differences in 2-D and 3-D amplitude responses

In order to show the differences between 3-D and 2-D VIV simulations, even at relatively low Reynolds numbers, direct 3-D numerical simulations were undertaken. Figure 20 illustrates the comparison between 2-D and 3-D simulations at ${\textit {Re}} = 300$ for $m^* = 2.546$. This particular Reynolds number choice is interesting because the previous transition analysis shows that the wake should be 3-D except for the highest-amplitude cases. Indeed, the predicted maximum amplitudes for the 3-D simulations centred in the neighbourhood of $U_r=4$ agree for predictions from 2-D simulations. A check on these cases does indicate that the wake remains 2-D there. However, there are significant differences between the amplitude response at smaller and larger reduced velocities. In particular, the predicted 3-D amplitudes between $4.5 \lesssim U_r \lesssim 8$ are significantly larger than their 2-D counterparts and the response branch is beginning to show the typical response observed in high-Reynolds-number experiments (where the lower branch maintains an almost constant amplitude until finally dropping at higher reduced velocity after desynchronisation (e.g. Khalak & Williamson Reference Khalak and Williamson1997). Although less noticeable, the initial branch amplitude is also over-predicted in the 2-D simulation. Both the 2-D and 3-D simulations fail to predict the upper branch, where the response amplitude is largest.

Figure 20. Amplitude response for 2-D and 3-D simulations with respect to $U_r$ at ${\textit {Re}} = 300$ for $m^* = 2.546$.

3.9. Amplitude response as the wake becomes increasingly turbulent

We also used 3-D simulations to examine the response curves as the flow becomes increasingly turbulent, especially to document the appearance of the upper branch and to further determine differences from 2-D simulation predictions. The response curves are shown in figure 21 for $Re=300$, 500 and 1000. These predictions are compared with the experimental response curve for $m^*=2.4$ of Khalak & Williamson (Reference Khalak and Williamson1997) for $2000 < Re < 12\,000$, and numerical simulation results from Zhao et al. (Reference Zhao, Cheng, An and Lu2014) at a slightly smaller mass ratio of $2$. There is excellent agreement between the current predictions and those of Zhao et al. (Reference Zhao, Cheng, An and Lu2014), noting the effect of the smaller mass ratio is to extend the reduced velocity range of the lower branch slightly. Both sets of simulations predict a noticeable upper branch. The amplitude in the upper branch is below that seen in experiments, which is consistent with the previous observed dependence of the upper branch amplitude with Reynolds number for VIV (Govardhan & Williamson Reference Govardhan and Williamson2006; Rajamuni et al. Reference Rajamuni, Thompson and Hourigan2018). On reducing the Reynolds number to $Re=500$, the upper branch disappears. This was also observed by Blackburn et al. (Reference Blackburn, Govardhan and Williamson2001) for simulations for $Re = O(500)$. That study focused on demonstrating that shedding in the lower branch was of type $2P$, consistent with experiments unlike the typical $2S$ wake seen over most of the resonant range in lower-$Re$ simulations.

Figure 21. Amplitude response for 2-D and 3-D simulations with respect to $U_r$ at ${\textit {Re}} = 300$ for $m^* = 2.546$.

In terms of the lower branch wake state observed in the current set of simulations, figure 22 provides a comparison of wake vortices at $Re=300$ for the 2-D (a) and 3-D (b) simulations for $U_r=7.5$. Thus, development of 3-D flow changes the wake from $2S$ (for 2-D) to $P+S$ (for 3-D), which can be associated with the higher lower-branch response amplitude in the latter case. On increasing the Reynolds number to $Re=500$, the lower-branch wake changes to a $2P$ mode, as shown in figure 23 and consistent with Blackburn et al. (Reference Blackburn, Govardhan and Williamson2001). The top two images display spanwise-averaged wake vorticity contours for $Re=300$ (left) and $Re=500$ (right). The bottom row of images shows a top down view of these wakes visualised using isosurfaces of constant $Q$-criterion ($Q=0.1$), which highlights rotational vorticity rather than shear. Clearly, the $Re=300$ wake is much less complex, although there is still significant distortion of the spanwise rollers. At $Re=500$, the wake is effectively turbulent while still showing the remnants of the instability modes examined in this paper.

Figure 22. Effect of 3-D transition on the lower branch. (a) Spanwise vorticity contours at $Re=300$, $m^*=2.546$, $U_r = 7.5$ for 2-D flow. (b) Same for 3-D flow at a representative spanwise cross-section.

Figure 23. Mode selection for the lower branch as Reynolds number is increased. (a) Span-averaged spanwise vorticity contours at $m^*=2.546$, $U_r = 7.5$ and $Re=300$ for 3-D flow. (b) $P+S$ wake mode for this case depicted using isosurfaces of $Q$-criterion ($Q = 0.1$). (c) Same as (a) for $Re=500$, showing the mode has switched to $2P$. (d) Corresponding $Q$-criterion contours for this case.