Numerical Investigation on Vortex-Induced Rotations of A Triangular Cylinder Using An Immersed Boundary Method

A numerical study of vortex-induced rotations (VIRs) of an equivalent triangular cylinder, which is free to rotate in the azimuthal direction in a uniform flow, is presented. Based on an immersed boundary method, the numerical model is established, and is verified through the benchmark problem of flow past a freely rotating rectangular body. The computation is performed for a fixed reduced mass of m*=2.0 and the structural stiffness and damping ratio are set to zero. The effects of Reynolds number (Re=25–180) on the characteristics of VIR are studied. It is found that the dynamic response of the triangular cylinder exhibits four distinct modes with increasing Re: a rest position, periodic rotational oscillation, random rotation and autorotation. For the rotational oscillation mode, the cylinder undergoes a periodic vibration around an equilibrium position with one side facing the incoming flow. Since the rotation effect, the outset of vortex shedding from cylinder shifts to a much lower Reynolds number. Further increase in Re leads to 2P and P+S vortex shedding modes besides the typical 2S pattern. Our simulation results also elucidate that the free rotation significantly changes the drag and lift forces. Inspired by these facts, the effect of free rotation on flow-induced vibration of a triangular cylinder in the in-line and transverse directions is investigated. The results show that when the translational vibration is coupled with rotation, the triangular cylinder presents a galloping response instead of vortex-induced vibration (VIV).


Introduction
Vortex-induced rotation (VIR) is a common phenomenon that we observe in nature and is closely related to many engineering fields, ranging from the torsion of suspension bridges in wind to the operation of marine current turbines. VIR is a double-edged sword: it not only induces structural instabilities, but also provides energy resource for human beings through rotation energy harvester. As a fluid-structure interaction (FSI) problem, the VIR behavior is quite complicated and has received considerable attention.
The study of this topic can date back to as far as the middle of 19th century when Maxwell (1853) first analyzed the motion of a falling card. Since then, researchers have attempted to develop a better understanding of the VIR re-sponse. Riabouchinsky (1935) studied the autorotation of a four-bladed rotor (autorotation is defined as continuous rotation of a freely rotatable object in an azimuthal direction in cross flow). Lugt (1980) conducted numerical simulations to study the autorotation of an elliptic cylinder. More works on the subject of autorotation were reviewed in Lugt (1983). Tatsuno et al. (1990) experimentally studied a freely rotatable square cylinder in a uniform flow at Re=3.5×10 4 . Zaki et al. (1994) recognized four distinct responses of a freely rotatable square cylinder in a range of Reynolds numbers from 10 3 to 10 4 : (i) a rest stable position, (ii) oscillation, (iii) reverse rotation and (iv) autorotation. In order to illustrate the flutter instability of the old Tacoma Narrows Bridge in the1940s, Hübner et al. (2001) proposed a model of H-shaped bridge deck which experienced coupled motion of vertical translation and rotation. Mittal et al. (2004) studied the effect of Reynolds number and thickness-tolength ratio on the autorotation of a plate. Srigrarom and Koh (2008) experimentally studied the self-excited rotational oscillation of an equilateral triangular cylinder. This rotation exists when the Strouhal number is in the range of 0.13<St<0.18. Lu et al. (2016) conducted a numerical study on the flow-induced rotary oscillation of circular cylinder with rigid splitter plate at Re=100. The results showed that the natural frequency of the rotating body in flow led to a wide regime of lock-in due to the combined effect of fluid moment, rotation response and phase difference. Ryu and Iaccarino (2017) numerically studied the vortex-induced rotations of a rigid square cylinder at 45≤Re≤150. The parametric investigation revealed six different dynamic responses of the cylinder and their coupled vortex patterns. More recently,  studied the vortex-induced rotation of an elliptic cylinder about its center with a variety of torsional frictions. The numerical results indicate that the elliptic cylinder can rotate with too large torsional friction or too small reduced velocity, while too small friction leads to transverse galloping instead of the desynchronized region.
As can be seen, the VIR characteristics are significantly influenced by various factors, such as Reynolds number, cross-section shape, structural natural frequency damping and so on. To the best of our knowledge, research on VIR is far from enough. On the other hand, previous work has shown that the rotation of non-circular cross section can give rise to the shift of separation point, resulting in drastic change of flow field and flow-induced forces on the body (Ryu and Iaccarino, 2017;Tu et al., 2014). Now comes the question: how the classical vortex-induced vibration (VIV) system will behave if it is added a rotational degree of freedom. This issue has not been studied systematically.
In this paper, vortex-induced rotations of an equilateral triangular cylinder in a uniform flow are numerically investigated, which is of important significance in physics and engineering applications (Srigrarom and Koh, 2008). Our objective is targeted at studying the VIR characteristics within a wide range of Reynolds number, as well as the effects of VIR on VIV. Simulations are performed by employing an immersed boundary (IB) method. Compared with conventional arbitrary Lagrangian-Eulerian (ALE) schemes (Zhu and Lin, 2018;Wang et al., 2019), the IB method has significant advantages, especially in the FSI simulations with topological changes.
The remainder of this paper is organized as follows. Section 2 presents the numerical model and its validation. In Section 3 we outline the setup of the problem. Section 4 provides the numerical results and discussions. Finally, concluding remarks are drawn in Section 5.

