Numerical Simulation of Cross-Flow Vortex-Induced Vibration of Hexagonal Cylinders with Face and Corner Orientations at Low Reynolds Number

: Vortex-induced vibrations (VIV) of hexagonal cylinders at Reynolds number of 1000 and mass ratio of 2 are studied numerically. In the numerical model, the Navier–Stokes equations are solved using finite volume method, and the fluid-structure interaction (FSI) is modelled using Arbitrary Lagrangian Eulerian (ALE) Scheme. The numerical model accounts for the cross-flow vibration of the cylinders, and is validated against published experimental and numerical results. In order to account for different angles of attack, the hexagonal cylinders are studied in the corner and face orientations. The results are compared with the published results of circular and square cylinders. Current results show that within the studied range of reduced velocities (up to 20), unlike circular and square cylinders, no lock-in response is observed in the hexagonal cylinders. The maximum normalized VIV amplitudes of the hexagonal cylinders are 0.45, and are significantly lower than those of circular and square cylinders. Vortex shedding regimes of the corner-oriented hexagons are mostly irregular. However, in the face-oriented hexagons, the shedding modes are more similar to the typical P + S and 2P modes.


Introduction
Vortex-induced vibration (VIV) is a major concern in design of slender structures exposed to the wind flow, such as: high-rise buildings [1,2] and bridges [3][4][5], or subject to water currents, such as: offshore platforms [6,7] marine risers [8][9][10] mooring elements [11][12][13], free spanning pipelines [14][15][16], or under action of other types of flow such as: heat exchanger tubes [17][18][19] and reactors [20,21]. As the flow passes around the bluff body, vortices are created in the wake region. These vortices generate pressure fluctuations because of a phenomenon known as "vortex shedding" that can cause bluff bodies to vibrate [22]. Vibrations may lead to fatigue damage in structures that are exposed to dynamic loads or flow-induced instabilities [23]. In the case of vortex-induced vibration, vibrations can be in two dimensions because of time-dependent drag and lift forces, as the vortices are shed from the bluff body [24].
Two fundamental experimental studies on cross-flow VIV of circular cylinders with one degree of freedom (1-DOF) were conducted by Feng [25] and Khalak and Williamson [26]. These pioneering studies showed that the VIV response is significantly related to the mass ratio, m * and the reduced velocity Vr. Mass ratio is defined as m * = m/md (the structural mass m divided by the displaced fluid mass md). The reduced velocity is defined as Vr = U/fnD, where U is the upstream velocity, D is the cylinder diameter, and fn is the natural frequency of the structure (in air or in the water). Feng's [25] study was conducted at very high mass ratio (m * = 248) and showed two branches of response, namely, the initial and the lower branch. In the works done by Williamson et al. [26][27][28] and at low mass ratios (m * < 10), in addition to the previously reported branches, an upper branch was observed. These experimental studies showed that the vortex patterns in the initial branch comprised of "2S" modes (two single vortices in one pair per vibration cycle), while those of upper and lower branches were "2P" (two vortex pairs are shed in each cycle). Moreover, it was understood that in 1-DOF, the amplitude of response in the upper branch gets as large as the diameter of the cylinder. It was shown that in the upper and lower branches the VIV frequency coincides with the vibration frequency of the cylinder, a phenomenon known as "lock-in." Limited research conducted on 2-DOF flow-induced vibration of circular cylinders [29,30], showed transverse vibration amplitudes of 1.5 times the diameter and 2P vortex mode.
The VIV response of circular cylinders has been investigated using numerical simulations as well [31][32][33]. Most of these studies use Reynolds-averaged Navier-Stokes (RANS) equations [34], and model an elastically mounted cylinder subjected to flow using two-dimensional simulations. The numerical simulations showed that the transition between the initial and upper branches is hysteretic. Also, in the transition from the upper to lower branch, a change in phase ϕ = 180° (the delay between the lift force and the cylinder response) is observed [35]. Prasanth et al. [36] studied the effect of blockage ratio (B) on the VIV response of a circular cylinder with m * = 10 at Re < 150 using finite element formulation. They carried out two sets of computations with B = 1% and 5%, where blockage B = D/H, is defined as the ratio of the diameter of the cylinder (D) to the distance between the lateral boundaries at the inlet (H). Their results showed that the hysteretic behavior between the lower and upper limits is developed at B = 5%, and vanishes at B = 1%.
In comparison with the VIV of circular cylinders, fewer studies are available for the VIV of bluff bodies with sharp edges. For a square cylinder, most noticeable are the experimental works of Bearman et al. [37], Amandolese and Hemon [38] and Zhao et al. [39], and the numerical study of Corless and Parkinson [40] and Zhao et al. [41]. The experimental study by Zhao et al. [39] of a cylinder with m * = 2.4 in water flow at angles of attack of α = 0° (square orientation) and α = 45° (diamond orientation), showed a galloping response in the square orientation and a similar response to that of a circular cylinder in the diamond orientation. Zhao et al. [41] conducted numerical simulations of VIV of a square cylinder with m * = 3, Re = 100, and at angles α = 0°, 22.5° and 45°, using Petrov-Galerkin finite element method. They reported 2S vortex shedding mode in all reduced velocities (1 < Vr < 30) for α = 0°. At α = 45°, P + S mode and 2P modes were observed, with a maximum normalized cross-flow vibration amplitude of 0.92 at Vr = 10.
Unlike square cylinders, the vortex-induced vibrations of polygonal cylinders have been marginally studied in wind engineering such as in bridge decks Kawatani et al. [42] or in marine structures such as bundled pipelines [43]. Limited studies are available on vortex-shedding from hexagonal cylinders [44][45][46]. This study aims to investigate the cross-flow VIV response of hexagonal cylinders with m * = 2 at Re = 1000. To do so, numerical models of one-degree-of freedom (1-DOF) vibration of the hexagonal cylinders based on the solution of two-dimensional Reynolds-averaged Navier-Stokes (RANS) equations and k-ω turbulence model are developed. As such, the in-line vibration and flow-induced instabilities such as galloping are not investigated herein. The lift and cross-flow vibrations as well as vortex shedding regimes at reduced velocities between 2 and 18 of face-oriented and corner-oriented hexagonal cylinders are discussed.

