Numerical Investigation of the Geometrical Effect on Flow-Induced Vibration Performance of Pivoted Bodies

: In this study, the Flow-Induced Vibration (FIV) of pivoted cylinders (at a distance) is numerically investigated as a potential source of energy harvesting. In particular, we investigate the effect of pivot point placement, arm length, and natural frequency on the FIV performance of six different cross sections in the Reynolds number of around 1000. All sections have similar mass, area, and moment of inertia to eliminate non-geometrical effects on the performance. Classical studies show that the synchronization phenomenon (lock-in) occurs when the vortex formation frequency is close enough to the body’s natural frequency. Due to the conﬁguration of the cylinder in this research (pivoted eccentrically), the natural frequency is also a function of the ﬂow velocity as well as the geometrical speciﬁcations of the system. The simulation is done for the arm lengths between − 3D and +3D for all cross sections. Results show that maximum output power is principally inﬂuenced more by the pivot location than the arm length. Although the box cross section has a higher amplitude of vibration, the circular cross section has the highest efﬁciency followed by the egg shape.


Introduction
In recent years, the development of flow-induced vibrations (FIV) energy harvesters has increased rapidly to offer a new source of energy. Due to the large strains and geometric deformations during FIV, they have traditionally been classified as a destructive phenomenon. One of the well-known examples of flutter-induced destruction is the Tacoma Narrows Bridge collapse in 1940, where torsional flutter at sufficiently large amplitudes caused catastrophic failure of the entire bridge. However, common and accessible FIV could be considered as a way to extract energy. Bernitsas et al. [1] have developed a device that uses the vortex induced vibration (VIV) phenomenon to generate electricity. Contrary to the VIV phenomenon, where significant oscillations develop in a small range of flow velocities and with limited oscillation amplitudes, other aeroelastic instabilities like flutter occur for an infinite range of flow velocities and without a self-limited response beyond the critical flow velocity which makes the flutter more promising for generating energy.
For instance, Hobbs and Hu [2] tested micro-watt energy harvesters inspired by tree trunks swaying in the wind. Their converter consists of four pivoted cylinders which affixed to the ground via a piezoelectric transducer. Yoshitake et al. [3] generated minuscule amounts of energy, using a device composed of Hula-Hoops and an electro magnetic transducer mechanism, in air flow. To study the aerodynamic efficiency of a drag assisted energy-harvesting device, Sung et al. [4] investigated the effects of the cylinder cross-sectional shape on the VIV. Their numerical simulations have demonstrated that an elliptical cylinder undergoes much larger displacements than a circular one. Nevertheless, their research mainly focuses on improving maximum displacement or amplitude rather than the angular velocity of the vibration.
In an attempt to study the performance of FIV, Arionfard and Nishi [5] carried out an experiment on a pivoted cylinder instead of a transitionally moving one. Being assisted by the drag force, the pivoting cylinder showed an increase in performance comparing to transitional VIV of a cylinder. As a result, different configurations with one and two cylinders were considered in the following researches to increase the performance by utilizing different mechanisms of vibrations [6][7][8]. However, an important way to improve the FIV performance is through the geometry of the bluff body and enhancement of the geometrical parameters is necessary in order to increase energy extraction performance. For several years great effort has been devoted to the study the effect of cross sections on FIV. However, common cross sections in aviation and civil engineering has attracted much more attention; Airfoil flutter [9], galloping of square, triangular, and semicircular sections [10], rectangular and D-sections [11,12] are some examples. However, to the author's best knowledge, very few publications are available in the literature that discuss the role of the geometrical parameters of the bluff body on the performance of the vibration. This paper reports geometrical effects on FIV performance of pivoted cylinders. Six cross section shapes are compared in which the circular cylinder is checked with our experimental data for validating the numerical simulation.
The present paper is organized as follows. The case is described in Section 2 followed by details of the numerical method, domains, and boundary conditions. Verification and validations is reported in Section 4, and the results are presented in in Section 5. We make conclusive remarks in Section 6.