Governing equations
The fluid flow is described by the continuity and the Navier-Stokes equations, assuming that the viscous Newtonian fluid is incompressible and homogeneous. The governing equations for fluid dynamics read as: where u is the velocity, t is the time, is the density of the fluid, p is the pressure, f denotes the body force, and is the kinematic viscosity associated with Reynolds number, defined as . U ∞ is the free stream velocity and D is the characteristic length.
Under the excitation of the alternate fluid forces due to the vortex shedding, the rigid cylinder in this study is only allowed to rotate freely with respect to its centroid in the azimuthal direction. Thus the structural dynamic equations can be described as: I θ τ where I is the moment of inertia, is the rotation angle, c is the damping coefficient, k is the stiffness coefficient, and is the torque.

Numerical method
To simulate the flow, we used an in-house computational fluid dynamics (CFD) C code called CgLes (Thomas and Williams, 1997;Ji et al., 2013), which has a proven high fidelity and parallelizing efficiency. In this code, the Navier-Stokes equations are discretized on a fixed staggered Cartesian grid and solved by the finite volume approach.
The immersed boundary method, first introduced by Peskin (1972) for the simulation of blood flow around the flexible leaflet of the human heart, has been incorporated into CgLes to model the fluid-structure interaction. By using the IB method and the second-order Adams-Bashforth time marching scheme, the Navier-Stokes equations can be discretized and rearranged as: comprises the convective and diffusive terms. The purpose of the IB method is to set the boundary conditions on the immersed boundary exactly. This can be expressed as: where U represents the velocity on the IB points interpolated from Cartesian grids, and V is the desired velocity on the IB points obtained by solving the governing equation of cylinder motion. The body force on the Cartesian grids is finally calculated by: where and are the interpolation and distribution functions, respectively. Note that represents the variables on the Cartesian grid x, while indicates the variables on the IB points X. Eqs. (4) and (6) are then solved iteratively in an implicit way until the boundary condition shown in Eq. (5) is satisfied. For the sake of conciseness, details of the iterative IB method are not presented here and readers can refer to previous work (Ji et al., 2012) for more information.
Through evaluating the hydrodynamic forces acting upon the cylinder by means of momentum balance over the corresponding fluid domain, the structural dynamic equations in the IB scheme can be written as: where x c is the center coordinates of the cylinder, N IBP is the total number of the IB points, and F(X i ) and ΔV i are the extra body force and discrete volume of the IB point X i . M denotes the total amount of torque added to the fluid, and R was evaluated as a sum over each grid cell with the cells volumetric solid fraction as a weight. In the case of evaluating the rate-of-change term from a full rigidity approximation, Eq. (7) is rearranged as: ρ s where is the density of the solid. Eq. (8) is solved by using Newmark's (1959) algorithm.