Numerical Method: Governing Equations and the Computational Model
Incompressible Reynolds-averaged Navier-Stokes (RANS) equations are used to simulate the flow. To avoid the mesh distortion, the moving boundaries of the cylinder surface are modelled with the Arbitrary Lagrangian Eulerian (ALE) scheme. Assuming u and p represent the time-average values of velocity and pressure, the RANS equations in a Cartesian coordinate system are written as where ν is the viscosity and Sij is the mean stress tensor. The small-scale fluctuations of velocity relating to the turbulence are reduced as ' ' i j u u , which is referred to as the Reynolds stresses.
Menter's [47] shear stress transport (SST) k-ω model, known to accurately capture the adverse pressure gradient flows (because of the concept of elliptic relaxation of ω-equation in the near wall region) is used to model the turbulence in Equation (1). The mean stress rate tensor Sij, and the Reynolds stress tensor where νt = k/ω is the turbulent viscosity and k is the turbulent energy. Menter's [47] formulation uses a k-ω model in the inner part of the boundary layer and a k-ε model in the free shear layer, where k is the turbulent kinetic energy, ω is the dissipation rate, and ε is the isotropic dissipation of the turbulence energy. The turbulence kinetic energy is given by where I is the turbulence intensity defined by The specific dissipation rate (ω) is where l = 0.07D is the turbulence length and Cμ is an empirical constant approximately equal to 0.09 [48].
The computational mesh with a corner-oriented hexagonal cylinder is shown in Figure 1. The length of the computational domain in the in-line and cross-flow directions are 50D and 40D, respectively [41]. The cylinder is located 20D downstream of the inlet. A free-stream velocity of u = U, v = 0 and in the x-direction is defined at the inlet. The pressure gradient and gradient of the fluid velocity at the outlet reduce to zero and the symmetry boundary conditions are imposed to the top and bottom surfaces. The fluid velocity along the cylinder wall surface is assumed to be the same as the structural speed of vibration of the cylinder (no-slip boundary condition). The equation of motion of the cylinder in the cross-flow direction (Y) is where C and K are the structural damping and lateral stiffness constants of the cylinder, respectively and FL is the fluid force on the cylinder in the cross-flow direction (lift force). To get the most adverse vibration responses, the damping constant is set to zero in the current study. In the computational model, a moving mesh zone (dynamic zone) with a diameter of 5D around the cylinder is adopted, and is allowed to move with the cylinder in the cross-flow direction as a rigid body. Transient ANSYS-Fluent [49] package is used to solve the fluid-structure interaction model (Equations (1)-(6)) using an iterative approach. After each time-step, the cylinder and the dynamic zone move in the cross-flow direction, and the displacements of the nodes are calculated based on the Laplace equation [50] ( ) 0 where Sy represents the nodal displacement in the cross-flow direction, γ is a parameter used to avoid excessive deformation of the near-wall elements, and A is the area of the element. No-slip boundary condition on the cylinder surface is adopted [51]. Selection of the time-step in a numerical simulation is a compromise. Trying to cut the time step to achieve a better temporal view would of course increase the computational resources demand. Trying to do so on a coarser mesh would mean fine scale structures would not be captured. Moreover, to get a Courant number less than or equal to 1, a normalized time step of 0.02 is calculated from