Case Description
The cases considered in this study are based on water channel tests performed by Arionfard and Nishi [5]. The channel length is 1 meter with a test section's dimension of 30 cm wide by 30 cm deep. For the numerical simulation, the submerged bluff body is defined as a cylindrical solid sub-domain which is pivoted at a specific distance l, enabling rotation around the Z-axis, where X is the streamwise coordinate and Y is the cross-stream coordinate. The variation of the arm length l is considered by using twelve different values from −3D to +3D, where D is the diameter of the cylinder, negative values of l represent a pivot on the downstream of the cylinder (like Figure 1) and positive values represent a pivot point on the upstream side of the cylinder. A torsional spring is defined at the pivot point shown in Figure 1 and provides a restoring moment during oscillation.   Table 1.
In all simulations, Reynolds number is calculated based on the uniform inlet velocity (0.01 m/s) and the vertical height of the cross section (D ≈ 10 cm).

Governing Equations
The unsteady flow field around the cylinder is numerically simulated by employing 3D Unsteady Reynolds-Averaged Navier-Stokes equations (URANS) in Cartesian coordinates. Although the average Reynolds number in this study is low (≈1000), the local Reynolds number increases in near wall boundaries. Therefore, a turbulent model is necessary to model the behavior of the vortexes on the wake side. There are many turbulence models available and the choice of model depends on many factors such as the physics of the problem, the accuracy required and the computational power available. The K-Omega-SST model is used in this study because it is favored for predicting the formation of the vortices and flow separation [13,14]. By applying Reynolds decomposition and taking the time-average of continuity and momentum equations yields the URANS equations for incompressible flows [15].
The equation of motion for a rigid body in the polar coordinate with linear torsional spring and damper is expressed as where the dot symbol stands for differentiation with respect to time t, I t is the moment of inertia of the moving cylinder and C t is the total damping coefficient consist of the structural damping and equivalent generator damping. K represents the torsional stiffness of the spring and θ is the rotational displacement. M h f is the hydrodynamic angular momentum applied on the cylinder about the CG (center of gravity) of the cylinder given by where F p and F v are normal pressure and tangential viscous contributions. r p and r v are corresponding arm lengths from the center of the oscillating body (CG) to the center of rotation defined as where ρ is the density, s f ,i the face area vector, p the pressure, µ the dynamic viscosity, and R DEV the deviatoric stress tensor. The hydrodynamic force coefficients are calculated by using the built-in (forcecoeffs) function in OpenFoam given by where p Dyn = 1 2 ρAU 2 , A is the cross section area (84.35 cm 2 ), l is the arm length, and lift and drag forces are calculated from the vertical and horizontal components of F p and F v given by Equation (3). The structural parameters used in this study are described in Table 2 (solid domain).
Equation (1) and URANS equations are strongly coupled with the following steps: First, based on initial and boundary conditions the pressure distribution is calculated. Second, the forces on the cylinder surface corresponding to the pressure are calculated. Third, the equation of motion is solved based on the acquired forces and the displacements are calculated, and finally the domain is re-meshed according to the new position of the cylinder. The algorithm used by solvers is discussed in more details in the following. Step size (s) 1 × 10 −5 * * Automatically adjusted during the simulation base on the Courant number.