Validation test
First, we verify the above numerical model by applying it for VIR of a rectangular cross-section body. This case has been numerically studied by Li et al. (2002) using a moving frame of reference algorithm. Fig. 1 shows the computational domain, discretized by a non-uniform Cartesian grid with a resolution of 768 × 384. A uniform grid with the spacing of Δx/D = Δy/D = 1/40 was applied in a region of 16D × 8D encompassing the cylinder. Properties of the mechanical system are specified below: Reynolds number 250, reduced velocity U R = U ∞ /(f N D) = 40 (f N is the structural torsion natural frequency), moment of inertia ratio Ik) = 0.6 , and damping ratio . To satisfy the courant-Friedrich-Levy (CFL) condition, the time increment is set to ΔtU ∞ /D = 0.005. Fig. 2 presents the vorticity field at the time point when the rectangular body approaches its maximum rotation. The time histories of the rotation angle and moment coefficient are plotted in Fig. 3, including both present and previous results (Li et al., 2002). It is seen that the rotational vibration takes a stable pattern after a while. The numerical results obtained here agree reasonably well with those of Li et al. (2002), and the maximum difference between two results of the rotation angle is smaller than 1.0%. This test case demonstrates that the present numerical model is adequate enough to settle the VIR problem.
In the next step, the coupled VIV and VIR response of an elliptic cylinder (Wang et al., 2019), is also simulated to assess the accuracy of the present numerical code. The computational configuration is described in Fig. 4, which has the same mesh system as the first validation case. The ellipse can be free to vibrate in both transverse and azimuthal directions. The area of the ellipse is , where a and b denote the x and y semi-axes, respectively. The ellipse's equivalent diameter is , and the aspect ratio is expressed as b/a=2.0. The Reynolds number based on the equivalent diameter is Re=150. The mass ratio of the ellipse is 10.0. Fig. 5 displays the transverse and rotational vibration amplitudes of the cylinder with respect to the reduced velocity. A favorable agreement is found between our results and the previous data.
After the validation of the numerical code, attention is firstly focused on the simulation of VIR of an equilateral triangular cylinder in a uniform flow, as shown in Fig. 6a. The two structural parameters, stiffness and damping coefficients, are set to be zero. The mass ratio of the cylinder is (m is the cylinder mass and D is the side length). The moment of inertia ratio is computed by I * = m * h 2 /(12D 2 ). Numerical simulations are conducted in a range of Reynolds numbers between 25 and 180. Further-more, in order to investigate the rotation effects on VIV response, calculations are also carried out for a triangular cylinder with three degrees of freedom (3DOF, in-line, transverse and azimuthal) at a typical Reynolds number of Re=100, as shown in Fig. 6b. The reduced velocities U r are identical in both in-line and transverse directions and varied from 2 to 18. In the vicinity of the cylinder, a uniform grid with the spacing of Δx/D=Δy/D=1/40 is arranged in a region of 16D×8D. The time step is set to ΔtU ∞ /D=0.005. A Dirichlet boundary condition is imposed upon the inlet boundary as u 1 =U ∞ and u 2 =0, while at the outlet the condition is p=0. A no-slip condition is employed on the cylinder surface. The top and bottom boundaries have a slip condition of and u 2 =0.

Flow configuration
A mesh independence test is performed on three mesh systems, i.e. coarse, medium, and fine, with the minimum mesh sizes of Δh/D = 1/16, 1/40 and 1/64, respectively. The tested parameters and the computed results are listed in Table 1, along with the percentage changes. The maximum    vibration amplitudes X max /D, Y max /D and , the mean drag coefficient C D,mean , and the root mean square (r.m.s.) value of the lift coefficient C L,rms are used to inspect the mesh quality. It is found that the maximum difference is 6.95% in the value of Y max /D when the mesh changes from the coarse to the medium, whereas no significant change is observed when comparing the medium and fine meshes (the maximum difference is 1.55% in the value). This test indicates that the medium mesh can guarantee the accuracy of the numerical results.