Validation of the Numerical Model
Numerical results from VIV of a circular cylinder with m * = 2.0 and at Re = 1000 using the model described in the previous section are shown in Figure 2. The normalized amplitudes of vibration (A/D) are plotted against the reduced velocities in Figure 2a. Current results agree well with the published numerical 2D and 3D results in the initial (Vr < 4) and upper branches (4 < Vr < 12). At Vr > 12, the amplitudes calculated from the current and published 2D numerical results [53,54] are lower than the experimental [28] and 3D numerical [32] results.
As can be seen in Figure 2a, most of the numerical models fail to capture the large amplitudes observed in the experimental study of Khalak and Williamson [28] within the upper branch. There are some exceptions to this. In a recent effort by Nikoo et al. [54] (not shown in Figure 2), the RANS-SST k-ω method with low-Re correction was used, and amplitudes of vibration closer to the experimental results within the upper branch were calculated. Pan et al. [51] used a fourth-order Runge-Kutta scheme to integrate the equation of motion (Equation (6)), and used a RNG k-ε turbulence model to solve the governing equations of Equation (1). They report that by reducing the normalized velocity Vr step to increments of 0.025, a normalized amplitude as high as A/D = 0.7 can be reached at Vr = 4.2 [53]. However, they argued that this maximum amplitude did not correspond to a stable state [53]. The current study aims to demonstrate the overall response of hexagonal cylinders with respect to cross-flow vibrations, and to compare it against that of a circular cylinder. Thus, Vr steps of one is used herein.  [28]. At 4 < Vr < 12 (the upper branch), the cylinder vibrates at a constant frequency, which is almost equal to its natural frequency (lock-in region). Outside this region, the vibration frequency follows the Strouhal law. It should be noted that using a smaller reduced velocity step (0.025), Pan et al. [51] observed the lock-in phenomenon at Vr between 4.4 and 11 (similar to the current results).

Mesh Sensitivity Analysis
In order to obtain the optimum mesh size of the computational model, mesh sensitivity analyses of a rigid circular cylinder at Re = 100, and rigid circular and corner-oriented hexagonal cylinders at Re = 1000 were conducted. Results are represented in Table 1 for four different mesh types. The number of divisions per the length parameter (Li), defined in Figure 1, and the total number of elements for each mesh are given in Table 1. As represented in Table 1 at both Re = 100 and Re = 1000, reasonable accuracy is observed using mesh 3. It can be seen that the root-mean-square lift, mean drag and Strouhal number of mesh 3, agree well with the previously reported results of circular cylinders at Re = 100 [55][56][57]. At Re = 1000 and using mesh 3, Strouhal number of the corner-oriented hexagon is only 0.3% different from the value reported by Khaledi and Andersson [46], using a DNS approach [58]. Moreover, St and CD results obtained with mesh 3 at Re = 1000 agree well with previously reported results. Thus, mesh 3 is adopted herein, and is used to model the vibration response of hexagonal cylinders in corner and face orientations. Table 1. Mesh dependence study of circular and hexagonal rigid cylinders using the computational model of Figure 1.