CFD Solver
The finite-volume-based open-source computational fluid dynamics library Open-FOAM is used to perform the numerical simulation of the flow field around the cylinder and solving the equation of motion. The governing equations were integrated over each control volume and the discrete values of the relevant quantities were determined at the center of the control volume. The diffusion term in the governing equations is discretized using second order central differencing scheme and for advection term, a second-order upwind scheme is utilized. To obtain a good resolution in time, time integration is performed by a second-order implicit scheme. Due to the unsteady nature of FIV, a PimpleDyM-Foam solver is used, which is a transient solver for turbulent incompressible flow on a moving mesh utilizing the PIMPLE (merged PISO-SIMPLE) algorithm. This solver is a modification of the pimpleFoam solver that supports meshes of class dynamicFvMesh. This class is a base class for meshes that can move and/or change topology. The built-in sixDoFRigidBodyMotion solver is utilized in the present study to model the rigid-body motion of the cylinder. One advantage of the sixDoFRigidBodyMotion is that the zone of dynamic mesh can be controlled with input parameters innerDistance and outerDistance, thus it is possible to fix the mesh near the cylinder wall. The fixed mesh moving with the cylinder ensures the large dynamical motion and computational accuracy of the flow near the cylinder wall. Otherwise, the finer mesh near the cylinder is vulnerable to be seriously distorted during the motion of the rigid body if the mesh near the cylinder wall is allowed to deform. Moreover, the fixed zone guarantees the accuracy of the outside boundary condition during the simulation.

Domain and Boundary Conditions
The mesh generation is performed by using the blockMesh and snappyHexMesh applications within the OpenFOAM package. A base hexahedral mesh is generated using blockMesh as a computational domain and the cylinder is snapped off the base mesh by using snappyHexMesh applications. Then, the remaining mesh is extruded to generate a 3D mesh.
The boundary condition on the cylinder is set to be a moving-wall, with no flux normal to the wall. The inlet boundary is defined as a velocity inlet with a uniform velocity of 0.01 m/s and zero pressure gradient was employed for the outlet. The top and bottom conditions defined as slip boundary while a no-slip condition is applied on the surfaces of the cylinder. The front and back walls are set to empty condition to simplify the simulation.
The initial conditions for the turbulence model were calculated from the inlet velocity and turbulence intensity at the inlet of the actual water channel, which was estimated by using PIV method. A summary of initial conditions is shown in Table 2.

Grid Independency and Validation
To reduce the computational cost and prevent mesh dependency, a preliminary study on necessary but sufficient resolution and domain size is done. To determine the domain size, six cases with different lengths and widths are simulated based on the CIR-3D conditions (CIR shape pivoted on the downstream with l = −3D). Then, the smallest size at which no further change is seen was selected. Similarly, the resolution of the background mesh (without refinement) is increased until the result did not change with increasing the mesh resolution. The most computationally efficient case is chosen based on the variation of C l , C d , and C m . According to the results of the domain size and resolution study shown in Figures 3 and 4, a refined domain size of 4D by 30D, with a resolution of 7680 elements and total domain size of 8D by 18D is chosen which leads to a blockage ratio of 0.125. Being aware of the limitations of this numerical study, we anticipate that the blockage potentially effects the sections in a similar way allowing comparison based on the difference in motion and hydrodynamic forces. An example of the mesh is shown in Figure 5.
The numerical model used in this study has been validated against our previous experimental results of a pivoted circular cylinder described in [5,6]. In the actual experiment, the cylinder is pivoted at a distance by using a connector arm and the Reynolds number is in the range of 2880 ≤ Re ≤ 22,300. A force moment sensor is used to measure the forces on the main shaft (at the pivot point) and then the measured forces and moments are used to calculate the hydrodynamic forces on the cylinder after dynamic and static tare. As the hydrodynamic forces are oscillating during the vibration, the corresponding amplitude to the peak frequency in the frequency domain is selected for evaluation after performing a Fast Fourier Transform (FFT). The numerical results are compared to the experiments done in the lowest Reynolds number in the experiment (Re = 2880). The numerical results are in good agreement with the experimental data according to Figure 6. Note that the experimental results are more accurate for Arm length ≈ 0 because for the smaller arm lengths the cylinder is more stationary and there is less turbulence induced noise on the cylinder as a result.