Results and discussions
4.1 Rotation modes ∆θ max π/2, 2π 5π π/2 2π Fig. 8 shows the time evolutions of rotation angle (is accumulated in counterclockwise direction) of the triangular cylinder for different Reynolds numbers. The variations of the maximum amount of rotation with respect to Re are shown in Fig. 9. It is clear that the rotation responses present four distinct modes: a rest position, periodic rotational oscillation, random rotation and autorotation. With the different cross-section shapes and Reynolds numbers, the four rotation modes are qualitatively consistent with those reported by Zaki et al. (1994) and Ryu and Iaccarino (2017) for a square cylinder. This is mainly because the rotational motion also depends on the moment of inertia in addition to the geometry and flow conditions (Ryu and Iaccarino, 2017). A little difference is that the present results exhibit three different periodic oscillations with asymptotic limits and , as shown in Fig. 9. Remarkably, the random rotation mode occurs between the periodic oscillation regimes with asymptotic limits and . The maximum amount of rotational oscillation is found to be closely related to the Reynolds number. As shown in Fig. 8, the triangular cylinder finally remains at rest at Re=25 because no vortex shedding occurs and the rotational motion of the cylinder is not excited (corresponding flow patterns will be shown in Section 4.3). This Reynolds number is considerably smaller than the critical value of Re cr ≈39, at which the steady and symmetric vortex patterns appear for flow over a fixed triangular cylinder (De and Dalal, 2006;Zielinska and Wesfreid, 1995). At 30≤ Re≤120, the rotation responses exhibit sinusoidal oscillations with small amplitudes. With the increase of Reynolds number in this regime, the maximum amount of rotation monotonically increases to the asymptotic limit , as shown in Fig. 9. This regime is marked as 'small-amplitude oscillation'. At Re=130, intermittent oscillations and rotations in a random fashion appear in the response curve, where the maximum amount of rotation cannot be defined precisely. We refer to this regime as 'random rotation'. However, such behavior was not found in previous experiments on self-excited rotation of a triangular cylinder (Srigrarom, 2003;Srigrarom and Koh, 2008), because only a few Reynolds numbers were tested in their studies. As Re increases to 140 and 150, the rotation mode recovers periodic oscillation again but has its own unique characteristics from that of previous regime. First, the cylinder rotates with large amplitude within asymptotic limits and . Second, the response curves in the 'small-amplitude oscillation' regime are sinusoidal, while they present a bell-like pattern here. When Re exceeds 150, the cylinder experiences clockwise autorotation. Instead of changing directly from oscillation to autorotation mode, the rotation response undergoes a transient state at tU ∞ /D < 10. The autorotation direction could depend on the path from the transient state to the dynamic equilibrium one. In the present simulation, only counterclockwise autorotation is observed.
The frequency analysis for the freely rotating triangular cylinder is shown in Fig. 10, where f v is the vortex shedding frequency, f CM represents torque frequency, is the rotation frequency, and St is the Strouhal number. It is evid- ent that the vortex shedding frequency of the cylinder seriously deviates from the Strouhal line, but is mainly dominated by the rotation frequency. At 30≤Re≤120, the three frequencies are found to be the same , corresponding to the sinusoidal oscillations with small amplitudes. At Re =140 and 150, f CM increases significantly, whereas and f v drop to a very low level. The nonidentity of f CM and leads to large rotational oscillation amplitudes. When Re exceeds 150, f v and increase and keep a stable level. The frequency relationships are found to be f CM /3 ≈ f v /2 ≈ , corresponding to the appearance of autorotation mode.

Dynamic equilibrium characteristics
In the previous section, we also note that for most Re, the triangular cylinder can adjust its orientation adaptively and finally undergoes a periodic rotational oscillation around an equilibrium position. Fig. 11 shows the variations of the averaged rotation angle corresponding to the dynamic equilibrium site. In the range of 30≤Re≤110, the averaged rotation angles are nearly identical and the value of is found to be about , i.e., the triangular cylinder rotates counterclockwise to a position with one side facing the incoming flow. Interestingly, at Re=120, 140 and 150, the averaged equilibrium angle presents dramatic fluctuations, such as , and . These values move up and down and seem in completely random patterns, but they can be reduced to , and by (n is an integer), which means that the triangular body obtains the same attack angle with one side facing the incoming flow after several circles of autorotation. Fig. 12 depicts the locations of the triangular cylinder approaching its dynamic equilibrium state, which is independent of the Reynolds number.

Vortex shedding modes
Previous studies have extensively concentrated on flow over the fixed and forced oscillating triangular cylinders (Alawadhi, 2013;De and Dalal, 2006;Iungo and Buresti, 2009;Tu et al., 2014), as well as flow-induced vibrations (Alonso et al., 2012;Wang et al., 2015). It is well known that the predominant vortex shedding pattern is the classic 2S mode (two single vortices are shed per cycle). However, if the cylinder is allowed to rotate freely in the azimuthal direction, the flow pattern may show different characteristics. Our results indicate that with the increase of Re, various flow modes, such as N (no vortex shedding), 2S, 2P (two vortex pairs are formed per cycle) and P+S (a single vortex and a pair of vortices are formed per cycle) patterns, appear behind the triangular cylinder. Some typical vorticity fields for flow past a freely rotating triangular cylinder are plotted in Fig. 13. Fig. 14 shows a summary of the distribution of vortex shedding modes for a freely rotating triangle. The results of a fixed counterpart are also provided for comparison. It is observed that the outset of vortex shedding from a triangular cylinder shifts to a much lower Reynolds number in the presence of rotational degree of freedom.    Fig. 15 shows the variations of the force coefficients of the triangular cylinder versus Re. As observed in Fig. 15a, in the subcritical range the mean drag coefficient C D,mean of the fixed cylinder decreases until Re reaches the critical value due to diminishing effect of viscous force (De and Dalal, 2006). Above Re cr , the drag coefficient presents an increasing function. In contrast, the C D,mean of the rotating one is significantly larger than that of the fixed case at Re<100. In this Re range, the cylinder undergoes a very small-amplitude rotational oscillation around the equilibrium position with one side facing the incoming flow, inducing a large-scale stagnation pressure pocket on the windward side. At Re>100, the C D,mean of the rotating cylinder is more or less reduced. This can be attributed to the largeamplitude rotational oscillation or autorotation which result in two physical changes: (i) asymmetric location of the stagnation pressure pocket on the windward side and (ii) a decrease of the shielding and suction effects on the leeward side.