One-Degree-of-Freedom Responses of the Hexagonal Cylinders
The VIV responses of circular and hexagonal cylinders with face and corner orientations are plotted in Figure 3 at reduced velocities from 2 to 12, at steps of one, and at steps of two for reduced velocities from 12 to 18. It should be noted that in a corner-orientation, the hexagon has smaller frontal area and therefore a lower corresponding Reynolds number (Re = 867). The difference in Re is reflected in the reduced velocity as well. Therefore, in the corner-oriented hexagon the reduced velocities do not correspond to the exact integers. The normalized amplitudes and frequencies are plotted in Figure 3. It is meaningful to compare the current amplitudes and frequencies with those of a square cylinder at flow approaching angles α of 0° and 45°. For this purpose, the two-degrees-offreedom response of a square cylinder with m * = 2.0 and Re = 100 conducted by Zhao et al. [41] is selected, and is plotted in Figure 3a.   Figure 3a. Current results show that compared to the corner-oriented square, at Vr < 16, the corner-oriented hexagon has diminished VIV response.
In the corner-oriented hexagon, the amplitude response is comprised of three branches: (a) A lower branch at 3 < Vr < 6, (b) an upper branch at Vr > 12, and (c) a transition zone at reduced velocities between 6 and 12. Normalized amplitudes range from 0.1 to 0.4 in the corner-oriented hexagon. In the face-oriented hexagon, normalized amplitudes as large as 0.2 are observed at reduced velocities as low as Vr = 2. At Vr > 3, the normalized amplitudes fluctuate between 0.3 and 0.45. At Vr < 9, the face-oriented hexagon has larger amplitudes than the corner-oriented hexagon. It is worth to mention that Kumar et al. [62] observed a new extended initial branch in the VIV response of elliptical cylinder at low Re and low mass ratio.
The normalized frequencies of the corner-and face-oriented hexagons are shown in Figure 3b,c, respectively. Vibration frequencies of the hexagonal and face-oriented square cylinders follow the Strouhal law (i.e., the cross-flow vibration frequency is equal to the vortex-shedding frequency) for the entire range of the reduced velocity, and no lock-in is observed. At Vr > 4, the normalized frequencies of the corner-oriented hexagon are above the Strouhal line, whereas, those of the faceoriented hexagon are below it. It should be noted that in the square study [41], lock-in response was only observed in the response of the corner-oriented square (α = 45°) and at 3 < Vr < 10. However, no lock-in response is seen in the hexagonal cylinders.
Time histories of the lift coefficients and vibration amplitudes of the corner-oriented hexagon cylinder are shown in Figure 4 for selected reduced velocities. The normalized time is hexagons is equal to D, and that of the corner-oriented hexagon is equal to √3 2 D based on a cylinder of unit length. As shown in Figure 4a, at small reduced velocities (Vr < 5), the mean value of the lift coefficient is zero throughout. At larger reduced velocities, although the cylinder configuration is symmetric, the mean value of the lift is non-zero. At Vr = 11.41, the non-zero lift fluctuates from one side of the cylinder to the other side. However, at Vr = 20.57, the mean lift has a negative value. Normalized amplitudes are shown in Figure 4b   The lift time histories of the face-oriented hexagonal cylinders are shown in Figure 5a. At Vr = 2 the mean value of the lift force is non-zero, and is on the positive side of the curve for the timespan considered. However, at Vr = 5 the non-zero lift changes from a positive mean value (32 < t * < 34) to a negative mean value (34 < t * < 36). At larger reduced velocities, the lift force becomes irregular. This irregular behavior is observed in the time history of the vibration amplitudes at larger reduced velocities (Vr = 10, 18), as well (Figure 5b). This phenomenon is related to the vibration modes and will be discussed later. The phase difference (ϕ) between the lift coefficient and the cross-flow displacement is determined by the FFT analysis of each cylinder at different reduced velocities, and is shown in Figure 6. The phase difference corresponds to the delay between the peak in the amplitude and the peak in the lift force component at the same frequency, and is calculated from the amplitude spectra of the lift force and cross-flow displacement (not shown herein). Phase differences of square cylinders are taken from [41], and are plotted as well. It should be noted that in the corner and face-oriented hexagons, the force and vibration responses are irregular at some reduced velocities (refer to Figures  4 and 5). Therefore, in calculating the phase angles, effort was made to capture a time domain that corresponded to a regular response, when possible. The most distinctive difference between the phase responses of the hexagonal cylinders and the rest, is the emergence of the out of phase response ( 0    ) at low reduced velocities. Similar behavior was reported in the response of a faceted cylinder (textured pipe) [63]. At reduced velocities of 5 and above, the lift force becomes out of phase ( 180    ) with the cross-flow amplitude. The out of phase response was reported to occur at Vr =7 and Vr =9, in the flat and corner-oriented squares, respectively ( Figure 6). Fluctuations in the phase angles of the hexagonal cylinders are observed at 10 < Vr <16. This is more evident in the face-oriented hexagon, where a phase angle of ϕ = 165° is seen at Vr = 14.  Figure 7 shows the vorticity contours of the corner-oriented hexagon at Vr = 2.28, at four instants in one vibration period during which the normalized vibration frequency is 0.346 and the normalized amplitude is 0.03. The trajectory of the cylinder in the xy plane, and its instant position are depicted next to each vortex contour. The legends are provided for comparison between the vorticity magnitudes. Vortex C is shed from the top of the cylinder as it is moving from the top to the bottom of the trajectory. At this instant two positive vortices (A and B) are already shed from the bottom of the cylinder, and are located between 2-4 D in the wake. As the cylinder moves down the trajectory, the vortices A and B merge and move further down the wake. While the cylinder is at its lowest trajectory position, three vortices (one positive and two negative) are being shed. At the end of the vibration period, vortices A and B are merged into a single vortex, and are followed by vortex C. As can be seen in Figure 7, the vortex regime of the corner-oriented hexagon at this small reduced velocity changes from a P + S mode (pair A + B and single C), to a 2S mode (C and A + B). Vortex regimes of the corner-oriented hexagon at Vr = 11.41, corresponding to frequency and amplitude of 2.4 and 0.35, respectively, are shown in Figure 8. One pair (B and C) and one single (A) vortex are shed in each vortex shedding period, describing a P + S mode. The shedding frequency is twice the vibration frequency. When the cylinder moves down the trajectory, the P + S is shed from the top, and when the cylinder moves upward, a new P + S is generated from the bottom. The P and S vortices have different sizes but similar strengths. As predicted from Figure 3b, the vibration frequency is smaller than the Strouhal frequency.  Figure 9 shows the vortex regime of the corner-oriented hexagon at Vr = 20.57. At this reduced velocity, the vibration frequency is significantly lower than the Strouhal frequency and is 3.5 times larger than the natural frequency of the cylinder. The vibration amplitude (Figure 3a) at Vr = 20.57 is not much different from the amplitude at Vr = 11.41. However, the vortex regimes at those reduced velocities are dissimilar. As shown in Figure 9, four vortices (labelled by A-D) are shed from the cylinder in each vibration period. These vortices have different sizes and magnitudes and do not pair up behind the cylinder. Zhao et al. [41] reported similar vortex regime in a square cylinder at α = 22.5°, with Vr = 7 (outside the lock-in region).  The P + 2S formation observed here is significantly different from that of the corner-oriented hexagon at the same reduce velocity. The four vortices progress slowly in the rear of the cylinder (almost 1D in a vibration period). The pair (A + B) has larger vorticity and is shed while the cylinder moves down the trajectory. The single vortices (C and D) are not as strong and move below y = 0 line, throughout the wake. This P + 2S mode is different from typical mode. As seen in Figure 10, another single vortex of smaller vorticity is detached from the pair (A + B), at half-period. The vortex regime of the face-oriented hexagon at Vr = 10, is shown in Figure 11, and clearly depicts a P + S mode in one vibration period. The already generated single vortex at the bottom (A), moves in the bottom rear wake of the cylinder, while a pair (C + B) is shed from the top. It can be seen that at the end of the period, a new pair of vortices are being developed at the bottom of the cylinder. Comparison between Figure 11 and Figure 8 at (10 < Vr <12) shows that in a P + S vortex regime mode, the amplitudes and frequencies of vibration of the hexagons in corner and face orientations are almost equal. At similar amplitudes (0.2 < A * < 0.3), Zhao et al. [41] reported a 2S mode for a circular cylinder with α = 0°. Figure 11. Vorticity contours at four instants in a vibration period of the face-oriented hexagonal cylinder at Vr = 10 at 21.32 < t* < 21.92. Figure 12 shows the vortex regime of the face-oriented hexagon at Vr =18. Four vortices in two pairs (2P mode) are shed within each vibration cycle. The current regime is a typical 2P mode, and the two pair of vortices that are shed from the cylinder in each vibration period, are symmetric about an imaginary y = 0 line. Toward the end of the vibration cycle, and at a distance of 8D in the rear wake, the pair of vortices still have high vorticity. This explains the high amplitude and frequency of vibration, and the irregular lift and amplitude responses of the face-oriented hexagon at large reduced velocities.