Results and Discussion
For simplicity, a Cartesian coordinate system is used for the discussion of results. The origin is the pivot point, and the X, Y, and Z axes are defined in the streamwise, transverse, and vertical direction, respectively. Figure 7 shows the maximum calculated power and the average amplitude of oscillation for each cross section including the corresponding arm ratio. The power spectrum density (PSD) of the angular velocity (which is widely being used for measuring the performance of random vibration converters) is used to calculate the power. The cumulative spectral power (CSP) of the PSD is then calculated by integrating over all frequencies base on the Parseval's theorem and used to estimate the dissipated power [16,17] where Q = √ KI t /C t is the quality factor and the PSD is calculated by using the fast Fourier transform of the angular velocity: whereθ(t) is the angular velocity of the vibration. By comparing the two charts, it is clear that the amplitude is not a proper performance metric even though it's been reported in many studies. For example, the highest power is achieved for the circular cross section while the box cross section oscillated with higher amplitude. According to the results, the angular velocity is lower near the ends of oscillation for the non-circular cross sections when pivoted on the downstream. There are two possible reasons for lower angular velocity in a cross section: First is the higher drag force in a higher angle of attacks in non-circular cross sections [18]. Higher drag force changes the stiffness nonlinearly and shifts the natural frequency f v out of the lock-in range based on the following equation derived from the equation of motion: where +l and −l correspond to the location of the pivot point on the upstream side or downstream side, respectively. A D0 = 1 2 ρDH w C D U 2 and DH w is the projected area of the cross section. For higher arm length, this increase in drag completely suppresses the vibration. The second reason is the lower spanwise correlation length. The first reason is discussed as vibration mechanism followed by a discussion over vorticity dynamic to understand the behavior of fluid around each section and its effect on correlation length.

Vibration Analysis
According to Figure 7, only the CIR and EGG cross sections produce reasonable power followed by BOX with a large difference (of around 50% lower power). The drag, lift, and moment ratio versus arm ratio are shown in Figure 8 for all sections along with the calculated power on a separate axis. These three sections show more power while the pivot is on the downstream side of the section as shown in Figure 8. However, the power is nearly zero for a BOX section pivoted on the upstream side (l > 0) regardless of the length of the arm. The opposite behavior is observed for DIA section: the power is nearly zero for a DIA section pivoted on the downstream side (l < 0) regardless of the length of the arm. The two remaining sections (RAU, TRI) show the lowest power with almost no effect of the arm length and the pivot location. The difference between the BOX and the rest of the sections is more clear by analyzing the vibration response shown in Figure 9.
The vibration frequency ( f v ) is far away from the natural frequency while the pivot is at the downstream of the section but it gradually goes up and close to the natural frequency. Even though aeroelastic instability is expected to be responsible for oscillation in this kind of cross section, the lock-in phenomena seem to improve the oscillation for sections with round edges. A similar change is seen for the Strouhal number St (= f s D/U), where = f s is the predominant vortex shedding frequency), as shown in Figure 9. The Strouhal number is very low for the BOX section while pivoted on the downstream. It eventually increases by the arm length and converges to 0.13 but for the rest of the sections, the Strouhal number is close to 0.2 which is considered in the lock-in range.
The maximum power depends largely on the natural frequency of the system which is a function of the pivot location and spring stiffness in our setup. Arionfard and Nishi [5] found that for a circular cross section the drag force assists the motion by reducing the natural frequency when the pivot is located at the downstream side of the cylinder (l < 0) based on Equation (7). As the moment of inertia (I t ), flow velocity (U), spring stiffness (K), and the projected area of the sections are constant, l A D0 ∝ l C D is responsible for changes in the natural frequency.
Note that the mathematical analysis provided in [5] is only valid for round shapes where the drag and lift coefficients are not a function of the angle of attack. This is with agreement with the results shown in Figure 8: The calculated power changes with l C D for CIR, EGG, and RAU shape while the calculated power for BOX, DIA, and TRI shapes shows less dependency to the drag coefficient or arm length.