Hydrodynamic forces
Owing to the Magnus effect, the root mean square (r.m.s.) lift coefficient C L,rms of the rotating cylinder rises to a higher level all over the Re range tested, comparing with the fixed case (as shown in Fig. 15b). C L,rms increases and obtains its peak in the small-amplitude oscillation regime, then falls to the bottom with the advent of random rotation mode. After that, C L,rms rises once again, which can be mainly attributed to the large-amplitude oscillation effect. Further increase in Re causes C L,rms to drop dramatically at the outset of autorotation, and subsequently remain stable. It is found that the C L,rms curve finally tends to collapse onto that of the fixed cylinder.

Rotation effect on translational VIV
In this section, the flow-induced motion of a triangular cylinder with three degrees of freedom (3DOF, i.e., in-line, transverse and azimuthal) is investigated at Re=100, where the lift fluctuation is significantly enhanced under the rotation effect (as shown in Fig. 15b). Fig. 16 shows the variations of the translational vibration amplitudes with reduced velocity U r . For comparison, previous results (Wang et al., 2015) for vibration amplitudes of a triangular cylinder with two degrees of freedom (2DOF, in-line and transverse) are also provided. It can be found that the 2DOF triangular cylinder clearly undergoes a VIV response and obtains its peak value in the lock-in region. The maximum vibration amplitudes in the in-line and transverse directions are X max /D=0.088 and Y max /D=0.241, respectively. By contrast, when the translational vibration is coupled with rota-tion, the triangular cylinder is subjected to galloping with large amplitude. The X max /D and Y max /D curves nearly have a monotonic increase with U r until reaching their peak values (X max /D =1.542, Y max /D =2.876) at the maximum reduced velocity U r =18. This behavior is very similar to that of a 2DOF (in-line and transverse) triangular cylinder with one side facing the incoming flow, for which readers can refer to Wang et al. (2015). It can be explained by reason that the free rotation of the triangular body can change its attack angle continuously to make the galloping response prone to happen. Fig. 17 shows the variations of the rotation angle for the 3DOF cylinder. Different from the periodic rotational oscillation at Re=100 in Section 4.1, the rotation responses here generally exhibit a random fashion, indicating that the coupled motion with 3DOF is highly nonlinear.

Conclusions
In this paper, vortex-induced rotation of an equilateral triangular cylinder is investigated numerically by using an immersed boundary method. Variations of rotation response, vortex shedding mode and hydrodynamic force are analyzed deeply within Reynolds number range of 25≤ Re≤180. Furthermore, simulations are also carried out for flow-induced vibration of a cylinder with 3DOF (in-line, transverse and azimuthal) at Re=100 to illustrate the rotation effect. The main findings can be summarized as follows.

π/2
The rotation responses of the triangular cylinder present four distinct modes: a rest position, rotational oscillation around the equilibrium position with one side facing the incoming flow, random rotation and autorotation. Below Re=25, no vortex shedding appears and the cylinder remains the initial rest position. In the periodic 'small-amplitude oscillation' regime (30≤Re≤120), with the increase of Re, the maximum amount of rotation monotonically increases to the asymptotic limit . For the 'large-amplitude oscillation' regime (Re=140 and 150), the cylinder ro-2π 5π tates within asymptotic limits and . Random rotation mode is observed between the two regimes above (Re=130). As Re exceeds 150, continuous rotation occurs in the clockwise direction.
Flow past a freely rotating triangular cylinder shows different characteristics compared with that for the fixed counterpart. Under the free rotation effect, the 2P and P+S modes appear behind the triangular cylinder in addition to the typical 2S mode, and the outset of vortex shedding from cylinder shifts to a much lower Reynolds number. The results also elucidate that the free rotation significantly changes the drag and lift forces.
When the vortex-induced vibrations in the in-line and transverse directions are coupled with rotation, the triangular cylinder experiences a galloping response with large amplitude. The translational amplitudes nearly have a monotonic increase with U r . This behavior is very similar to that of a 2DOF (in-line and transverse) triangular cylinder with one side facing the incoming flow. Because of the complex coupled motion, the cylinder rotation presents a random fashion.