Conclusions
One degree of freedom VIV response of hexagonal cylinders in corner and face orientations at Re = 1000 and at low mass ratio of m * = 2 were investigated using numerical simulations at reduced velocities up to 20. Results were compared against a published numerical work on VIV of square cylinders at low Reynolds number and mass ratio [41]. Unlike square cylinders, no lock-in region was observed in the VIV responses of the hexagonal cylinders at the studied reduced velocities. Normalized frequencies (vibration frequency divided by the natural frequency) of the corneroriented hexagon were found to be dominantly below the Strouhal frequency (fs), whereas, in the face-oriented hexagon, the normalized frequencies were slightly larger than fs.
Within the studied reduced velocity range, the maximum normalized cross-flow amplitude of vibration of the hexagonal cylinders were shown to be around 0.45. This is significantly lower than the reported maximum amplitudes of 0.96 for circular cylinder [28] and 0.925 for corner-oriented square cylinder [41]. At large reduced velocities, similar amplitudes (around 0.4) were observed in the corner-oriented square and hexagonal cylinders. Unlike other cylinders, phase differences between amplitude and lift forces ( 0    ) were seen in the hexagonal cylinders, at low reduced velocities. The phase difference settled at 180    in the corner-oriented hexagon at large reduced velocities. However, fluctuations were seen in the phase differences of the face-oriented hexagon.
The vortex shedding regimes of the face-oriented hexagons were more similar to the typical modes reported in the literature, compared with the corner-oriented hexagon. At low reduced velocities (and low amplitudes), P + S modes were observed in the corner-oriented hexagon. At large amplitudes, four non-pairing vortices of different sizes and magnitudes were observed in the corneroriented hexagon, similar to observations of Zhao et al. [41] of a corner-oriented square at large amplitudes. P + 2S modes were observed in the face-oriented hexagon, and explained the existence of sizeable vibration amplitudes at low reduced velocities. At larger reduced velocities, P + S modes were observed, and clarified the larger amplitudes of face-oriented hexagon compared to flat square cylinders with 2S modes reported in [41]. At very large reduced velocities distinct 2P modes were seen, and were shown to have significant vorticity magnitudes in a distance as far as 8D in the rear wake of the face-oriented cylinder.