**2. 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.

**Figure 1.** Schematic diagram of a typical computational domain and boundaries. Here, the pivot point is located at the downstream of the cylinder. The solid sub-domain is Ω*S*, the fluid sub-domain is Ω*F*, and Ω*RF* is the refined part of the fluid sub-domain. The arm length *l* is the distance between the center of the bluff body and the pivot point

Six cross section shapes are chosen in this study as solid sub-domains. The area of all cross sections are equal and the height of sections (D) is as similar as possible to keep the Reynolds number within the same range in all simulations. More details of the geometrical parameters of each cross-section is shown in Figure 2 and described in Table 1.

**Figure 2.** Geometry of the bluff bodies. Dimensions are given in 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).

**Table 1.** Geometric parameters of the cross sections.


#### **3. Numerical Method**

#### *3.1. 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

$$M\_t \ddot{\theta} + \mathbb{C}\_t \dot{\theta} + K\theta = M\_{hf} \tag{1}$$

where the dot symbol stands for differentiation with respect to time *t*, *It* is the moment of inertia of the moving cylinder and *Ct* 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. *Mh f* is the hydrodynamic angular momentum applied on the cylinder about the CG (center of gravity) of the cylinder given by

$$M\_{hf} = F\_{p}r\_{p} + F\_{\upsilon}r\_{\upsilon} \tag{2}$$

where *Fp* and *Fv* are normal pressure and tangential viscous contributions. *rp* and *rv* are corresponding arm lengths from the center of the oscillating body (CG) to the center of rotation defined as

$$\begin{aligned} F\_p &= \sum\_i \rho\_i s\_{f,i} (p\_i - p\_{ref}) \\ F\_v &= \sum\_i s\_{f,i} (\mu R\_{DEV}) \end{aligned} \tag{3}$$

where *ρ* is the density, *sf* ,*<sup>i</sup>* the face area vector, p the pressure, *μ* the dynamic viscosity, and *RDEV* the deviatoric stress tensor. The hydrodynamic force coefficients are calculated by using the built-in (forcecoeffs) function in OpenFoam given by

$$\begin{aligned} \mathbf{C}\_{l} &= \frac{liftForce}{p\_{Dyn}}\\ \mathbf{C}\_{d} &= \frac{dragForce}{p\_{Dyn}}\\ \mathbf{C}\_{m} &= \frac{M\_{lf}}{p\_{Dyn}l} \end{aligned} \tag{4}$$

where *pDyn* = <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>AU*2, *<sup>A</sup>* is the cross section area (84.35 cm2), *<sup>l</sup>* is the arm length, and lift and drag forces are calculated from the vertical and horizontal components of *Fp* and *Fv* 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.


**Table 2.** Initial conditions for solid and fluid domain.

\* Automatically adjusted during the simulation base on the Courant number.

### *3.2. 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.

#### *3.3. 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.

#### **4. 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* = −3*D*). 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 *Cl*, *Cd*, and *Cm*. 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.

**Figure 3.** The variation of hydrodynamic coefficients with the domain size.

**Figure 4.** The variation of hydrodynamic coefficients with the mesh size.

**Figure 5.** An example of the mesh with reuleaux shape snapped off of the grid.

**Figure 6.** Comparison between the numerical results and experiment.

#### **5. 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]

$$CSP = P = \frac{I\_t}{Q} \int\_0^\infty (PSD) df\tag{5}$$

where *<sup>Q</sup>* <sup>=</sup> <sup>√</sup>*K It*/*Ct* is the quality factor and the PSD is calculated by using the fast Fourier transform of the angular velocity:

$$PSD = \left| FFT(\dot{\theta}(t)) \right|^2 \tag{6}$$

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 *fv* out of the lock-in range based on the following equation derived from the equation of motion:

$$f\_N = \frac{1}{2\pi} \sqrt{\frac{K \pm I \, A\_{D0} l I^2}{I\_t}} \tag{7}$$

where +*l* and −*l* correspond to the location of the pivot point on the upstream side or downstream side, respectively. *AD*<sup>0</sup> = <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>DHwCDU*<sup>2</sup> and *DHw* 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.

**Figure 7.** The maximum power output (**a**) and the mean amplitude (**b**) for each cross section.