Vorticity Analysis
The steady-state vorticity field for the cases with the highest performance is shown in Figure 10. For sufficient oscillation amplitudes, symmetrical shedding with 2S mode is triggered in all cases as expected due to the low Reynolds number. The 2S mode is associated with the initial branch [19] where two single vortices shed per cycle, one by the top shear layer and another one by the bottom shear layer. The vorticity field animations can be found in the Supplementary Videos S1-S6.
To compare the correlation length (which is a measure of the span-wise length, that the vortices remain in phase) for each section, the three-dimensional state of the wake for each simulated case is visualized in Figures 11 and 12. The wakes are extracted by using a threshold filter the way that the pressure lies within 10 to 100 Pascal for all cases. A few factors influence the correlation length in FIV, including the amplitude of vibration, aspect ratio, surface roughness and the Reynolds number [15]. Here, the Reynolds number and surface roughness are similar for all cases while the amplitude of vibration and aspect ratio (which is a function of geometry) are changing.
The correlation length is higher when the pivot is at the upstream side (l > 0) for all sections except for the BOX. It is well known that body motion at a frequency close to that of the natural vortex shedding has a strong organizing effect on the shedding wake, which is manifested by a sharp increase in the spanwise correlation of the flow and forces on the body. However, the increase in three dimensionality of the flow behind the BOX section contradicts this pattern. A similar increase in three dimensionality is observed for the TRI section as well, but it is due to smaller vibration amplitude for all lengths in this section. The formation of the vortex line for the cross sections with the highest calculated power is more evident in the animations provided in the online Supplementary Videos V1-V6.

Conclusions
3D numerical simulation of fluid-induced vibration has been reported for a series of cylinders with different cross sections including circular, rectangle, diamond, triangle, reuleaux, and egg shape. The cylinders are pivoted at distance from the centre to study the geometrical effect of the FIV performance and to compare the results with our previous experimental study. The cross-sectional area, moment of inertia, spring stiffness, inlet velocity, and damping coefficient are set to be similar for all cases to eliminate the effect of non-geometrical parameters. According to the results, the circular and egg shape cross sections are the most efficient shapes regardless of the pivot location followed by the box, diamond, reuleaux, and triangle shapes. The vorticity field shows that the 2S mode is triggered for all cases mainly due to the low Reynolds number; thus, the vibrations are expected to be in the initial branch. Moreover, 3D visualization of the wake for each section shows that the correlation length is higher for round shapes especially when the pivot is at the upstream side while for the shapes with sharp edges, the three-dimensionality of the wake is higher.
There are two major limitations in this study that could be addressed in future research. First, the domain size: even though a grid independency study is done for the circular cylinder and there is a good agreement with the experiment, similar results are not necessarily expected for other cross sections or arm lengths. This applies to the blockage ratio as well. It is assumed that the blockage has a similar effect on all cases if kept constant for all cross sections. Second, the Reynolds number: the results are compared to the experiments done with Reynolds number of around 2800 assuming both numerical and experimental tests are in the same flow regime (1000 ≥ Reynolds ≥ 3000). Moreover, the Reynolds number in this study is much smaller than that of actual operating conditions. Being aware of the limitations of this numerical study, we concluded that the hydrodynamic forces, displacement and calculated power of the cross sections are still comparable with each other if not to the experiment.

Supplementary Materials:
The following are available online at https://www.mdpi.com/1996-107 3/14/4/1128/s1, Video V1: The velocity field on the wake side of the BOX (l = −2D), Video V2: The velocity field on the wake side of the CIR (l = −3D), Video V3: The velocity field on the wake side of the DIA (l = +2.5D), Video V4: The velocity field on the wake side of the EGG (l = −2.5D), Video V5: The velocity field on the wake side of the RAU (l = +3D), Video V6: The velocity field on the wake side of the TRI (l = +2.5D), Video S1: The vorticity field for the CIR (l = −3D), Video S2: The vorticity field for the DIA (l = +2.5D), Video S3: The vorticity field for the EGG (l = −2.5D), Video S4: The vorticity field for the RAU (l = +3D), Video S5: The vorticity field for the TRI (l = +2.5D), Video S6: The vorticity field for the BOX (l = −2D).