# Study on Fluid-Induced Vibration Power Harvesting of Square Columns under Different Attack Angles

Authors: Meng Zhang; Guifeng Zhao; Junlei Wang
Journal: Geofluids
 (2017)

---

## Abstract

A model of the flow-vibration-electrical circuit multiphysical coupling system for solving square column vortex-induced vibration piezoelectric energy harvesting (VIVPEH) is proposed in this paper. The quasi steady state theory is adopted to describe the fluid solid coupling process of vortex-induced vibration based on the finite volume method coupled Gauss equation. The vibrational response and the quasi steady state form of the output voltage are solved by means of the matrix coefficient method and interactive computing. The results show that attack angles play an important role in the performance of square column VIVPEH, of whichα=45° is a relatively ideal attack angle of square column VIVPEH.


---

## Body

## 1. Introduction

Recently, the development and utilization of new energy sources have become a research hotspot, among which the study of capturing energy from environment has received much more attention. One of the most important ways of environmental energy harvesting is capturing energy from fluid, which can be divided into two kinds, wind energy and water energy. Most of the traditional wind power and hydroelectric power facilities use the rotating turbine device to harvest energy with large volume device and low energy density. The micro energy technology, which can extract energy from environment and convert it into electric energy [1–4], has the features of functional continuity, small volume, and high energy density. In the late 1990s, the vibration piezoelectric energy harvesting technology has been widely used to harvest environmental flow energy and convert it into vibration energy [5–9], which is a kind of micro energy technology with continuous and nonconsuming energy supply. Therefore, it is an effective method for the microminiaturization of flow induced vibration energy harvesting device.In fluid dynamics, there is a potential physical phenomenon that can be used for energy harvesting called vortex-induced vibration. Vortex-induced vibration is that when the fluid flows through the bluff body, the formation and periodic shedding of the vortex will cause the vibration of the bluff body. Once the vibration intensity reaches a certain level, the flow field shedding will be locked, which results in large vibration energy. In other words, vortex-induced vibration is a kind of periodic, steady, or unsteady fluid structure interaction phenomenon, which has the characteristics of continuity and easy excitation [10, 11].It is a challenging work to solve the problem of vortex-induced vibration energy harvesting. The key problem here is how to transform the flow energy into vibration energy efficiently. In recent years, many meaningful research works have been carried out on using vortex-induced vibration to collect ocean energy and wind energy. Among them, the energy conversion of circular bluff body piezoelectric vortex-induced vibration is mostly concerned. Allen and Smits [12] have studied the theory of energy harvesting of piezoelectric materials and designed an “eel” energy harvesting model which can be used to harvest the fluid kinetic energy in the water tank. On the basis of film theory, the “eel” model device can also be used to harvest the vortex shedding energy of “lock-in” phenomenon. Taylor et al. [13] used the “eel” device of a PVDF polymer to harvest marine energy, which was placed in a water tank with the length of 241.3 mm, the width of 76.2 mm, and the thickness of 150 μm. The research shows that when the flapping frequency of PVDF closes to the vortex shedding frequency, the energy collection performance will be improved; the maximum voltage 3 V appears in the water flow velocity of 0.5 m/s. In Marine Renewable Energy Laboratory (MRELab) at the University of Michigan, Bernitsas et al. [14, 15] have presented a vortex-induced vibration for aquatic clean energy (VIVACE) to utilize the VIV phenomenon to generate power. The latter studies have been conducted in support of model tests for VIVACE converter, which harnesses hydrokinetic energy enhancing flow induced motions (FIM) and particularly VIV and various forms of galloping. Lee and Bernitsas [16, 17] built a device/system Vck to replace the physical damper/springs of the VIVACE with virtual elements. The testing was performed in the Low Turbulence Free Surface Water Channel of the University of Michigan at 40000<Re<120000 and damping 0<ξ<0.16. To continue and improve the work of VIVACE, which is used to convert hydrokinetic energy from ocean/river currents to electricity. Raghavan and Bernitsas [18, 19] conducted experiments in the regime right before transition from laminar to turbulent flow. Hysteresis was not observed in the experiments because the parameters remain in the 2P domain. Higher VIV amplitude of 2.1 to 2.7 diameters was obtained and the range of synchronization was increased by increasing the damping in VIVACE. Chang et al. [20] put a PTC device in the surface of the VIVACE converter to enhance the viscous damping; the results show that proper PTC could enhance the conversion efficiency of hydrokinetic energy to mechanical in VIVACE. For solving the problem of VIV, some excellent researches have also been conducted in MRELab. With the help of OpenFOAM software; Wu et al. [21] have developed a CFD code to solve the VIV problem of a single cylinder with PTC. Ding et al. [22] have also developed a CFD code in OpenFOAM to solve the VIV problem of multiple circular cylinders. Abdelkefi et al. [23, 24] and Mehmood et al. [25] investigated the possibility of harvesting energy from elastically mounted circular cylinders. They used, respectively, a wake oscillator model and direct numerical simulations to determine the fluctuating lift force. Abdelkefi et al. [23] reported that the aerodynamic nonlinearity results in the presence of hardening behavior. Both research studies determined the effects of the load resistance on the synchronization region and level of the harvested power. They demonstrated that maximum level of the harvested power can be associated with minimum values in the transverse displacement due to the shunt damping effect. To investigate the possibility of harvesting energy from other flow induced vibrations, different investigations have been performed. Abdelkefi and Nuhait [26] investigated the effects of cambered wing-based piezoaeroelastic energy harvesters on the flutter speed and the performance of the harvester. Dai et al. [27] and Yan and Abdelkefi [28] investigated the potential of harvesting energy from a combination of base excitations and VIV or galloping, respectively. Robbins et al. [29] proposed a device similar to the “eel” model, consisting of piezoelectric film arrays (single size 203 × 279 × 0.5 mm) with different configurations to collect wind energy generated by induced draft fan; the maximum output power is up to 10 mW. Akaydin et al. [30] proposed a structure (main size 30 × 16 × 0.2 mm) which can separate the piezoelectric film from the bluff body, in which the maximum output power is 4 μW. Gao et al. [31] carried out an experiment model (unit piezoelectric element size is 58 × 10 mm, cylinder diameter is 29.1 mm, and cylinder length is 36 mm), in which the cylinder is connected with the PZT cantilever beam, and the PZT beam is fixed at one end in the wind tunnel. Thus, electrical charges were accompanied by the wind-induced vibration of the PZT cantilever beam. The “lock-in” phenomenon is observed in the wind-induced vibration and a high voltage output is obtained in this region. Akaydin et al. [32] presented an energy harvesting experiment in a wind tunnel, in which the cylinder is connected with the free end of a cantilever beam with aluminum shim device. The results show that when the circuit resistance is changed, the natural frequency of the system will be changed also, so that the maximum value of the electromechanical coupling damping can be obtained. The above researches are mainly focused on the energy collection models of the cylindrical bluff body under vortex-induced vibration. Studies on the vortex-induced vibration and piezoelectric energy harvesting of noncylindrical section, such as square column and triangular column, are relatively less.It is very important to consider the different attack angles in the analysis of square column vortex-induced vibration because of the asymmetry. The motion of the fluid at the edge of the square column is very irregular, which will greatly affect the vortex shedding near the column wall. The square column vortex shedding process is different from that of smooth circular cylinder; there is a direct vortex shedding without a smooth transition region. The boundary layer separation point of the square column is fixed at the corner of the leading edge without changing its position with the varying of Reynolds number. Therefore, it is of great practical significance to study on the influence of different attack angles of square column vortex-induced vibration piezoelectric energy harvesting (VIVPEH).In this paper, a model of square column VIVPEH is presented, which can be used to investigate the mechanism of different attack angles on the energy capturing efficiency. The quasi steady state theory is adopted to describe the fluid solid coupling process of vortex-induced vibration based on the finite volume method coupled Gauss equation. Finally, the performance of different attack angles on the VIVPEH is carried out with considering the flow field parameters, the vibration characteristics, and the energy collection parameters.


## 2. Physical Model

The physical model of square column VIVPEH is shown in Figure1, in which U stands for the inflow velocity, α stands for the angle between the square column and the inflow velocity, M is the mass of the column, K is the elastic coefficient of the system, C is the damping of the system, and R is the external resistance load.Figure 1
Physical model of square column energy harvesting system.
(a)
Side view of the system
(b)
Mass-Spring-Damping-Resistance load system


## 3. Mathematical Model of Energy Harvesting

To analyze the coupling process of three different fields, external flow field, vortex-induced vibration, and circuit, this paper uses theNavier-Stokes equation to describe the vortex-induced vibration, uses a linear second-order differential equation to describe the vortex-induced vibration of the single degree of freedom M-C-K (mass spring damping) system, and finally uses coupledGauss law and vibration equations to describe the electromechanical coupling system.


### 3.1. Fluid Solid Coupling Model

The external flow field is calculated by the continuity equation and theNavier-Stokes equation. Flow simulations presented in this paper are produced by open source CFD tool OpenFOAM, which is composed of C++ libraries solving continuum mechanics problems with a finite volume discretization method. Suppose the external flow field is 2D and unsteady. The time-dependent viscous flow solutions can be obtained by numerical approximation of the incompressible unsteady Reynolds-AveragedNavier-Stokes (URANS) equations in conjunction with the one-equationSpalart-Allmaras(S-A) turbulence model [33], where a second-orderGauss integration scheme with a linear interpolation is used in the governing equations for the divergence, gradient, and Laplacian terms. For time integration, second-order backwardsEuler method is employed. The numerical discretization scheme has second-order accuracy in space and time. A pressure implicit with splitting of operators (PISO) algorithm is used for solving momentum and continuity equations together in a segregated way. The equations of motion for the square column are solved using a second-order mixed implicit and explicit time integration scheme. The basic equations of URANS are(1)∂Ui∂xi=0,∂Ui∂t+Uj∂Ui∂xj=-1ρ∂p∂xi+∂∂xj2νSij-uj′ui′¯,where p is pressure, ρ is fluid density, ν is dynamic viscosity, Ui is the mean flow velocity vector, and Sij is the strain rate tensor,(2)Sij=12∂Ui∂xj+∂Uj∂xi.To solve the URANS equations for mean flow properties and potential turbulence flow, the Boussinesq eddy-viscosity approximation is adopted here, which relates to the Reynolds stress and the velocity gradient. The quantityρuj′ui′¯ is the Reynolds stress tensor and can be modeled as ρuj′ui′¯=2μtSij, where μt is the turbulence eddy viscosity.The S-A model is widely used for turbulence closure. The eddy-viscosity coefficient can be calculated from the following transport equation:(3)∂ν~∂t+uj∂ν~∂χj=cb1S~ν~-cw1fwν~d2+1σ∂∂χjν+ν~∂ν~∂χj+cb2∂ν~∂χi·∂ν~∂χj.In S-A model, the turbulence eddy viscosityμt can be obtained by(4)μt=ρν~fv1,in which fv1=χ3/χ3+c3ν1, and χ≡ν~/ν, where χ is an intermediate value. ν~ is working variable of the turbulence model and depends on the transport equation (3).The details of the transport equation are in Spalart and Allmaras [33], and the trip terms ft1 and ft2 are switched off and a “trip-less” initial condition was added for solving working variable ν~. This approach was successfully used in the work of Wu et al. [21] and Ding et al. [22] for circular cylinders.The motion equations of a single degree of freedomM-C-K system can be represented by the following linear second-order differential equation:(5)MY¨+CY˙+KY=Fy.The relationships betweenM, K, and C are as follows:(6)ωn=KM,C=2ξMωn,where Fy stands for the force on unit volume of the flow field, which is perpendicular to the flow direction,Y stands for the column vibration displacement, Y˙ and Y¨ represent the first- or second-order derivative of the column vibration displacement, respectively, ωn is the natural circular frequency, and ξ is the dimensionless damping ratio. The dynamic response of the column can be obtained by solving the motion equations and the fluid governing equations simultaneously.


### 3.2. Electromechanical Coupling Model

In order to describe the relationship between the amplitude and the voltage in the vortex-induced vibration circuit,Gauss lawis adopted in this paper. Theoretical derivations are as follows:(7)MY¨+CY˙+KY-θV=Fy,(8)θY˙+CpV˙+VR=0,where θ is the electromechanical coupling coefficient, Cp is the capacitance coefficient, and V is the voltage.The influences of the vortex-induced vibration system on the circuit output voltage have been considered in (7) and (8), respectively. At the same time, the negative feedback effect of the circuit on the vibration system is also taken into account; that is to say, the influence of electromechanical coupling is considered. Combined with the flow field calculation results, we can carry out the flow-mechanical-electrical coupling analysis.In order to solve the damping and natural frequency of the system, we use the matrix method to calculate the two-order nonhomogeneous ordinary differential equation (7). The homogeneous equation of (7) is as follows:(9)MY¨+CY˙+KY-θV=0.Let X1=Y, X2=Y˙, and X3=V; substitute (6) into (8) and (9); we can get (10)X˙1=X2X˙2=-ωn2X1-2ξωnX2+θMX3X˙3=-θCpX2-X3RCp.The above equations can be expressed in the following matrix form:(11)X˙=BRX,where(12)X=X1,X2,X3T,BR=010-ωn2-2ξωnθM0-θCp-1RCp.The matrixB(R) has three different eigenvalues of ki, of which i=1,2,3. Spalart and Allmaras [33] have pointed out that the first two eigenvalues are similar to the ones of vibration system without circuit, and yet the third eigenvalue is associated with the electromechanical coupling effect, such as piezoelectric system affected by the foundation or the aeroelastic excitation, and is negative constant. There are conjugate relations between k1 and k2, in which the real part and the imaginary part of the conjugate solution stand for the damping and natural frequency of the electromechanical coupling system, respectively. Given that k3 is negative constant, we consider only the real part of k1 and k2, when computing the trivial solution of the matrix B(R).


## 3.1. Fluid Solid Coupling Model

The external flow field is calculated by the continuity equation and theNavier-Stokes equation. Flow simulations presented in this paper are produced by open source CFD tool OpenFOAM, which is composed of C++ libraries solving continuum mechanics problems with a finite volume discretization method. Suppose the external flow field is 2D and unsteady. The time-dependent viscous flow solutions can be obtained by numerical approximation of the incompressible unsteady Reynolds-AveragedNavier-Stokes (URANS) equations in conjunction with the one-equationSpalart-Allmaras(S-A) turbulence model [33], where a second-orderGauss integration scheme with a linear interpolation is used in the governing equations for the divergence, gradient, and Laplacian terms. For time integration, second-order backwardsEuler method is employed. The numerical discretization scheme has second-order accuracy in space and time. A pressure implicit with splitting of operators (PISO) algorithm is used for solving momentum and continuity equations together in a segregated way. The equations of motion for the square column are solved using a second-order mixed implicit and explicit time integration scheme. The basic equations of URANS are(1)∂Ui∂xi=0,∂Ui∂t+Uj∂Ui∂xj=-1ρ∂p∂xi+∂∂xj2νSij-uj′ui′¯,where p is pressure, ρ is fluid density, ν is dynamic viscosity, Ui is the mean flow velocity vector, and Sij is the strain rate tensor,(2)Sij=12∂Ui∂xj+∂Uj∂xi.To solve the URANS equations for mean flow properties and potential turbulence flow, the Boussinesq eddy-viscosity approximation is adopted here, which relates to the Reynolds stress and the velocity gradient. The quantityρuj′ui′¯ is the Reynolds stress tensor and can be modeled as ρuj′ui′¯=2μtSij, where μt is the turbulence eddy viscosity.The S-A model is widely used for turbulence closure. The eddy-viscosity coefficient can be calculated from the following transport equation:(3)∂ν~∂t+uj∂ν~∂χj=cb1S~ν~-cw1fwν~d2+1σ∂∂χjν+ν~∂ν~∂χj+cb2∂ν~∂χi·∂ν~∂χj.In S-A model, the turbulence eddy viscosityμt can be obtained by(4)μt=ρν~fv1,in which fv1=χ3/χ3+c3ν1, and χ≡ν~/ν, where χ is an intermediate value. ν~ is working variable of the turbulence model and depends on the transport equation (3).The details of the transport equation are in Spalart and Allmaras [33], and the trip terms ft1 and ft2 are switched off and a “trip-less” initial condition was added for solving working variable ν~. This approach was successfully used in the work of Wu et al. [21] and Ding et al. [22] for circular cylinders.The motion equations of a single degree of freedomM-C-K system can be represented by the following linear second-order differential equation:(5)MY¨+CY˙+KY=Fy.The relationships betweenM, K, and C are as follows:(6)ωn=KM,C=2ξMωn,where Fy stands for the force on unit volume of the flow field, which is perpendicular to the flow direction,Y stands for the column vibration displacement, Y˙ and Y¨ represent the first- or second-order derivative of the column vibration displacement, respectively, ωn is the natural circular frequency, and ξ is the dimensionless damping ratio. The dynamic response of the column can be obtained by solving the motion equations and the fluid governing equations simultaneously.


## 3.2. Electromechanical Coupling Model

In order to describe the relationship between the amplitude and the voltage in the vortex-induced vibration circuit,Gauss lawis adopted in this paper. Theoretical derivations are as follows:(7)MY¨+CY˙+KY-θV=Fy,(8)θY˙+CpV˙+VR=0,where θ is the electromechanical coupling coefficient, Cp is the capacitance coefficient, and V is the voltage.The influences of the vortex-induced vibration system on the circuit output voltage have been considered in (7) and (8), respectively. At the same time, the negative feedback effect of the circuit on the vibration system is also taken into account; that is to say, the influence of electromechanical coupling is considered. Combined with the flow field calculation results, we can carry out the flow-mechanical-electrical coupling analysis.In order to solve the damping and natural frequency of the system, we use the matrix method to calculate the two-order nonhomogeneous ordinary differential equation (7). The homogeneous equation of (7) is as follows:(9)MY¨+CY˙+KY-θV=0.Let X1=Y, X2=Y˙, and X3=V; substitute (6) into (8) and (9); we can get (10)X˙1=X2X˙2=-ωn2X1-2ξωnX2+θMX3X˙3=-θCpX2-X3RCp.The above equations can be expressed in the following matrix form:(11)X˙=BRX,where(12)X=X1,X2,X3T,BR=010-ωn2-2ξωnθM0-θCp-1RCp.The matrixB(R) has three different eigenvalues of ki, of which i=1,2,3. Spalart and Allmaras [33] have pointed out that the first two eigenvalues are similar to the ones of vibration system without circuit, and yet the third eigenvalue is associated with the electromechanical coupling effect, such as piezoelectric system affected by the foundation or the aeroelastic excitation, and is negative constant. There are conjugate relations between k1 and k2, in which the real part and the imaginary part of the conjugate solution stand for the damping and natural frequency of the electromechanical coupling system, respectively. Given that k3 is negative constant, we consider only the real part of k1 and k2, when computing the trivial solution of the matrix B(R).


## 4. Quasi Steady State Model for Output Voltage

In this section, we adopt the quasi steady state model proposed by Barrero-Gil et al. [34] to describe the amplitude of vortex-induced vibration and calculate the time-varying vibration energy harvesting. When the vortex-induced vibration is in the synchronization region, the vibration amplitude of the system can be expressed as the following sine function:(13)Y=Ymaxsin⁡ωnt,where Ymax is the maximum column vibration displacement.It is worth noting that the voltage time-history curve and the vibration amplitude time-history curve are synchronous; that is to say, there is no phase difference. Substituting (13) into (11), we can get the analytical solution of the quasi steady state voltage with MATLAB.(14)Vt=θωnRYmax1+ωn2R2Cp2e-t/RCp-cos⁡ωnt-ωnRCpsin⁡ωnt.The corresponding output power can be obtained by the following equation:(15)Pt=V2tR.The above mathematical expressions include the solving process of fluid solid coupling and electromechanical coupling.First of all, we can get the flow field pressuresP and Fy by solving the URANS equations in OpenFOAM. Secondly, we can obtain the column vibration displacement Y, the damping C, and natural circular frequency ωn, by solving (5) to (11). It should be noted that the motion of the column is influenced by the pressure of the flow field. At the same time, the vibration of the column gives feedback to the flow field and causes the change of flow field distribution. Thus, the fluid solid coupling problem can be solved by interactive computing. Finally, the time-history curve of output voltage and output power can be obtained by (14) and (15). The above is the whole computing process of flow-vibration-electromechanical coupling system.


## 5. Numerical Results

### 5.1. Analysis Parameters and Cases

To carry out the numerical calculation of vortex-induced vibration, this paper provides the parameters of the vibration energy harvesting system, as shown in Table1, in which Dsqu is the side length of the square column, DNor stands for the dimensionless characteristic length of the square column, and the relationship between DNor and Dsqu is given in(16)DNor=Dsqusin⁡α+cos⁡α.Based on the parameters of Table 1 and formula (16), we can calculate the numerical values of DNor as 0.0016 under 0° attack angle, 0.00196 under 15° attack angle, 0.00219 under 30° attack angle, 0.00226 under 45° attack angle, 0.00219 under 60° attack angle, and 0.00196 under 75° attack angle.Table 1
Calculation parameters for square column VIVPEH system.
Symbols and units
Physical meanings
Numerical values
M [kg]
Square column mass
0.2979
K [N/m]
Elastic coefficient of the system
579~584
C [N·s/m]
Damping of the system
0.0325~0.45
D
s
q
u [m]
Square column side length
0.0016
α [degree]
Angle between the square column and the inflow velocity
0~75
ζ
Damping ratio of the system
0.00121
f
n [Hz]
Natural frequency of the system
7.012~7.132 (water)
C
p [nF]
Capacitance
120
θ
Electromechanical coupling coefficient
1.55 × 10−3
μ [Pa·s]
Dynamic viscosity
0.0011379
ν [m2/s]
Kinematic viscosity
1.139 × 10−6 (water)
ρ [kg/m3]
Density
999.1026 (water)The computational grid and boundary conditions of square column vortex-induced vibration under different attack angles are given in Figure2, where the computational domain is 20 × 20D, and the entire domain includes five boundaries: velocity inlet, velocity outlet, top, bottom, and a column wall. The inlet velocity is considered as uniform and constant velocity. For outlet boundary, a zero gradient condition is specified for velocity. The top and bottom condition are defined as a wall boundary. In present numerical study, a moving wall boundary condition is applied for the square column when the column is in VIV. The two-dimensional, structured grids were generated with the help of “Gambit” software. The grid domain size is 20 × 20D. The square column was set in the center in the domain to ensure that the results of the numerical model are accurate. The conditions at the outlet are close to the assumed conditions. The computational domain in the vicinity of each cylinder is a 3 × 3D square where the grid density for the near-wall region is enhanced to solve for high resolution in flow properties.Figure 2
Computational grid and boundary conditions.
(a)
α
=
0°
(b)
α
=
15°
(c)
α
=
30°
(d)
α
=
45°
(e)
α
=
60°
(f)
α
=
75°
(g)When the external resistance value isR = 1 × 106 Ω, the spring stiffness value of the square column VIVPEH system is Ksqu=580 N/m, and the corresponding damping value is Csqu=0.2 Ns/m. Taking the range values of the flow velocity as 0.03927 m/s to 0.08975 m/s, we can calculate the numerical values of Ursqu under different attack angles, as shown in Table 2, in which Ursqu is the reduction velocity, defined as follows:(17)Ursqu=UfnDsqu,where fn=ωn/2π.Table 2
Run case of square column in OpenFOAM.
U (m/s)
Ursqu
α
=
0°
α
=
15
° (α=75°)
α
=
30
° (α=60°)
α=45°
0.03927
3.5
2.85784
2.56223
2.47525
0.04488
4
3.26611
2.92826
2.82885
0.05049
4.5
3.67437
3.29429
3.18246
0.0561
5
4.08263
3.66032
3.53607
0.05834
5.2
4.24594
3.80673
3.67751
0.06059
5.4
4.40924
3.95315
3.81895
0.06283
5.6
4.57255
4.09956
3.9604
0.06732
6
4.89916
4.39239
4.24328
0.07293
6.5
5.30742
4.75842
4.59689
0.07854
7
5.71569
5.12445
4.9505
0.08415
7.5
6.12395
5.49048
5.3041
0.08976
8
6.53221
5.85652
5.65771
0.09537
8.5
6.94048
6.22255
6.01132In this paper, the matrix method is used to evaluate the influence of the external resistance load on the damping and natural frequency of the electromechanical coupling system by MATLAB. Then, the damping and natural frequency values can be used for initial conditions in OpenFOAM to compute the vibration amplitude of the vortex-induced vibration system under different Reynolds numbers(94<Re<115). The main focus of this paper is to investigate the influence of external resistance load (R = 1 × 103 Ω, 1 × 104 Ω, 1 × 105 Ω, 1 × 106 Ω, 1 × 107 Ω) on the energy conversion of the square column VIVPEH system, including the effect on vibration amplitude, output voltage, and power.


### 5.2. System Damping and Natural Frequency Characteristics

The damping and natural frequency of the electromechanical coupling system can be obtained by MATLAB, and the real and imaginary parts of the circuit conjugate solutions are shown in Figure3.Figure 3
Real and imaginary parts of the circuit conjugate solutions.According to Figure3, we can see that the total damping of the system is small when the resistance load is small. In the case of R<1×105 Ω, the system total damping is increased with the resistance load increasing. Once the resistance load R reaches 1 × 105 Ω, the system total damping reaches the maximum value. It is noteworthy that the total damping of the system decreases instead of increasing when the resistance load keeps increasing; that is, R>1×105 Ω. For the natural circular frequency, when R<3×104 Ω, the frequency value remains at 44 rad/s; when R>2×106 Ω, the frequency value approximately remains at 50 rad/s. Generally, the natural circular frequency value of the system is relatively stable, which is kept in the range of 44–50 rad/s.


### 5.3. Vibration Characteristics of Square Column VIVPEH under Different Attack Angles

In the following numerical simulations, we take the resistance loadR = 1 × 106 Ω and compute the amplitude response of the square column VIVPEH under different attack angles and different flow velocities, as shown in Figure 4. Note here that the natural frequency of square column is fn=6.98 Hz, when the resistance load is R = 1 × 106 Ω. The dimensionless vibration amplitude Ysqu/Dsqu is adopted to indicate the vibration response of the square column. The maximum value of Ysqu can be obtained by means of averaging the peak value of at least 60 displacement time-history response curves. In this paper, the vibration amplitude of the smooth circular column is provided to compare with that of the square column.Figure 4
Dimensionless vibration amplitude of square column VIVPEH under different attack angles and different flow velocities.As we can see in Figure4, different attack angles obviously have an effect on the peak vibration amplitude and lock-in region. Moreover, we can observe the “presynchronization,” “synchronization,” and “postsynchronization” of vortex-induced vibration curves of square column VIVPEH under different attack angles.


#### 5.3.1.α=0°

The lock-in region is fromUrsqu=6.3 to the end of Ursqu=6.7, and the maximum amplitude is Ysqumax/Dsqu = 0.05, which appeared at Ursqu=6.5, as shown in Figure 4. In order to see more clearly about the vibration amplitude Ysqu/Dsqu, some displacement time-history curves and FFT analysis results for α=0° are given in Figure 5. It can be seen that when Ursqu is small, the amplitude of time-history curve presents a stable sinusoidal curve, and the maximum amplitude is so small that it can be ignored, which means that there is almost no vortex-induced vibration in the system. With the increases of Ursqu, the system enters “presynchronization” phase and the vibration amplitude Ysqu/Dsqu gradually increases. When the vortex shedding frequency approaches the natural frequency of square column VIVPEH, the oscillation amplitude of the square column VIVPEH increases significantly and the system enters synchronous phase, which means the phenomenon of “lock-in” occurs. Under synchronous phase, the vortex shedding frequency is kept constant and the amplitude of the system will remain at a higher value when the flow velocity Ursqu is increased from Ursqu=6.3 to the end of Ursqu=6.7.Figure 5
Displacement time-history and FFT analysis results forα=0°.
(a)
Nondimensional displacement forUrsqu=6.0
(b)
FFT analysis forUrsqu=6.0
(c)
Nondimensional displacement forUrsqu=6.5
(d)
FFT analysis forUrsqu=6.5


#### 5.3.2.α=15°

It can be seen from Figure2(b) that there is a corner at the upper windward of the square column, which can obviously affect the flow field. As is shown in Figure 5, the maximum amplitude decreases and the lock-in region is narrow (from Ursqu=5.4 to Ursqu=5.7). The maximum amplitude is Ysqumax/Dsqu = 0.148, which appeared at Ursqu=5.7. The displacement time-history curves and FFT analysis results for α=15° are given in Figure 6. Similar with the results of α=0°, the displacement time-history curve can also be divided into the following four stages: “unsynchronization,” “presynchronization,” “synchronization,” and “postsynchronization.”Figure 6
Displacement time-history and FFT analysis results forα=15°.
(a)
Nondimensional displacement forUrsqu=5.307
(b)
FFT analysis forUrsqu=5.307
(c)
Nondimensional displacement forUrsqu=5.716
(d)
FFT analysis forUrsqu=5.716


#### 5.3.3.α=30°

The results in Figure4 show that the maximum amplitude of the square column is obviously increased with a wider range of lock-in region (from Ursqu=4.5 to Ursqu=5.7); the maximum amplitude Ysqumax/Dsqu reaches 0.41. Similarly, some of the displacement time-history curves and FFT analysis results for α=30° are given in Figure 7. The results in Figure 7(a) show that the whole time-history curve appears as a parabolic shape and the growth rate decreased gradually to a stable level at late stage, which indicates the coupling of the vortex shedding frequency and natural frequency; that is, the lock-in phenomenon occurs. It can also be seen that there is a “beat” phenomenon in the amplitude curves, as shown in Figure 7(a). There are three peaks in the frequency spectrum curve, as shown in Figure 7(b), which indicates that three harmonics occur in the system. The explaining of the above phenomena is as follows. In the process of vortex-induced vibration of the square column, the vortex shedding frequency couples with the natural frequency when the flow velocity Ursqu reaches a certain value. The flow field near the square column surface is strongly disturbed because of the corner point of the square column, which results in the superposition of three vibration frequencies. This is because the unbalance of the upper and lower aerodynamic force is caused by the obvious asymmetry under attack angle α=30°. When the flow velocity is about Ursqu=5.49, as shown in Figures 7(c) and 7(d), the amplitude curve appears as a complete sine curve without noise. The maximum amplitude Ysqumax/Dsqu is increased to about 0.41 and the vibration frequency is stabilized at 6.98 Hz, which can be regarded as the best working condition of square column VIVPEH.Figure 7
Displacement time-history and FFT analysis results forα=30°.
(a)
Nondimensional displacement forUrsqu=4.758
(b)
FFT analysis forUrsqu=4.758
(c)
Nondimensional displacement forUrsqu=5.49
(d)
FFT analysis forUrsqu=5.49


#### 5.3.4.α=45°

As is shown in Figure4, the displacement time-history curve can also be divided into the following four stages: “unsynchronization,” “presynchronization,” “synchronization,” and “postsynchronization.” The maximum amplitude is Ysqumax/Dsqu=0.28 and the lock-in region is from Ursqu=4.2 to the end of Ursqu=5.4, which appears earlier than that of other attack angles. The displacement time-history curves and FFT analysis results for α=45° are given in Figure 8. It can be seen that there is no severe aerodynamic disturbance, which shows that the vortex-induced vibration response of symmetric bluff body is relatively stable.Figure 8
Displacement time-history and FFT analysis results forα=45°.
(a)
Nondimensional displacement forUrsqu=3.678
(b)
FFT analysis forUrsqu=3.678
(c)
Nondimensional displacement forUrsqu=4.951
(d)
FFT analysis forUrsqu=4.951


#### 5.3.5.α=60°

The displacement time-history curves and FFT analysis results forα=60° are given in Figure 9, which are similar in shape to the one for α=30° with only slight differences in values. According to the displacement time-history curves, as shown in Figures 7 and 9, it can be seen that the amplitude results are almost exactly the same as the one for α=30° in the lock-in region. The same conclusion can also be obtained from the spectral analysis results; that is, the phenomena of harmonic and noise are similar. This is because the spring force of the system and the gravity of the column provide a balance in the flow field. So it can be considered that the above two cases of α=30° and α=60° are symmetric.Figure 9
Displacement time-history and FFT analysis results forα=60°.
(a)
Nondimensional displacement forUrsqu=4.758
(b)
FFT analysis forUrsqu=4.758
(c)
Nondimensional displacement forUrsqu=5.49
(d)
FFT analysis forUrsqu=5.49


#### 5.3.6.α=75°

The results in Figure10 show the displacement time-history curves and FFT analysis results for α=75°. It is easy to see that both the displacement time-history curves and the spectral analysis results are similar to those of α=15° case, which indicates that the cases of α=15° and α=75° are also symmetric.Figure 10
Displacement time-history and FFT analysis results forα=75°.
(a)
Nondimensional displacement forUrsqu=5.307
(b)
FFT analysis forUrsqu=5.307
(c)
Nondimensional displacement forUrsqu=5.716
(d)
FFT analysis forUrsqu=5.716In order to verify the above conclusions, theStrouhal Numbers for the following four cases, α=15°, α=30°, α=60°, and α=75°, are compared with each other, as shown in Table 3. It can be seen that theStrouhal Numberresults for the case of α=15° (α=30°) are almost equal to those for the case of α=75° (α=60°), which indicates that the above analysis is correct.Table 3
Comparison ofStrouhal Number with different attack angle.
U (m/s)
Strouhal Number
15°
75°
30°
60°
0.05049
0.1682
0.1682
0.1832
0.1831
0.0561
0.1714
0.1714
0.1873
0.1872
0.06278
0.1743
0.1742
0.1896
0.1896
0.07293
0.1746
0.1746
0.1922
0.1922To show the difference more clearly of the lock-in region of square column VIVPEH under different attack angles, we choose the dimensionless frequencyfsqu/fnsqu to indicate the frequency characteristics of square column VIVPEH, where fsqu is the vortex shedding frequency of the system, which can be obtained by Fast Fourier Transform (FFT) of the displacement time-history curve, shown in Figures 5–10; fnsqu is the natural frequency of the square column. The results in Figure 11 show the different lock-in regions of square column VIVPEH under different attack angles. For the case of α=0°, the range of values for the lock-in region is 6.5 to 7.0 (6.5≤Ursqu≤7.0); the corresponding bandwidth value is 0.5. For the case of α=15° and α=75°, the range of values for the lock-in region is 5.4 to 5.7 (5.4≤Ursqu≤5.7); the corresponding bandwidth value is just 0.3. For the case of α=30° and α=60°, the range of values for the lock-in region is 4.6 to 5.5 (4.6≤Ursqu≤5.5), the corresponding bandwidth value is 0.9, which is significantly larger than that of the above two cases. When the attack angle is α=45°, the range of values for the lock-in region is 4.2 to 5.4 (4.2≤Ursqu≤5.4); the corresponding bandwidth value increases to 1.2. In addition, the three different branch types of square column VIVPEH can be observed in Figure 11, such as “presynchronization,” “synchronization,” and “postsynchronization.”Figure 11
Nondimensional frequency of square column VIVPEH under different attack angles.Based on the above analysis results, we believe that the calculation parameters and response parameters for the case ofα=15° and α=30° are about the same as those for the case of α=75° and α=60°.


### 5.4. Phase Angle Analysis of Square Column VIVPEH under Different Attack Angles

The analysis results of the initial phase angle and the synchronous phase angle of square column VIVPEH under different attack angles are given in Figures12–15. It can be seen that the vibration is stable in the initial stage with a single amplitude and frequency, which is completely determined by the vortex shedding frequency. Therefore, the displacement time-history curve is almost synchronous with the lift coefficient curve without phase delay. This is because the amplitude is small, so that there is almost no stagnation when the vibration amplitude reaches its maximum or minimum value. Different from the initial vibration state, the vortex shedding frequency is locked again in the synchronous state, resulting in a multiple relationship between the vibration frequency and the natural frequency. That means, in synchronization region, there exists no phase difference between the vortex shedding frequency and vibration frequency. Accordingly, the displacement time-history curve is completely synchronized with the lift coefficient curve with some phase delay.Figure 12
Results of phase angle in initial and synchronization state forα=0°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure 13
Results of phase angle in initial and synchronization state forα=15°/α=75°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure 14
Results of phase angle in initial and synchronization state forα=30°/α=60°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure 15
Results of phase angle in initial and synchronization state forα=45°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure16 shows that the vortex-induced vibration of square column VIVPEH will stagnate for some time at the wave crest or the wave trough, because of the buffering effect of spring. As is shown in Figure 16(a), the vortices of V1 and V2 move forward to a distance of S1 when the wave crests of the two adjacent steps appear. At the moment, V1 obviously becomes thinner and longer, while the vibration of the square column is still in the wave crest. Similarly, the vortices of V3 and V4 move forward to a distance of S2 when the wave troughs of the two adjacent steps appear. It should be pointed out that, the vibration of the square column is always in the wave crest or the wave trough, which leads to the generation of the phase angle.Figure 16
Diagram of phase angle between displacement and lift force development history.
(a)
Y
m
a
x
+ (wave crest)
(b)
Y
m
a
x
- (wave trough)


### 5.5. Analysis of Near Wake Vortex Shedding of Square Column VIVPEH under Different Attack Angles

The shapes of wake vortices in synchronized state of square column VIVPEH under different attack angles are given in Figures17–20, in which T is the vibration period with subscript representing the value of an attack angle. The direction of the negative vorticity region is counterclockwise, which is expressed in blue; while the direction of the positive one is clockwise, which is expressed in red. It can be seen that, for the case of α=0°, the flow field around the square column is stable, the vibration amplitude is small, and the wake vortex structure presents a regular 2S shape. In other words, a positive and negative vortex pair sheds in a cycle. For α=15°/α=75°, the position of near wake vortices shedding of square column VIVPEH began moving forward to the near column wall with increasing vibration amplitudes. Correspondingly, the wake vortex shedding array gradually changes its shape from 2S shape to a circle. As for α=30°/α=60°, it can be observed that the wake vortex shedding mode changes obviously with an increasing width of the vortex array, due to the increase of vibration amplitude and the shape of the vortex pair is changed from a flat shape to a regular elliptical shape. When the attack angle is α=45°, the vibration amplitude decreases, which results in the wake vortex shedding mode changing back to the stable 2S mode.Figure 17
Vortex structures in synchronization region forα=0°.Figure 18
Vortex structures in synchronization region forα=15°/α=75°.Figure 19
Vortex structures in synchronization region forα=30°/α=60°.Figure 20
Vortex structures in synchronization region forα=45°.


### 5.6. Voltage Output and Power Output of Square Column VIVPEH under Different Attack Angles

Due to the influence of different attack angles, the vortex separation point of square column VIVPEH is different from that of cylinder, which results in the following analysis being more complicated. Considering the symmetric nature of the caseα=15°/α=75° and the case α=30°/α=60°, we take only the case of α=0°, α=15°, α=30°, and α=45° in the following analysis of voltage output. When the resistance load is R = 1 × 106 Ω, the output voltage of the system can be calculated by (14).The purpose of this section is to investigate the maximum value of the output voltage and the lock-in region to select the optimal attack angle of square column VIVPEH. The results of the maximum output voltage of square column VIVPEH under different attack angles and the effective working area of synchronization are given in Figure21. It can be seen that the maximum output voltage of square column VIVPEH appears at the case of α=45° and the corresponding value is 6.732 V, while the minimum value appears at the case of α=15°/α=75°, which indicates that these two cases are not suitable for energy harvesting. Similarly, the output voltage of the case of α=0° is slightly higher than that of the case of α=15°/α=75°; however, its value is still small. When the attack angle increases from α=15° to α=45°, the output voltage of square column VIVPEH system increases to its maximum value. The results of the working region of synchronization in Figure 21 show that the total bandwidth of the above six cases is (4.2–5.7, 5.3–5.7), and there are about 0.6 bandwidth of the region of nonsynchronization. The maximum bandwidth is (4.2–5.4) under attack angle α=45°, which is higher than that of all other attack angle cases. Accordingly, the output power of the system can be calculated by (15).Figure 21
Comparison of maximum value of voltage output and synchronization effective area between VIVPEH with different shapes.The results of the maximum output power of square column VIVPEH under different attack angles are given in Figure22. It can also be seen that the maximum value of the output power appears at the case of α=45° and the corresponding value is 4.5 × 10−5 W. Therefore, square column VIVPEH under attack angle α=45° is an ideal PEH, because of its higher voltage output value and larger working bandwidth.Figure 22
Comparison of maximum value of power output between VIVPEH with different shapes.


## 5.1. Analysis Parameters and Cases

To carry out the numerical calculation of vortex-induced vibration, this paper provides the parameters of the vibration energy harvesting system, as shown in Table1, in which Dsqu is the side length of the square column, DNor stands for the dimensionless characteristic length of the square column, and the relationship between DNor and Dsqu is given in(16)DNor=Dsqusin⁡α+cos⁡α.Based on the parameters of Table 1 and formula (16), we can calculate the numerical values of DNor as 0.0016 under 0° attack angle, 0.00196 under 15° attack angle, 0.00219 under 30° attack angle, 0.00226 under 45° attack angle, 0.00219 under 60° attack angle, and 0.00196 under 75° attack angle.Table 1
Calculation parameters for square column VIVPEH system.
Symbols and units
Physical meanings
Numerical values
M [kg]
Square column mass
0.2979
K [N/m]
Elastic coefficient of the system
579~584
C [N·s/m]
Damping of the system
0.0325~0.45
D
s
q
u [m]
Square column side length
0.0016
α [degree]
Angle between the square column and the inflow velocity
0~75
ζ
Damping ratio of the system
0.00121
f
n [Hz]
Natural frequency of the system
7.012~7.132 (water)
C
p [nF]
Capacitance
120
θ
Electromechanical coupling coefficient
1.55 × 10−3
μ [Pa·s]
Dynamic viscosity
0.0011379
ν [m2/s]
Kinematic viscosity
1.139 × 10−6 (water)
ρ [kg/m3]
Density
999.1026 (water)The computational grid and boundary conditions of square column vortex-induced vibration under different attack angles are given in Figure2, where the computational domain is 20 × 20D, and the entire domain includes five boundaries: velocity inlet, velocity outlet, top, bottom, and a column wall. The inlet velocity is considered as uniform and constant velocity. For outlet boundary, a zero gradient condition is specified for velocity. The top and bottom condition are defined as a wall boundary. In present numerical study, a moving wall boundary condition is applied for the square column when the column is in VIV. The two-dimensional, structured grids were generated with the help of “Gambit” software. The grid domain size is 20 × 20D. The square column was set in the center in the domain to ensure that the results of the numerical model are accurate. The conditions at the outlet are close to the assumed conditions. The computational domain in the vicinity of each cylinder is a 3 × 3D square where the grid density for the near-wall region is enhanced to solve for high resolution in flow properties.Figure 2
Computational grid and boundary conditions.
(a)
α
=
0°
(b)
α
=
15°
(c)
α
=
30°
(d)
α
=
45°
(e)
α
=
60°
(f)
α
=
75°
(g)When the external resistance value isR = 1 × 106 Ω, the spring stiffness value of the square column VIVPEH system is Ksqu=580 N/m, and the corresponding damping value is Csqu=0.2 Ns/m. Taking the range values of the flow velocity as 0.03927 m/s to 0.08975 m/s, we can calculate the numerical values of Ursqu under different attack angles, as shown in Table 2, in which Ursqu is the reduction velocity, defined as follows:(17)Ursqu=UfnDsqu,where fn=ωn/2π.Table 2
Run case of square column in OpenFOAM.
U (m/s)
Ursqu
α
=
0°
α
=
15
° (α=75°)
α
=
30
° (α=60°)
α=45°
0.03927
3.5
2.85784
2.56223
2.47525
0.04488
4
3.26611
2.92826
2.82885
0.05049
4.5
3.67437
3.29429
3.18246
0.0561
5
4.08263
3.66032
3.53607
0.05834
5.2
4.24594
3.80673
3.67751
0.06059
5.4
4.40924
3.95315
3.81895
0.06283
5.6
4.57255
4.09956
3.9604
0.06732
6
4.89916
4.39239
4.24328
0.07293
6.5
5.30742
4.75842
4.59689
0.07854
7
5.71569
5.12445
4.9505
0.08415
7.5
6.12395
5.49048
5.3041
0.08976
8
6.53221
5.85652
5.65771
0.09537
8.5
6.94048
6.22255
6.01132In this paper, the matrix method is used to evaluate the influence of the external resistance load on the damping and natural frequency of the electromechanical coupling system by MATLAB. Then, the damping and natural frequency values can be used for initial conditions in OpenFOAM to compute the vibration amplitude of the vortex-induced vibration system under different Reynolds numbers(94<Re<115). The main focus of this paper is to investigate the influence of external resistance load (R = 1 × 103 Ω, 1 × 104 Ω, 1 × 105 Ω, 1 × 106 Ω, 1 × 107 Ω) on the energy conversion of the square column VIVPEH system, including the effect on vibration amplitude, output voltage, and power.


## 5.2. System Damping and Natural Frequency Characteristics

The damping and natural frequency of the electromechanical coupling system can be obtained by MATLAB, and the real and imaginary parts of the circuit conjugate solutions are shown in Figure3.Figure 3
Real and imaginary parts of the circuit conjugate solutions.According to Figure3, we can see that the total damping of the system is small when the resistance load is small. In the case of R<1×105 Ω, the system total damping is increased with the resistance load increasing. Once the resistance load R reaches 1 × 105 Ω, the system total damping reaches the maximum value. It is noteworthy that the total damping of the system decreases instead of increasing when the resistance load keeps increasing; that is, R>1×105 Ω. For the natural circular frequency, when R<3×104 Ω, the frequency value remains at 44 rad/s; when R>2×106 Ω, the frequency value approximately remains at 50 rad/s. Generally, the natural circular frequency value of the system is relatively stable, which is kept in the range of 44–50 rad/s.


## 5.3. Vibration Characteristics of Square Column VIVPEH under Different Attack Angles

In the following numerical simulations, we take the resistance loadR = 1 × 106 Ω and compute the amplitude response of the square column VIVPEH under different attack angles and different flow velocities, as shown in Figure 4. Note here that the natural frequency of square column is fn=6.98 Hz, when the resistance load is R = 1 × 106 Ω. The dimensionless vibration amplitude Ysqu/Dsqu is adopted to indicate the vibration response of the square column. The maximum value of Ysqu can be obtained by means of averaging the peak value of at least 60 displacement time-history response curves. In this paper, the vibration amplitude of the smooth circular column is provided to compare with that of the square column.Figure 4
Dimensionless vibration amplitude of square column VIVPEH under different attack angles and different flow velocities.As we can see in Figure4, different attack angles obviously have an effect on the peak vibration amplitude and lock-in region. Moreover, we can observe the “presynchronization,” “synchronization,” and “postsynchronization” of vortex-induced vibration curves of square column VIVPEH under different attack angles.


### 5.3.1.α=0°

The lock-in region is fromUrsqu=6.3 to the end of Ursqu=6.7, and the maximum amplitude is Ysqumax/Dsqu = 0.05, which appeared at Ursqu=6.5, as shown in Figure 4. In order to see more clearly about the vibration amplitude Ysqu/Dsqu, some displacement time-history curves and FFT analysis results for α=0° are given in Figure 5. It can be seen that when Ursqu is small, the amplitude of time-history curve presents a stable sinusoidal curve, and the maximum amplitude is so small that it can be ignored, which means that there is almost no vortex-induced vibration in the system. With the increases of Ursqu, the system enters “presynchronization” phase and the vibration amplitude Ysqu/Dsqu gradually increases. When the vortex shedding frequency approaches the natural frequency of square column VIVPEH, the oscillation amplitude of the square column VIVPEH increases significantly and the system enters synchronous phase, which means the phenomenon of “lock-in” occurs. Under synchronous phase, the vortex shedding frequency is kept constant and the amplitude of the system will remain at a higher value when the flow velocity Ursqu is increased from Ursqu=6.3 to the end of Ursqu=6.7.Figure 5
Displacement time-history and FFT analysis results forα=0°.
(a)
Nondimensional displacement forUrsqu=6.0
(b)
FFT analysis forUrsqu=6.0
(c)
Nondimensional displacement forUrsqu=6.5
(d)
FFT analysis forUrsqu=6.5


### 5.3.2.α=15°

It can be seen from Figure2(b) that there is a corner at the upper windward of the square column, which can obviously affect the flow field. As is shown in Figure 5, the maximum amplitude decreases and the lock-in region is narrow (from Ursqu=5.4 to Ursqu=5.7). The maximum amplitude is Ysqumax/Dsqu = 0.148, which appeared at Ursqu=5.7. The displacement time-history curves and FFT analysis results for α=15° are given in Figure 6. Similar with the results of α=0°, the displacement time-history curve can also be divided into the following four stages: “unsynchronization,” “presynchronization,” “synchronization,” and “postsynchronization.”Figure 6
Displacement time-history and FFT analysis results forα=15°.
(a)
Nondimensional displacement forUrsqu=5.307
(b)
FFT analysis forUrsqu=5.307
(c)
Nondimensional displacement forUrsqu=5.716
(d)
FFT analysis forUrsqu=5.716


### 5.3.3.α=30°

The results in Figure4 show that the maximum amplitude of the square column is obviously increased with a wider range of lock-in region (from Ursqu=4.5 to Ursqu=5.7); the maximum amplitude Ysqumax/Dsqu reaches 0.41. Similarly, some of the displacement time-history curves and FFT analysis results for α=30° are given in Figure 7. The results in Figure 7(a) show that the whole time-history curve appears as a parabolic shape and the growth rate decreased gradually to a stable level at late stage, which indicates the coupling of the vortex shedding frequency and natural frequency; that is, the lock-in phenomenon occurs. It can also be seen that there is a “beat” phenomenon in the amplitude curves, as shown in Figure 7(a). There are three peaks in the frequency spectrum curve, as shown in Figure 7(b), which indicates that three harmonics occur in the system. The explaining of the above phenomena is as follows. In the process of vortex-induced vibration of the square column, the vortex shedding frequency couples with the natural frequency when the flow velocity Ursqu reaches a certain value. The flow field near the square column surface is strongly disturbed because of the corner point of the square column, which results in the superposition of three vibration frequencies. This is because the unbalance of the upper and lower aerodynamic force is caused by the obvious asymmetry under attack angle α=30°. When the flow velocity is about Ursqu=5.49, as shown in Figures 7(c) and 7(d), the amplitude curve appears as a complete sine curve without noise. The maximum amplitude Ysqumax/Dsqu is increased to about 0.41 and the vibration frequency is stabilized at 6.98 Hz, which can be regarded as the best working condition of square column VIVPEH.Figure 7
Displacement time-history and FFT analysis results forα=30°.
(a)
Nondimensional displacement forUrsqu=4.758
(b)
FFT analysis forUrsqu=4.758
(c)
Nondimensional displacement forUrsqu=5.49
(d)
FFT analysis forUrsqu=5.49


### 5.3.4.α=45°

As is shown in Figure4, the displacement time-history curve can also be divided into the following four stages: “unsynchronization,” “presynchronization,” “synchronization,” and “postsynchronization.” The maximum amplitude is Ysqumax/Dsqu=0.28 and the lock-in region is from Ursqu=4.2 to the end of Ursqu=5.4, which appears earlier than that of other attack angles. The displacement time-history curves and FFT analysis results for α=45° are given in Figure 8. It can be seen that there is no severe aerodynamic disturbance, which shows that the vortex-induced vibration response of symmetric bluff body is relatively stable.Figure 8
Displacement time-history and FFT analysis results forα=45°.
(a)
Nondimensional displacement forUrsqu=3.678
(b)
FFT analysis forUrsqu=3.678
(c)
Nondimensional displacement forUrsqu=4.951
(d)
FFT analysis forUrsqu=4.951


### 5.3.5.α=60°

The displacement time-history curves and FFT analysis results forα=60° are given in Figure 9, which are similar in shape to the one for α=30° with only slight differences in values. According to the displacement time-history curves, as shown in Figures 7 and 9, it can be seen that the amplitude results are almost exactly the same as the one for α=30° in the lock-in region. The same conclusion can also be obtained from the spectral analysis results; that is, the phenomena of harmonic and noise are similar. This is because the spring force of the system and the gravity of the column provide a balance in the flow field. So it can be considered that the above two cases of α=30° and α=60° are symmetric.Figure 9
Displacement time-history and FFT analysis results forα=60°.
(a)
Nondimensional displacement forUrsqu=4.758
(b)
FFT analysis forUrsqu=4.758
(c)
Nondimensional displacement forUrsqu=5.49
(d)
FFT analysis forUrsqu=5.49


### 5.3.6.α=75°

The results in Figure10 show the displacement time-history curves and FFT analysis results for α=75°. It is easy to see that both the displacement time-history curves and the spectral analysis results are similar to those of α=15° case, which indicates that the cases of α=15° and α=75° are also symmetric.Figure 10
Displacement time-history and FFT analysis results forα=75°.
(a)
Nondimensional displacement forUrsqu=5.307
(b)
FFT analysis forUrsqu=5.307
(c)
Nondimensional displacement forUrsqu=5.716
(d)
FFT analysis forUrsqu=5.716In order to verify the above conclusions, theStrouhal Numbers for the following four cases, α=15°, α=30°, α=60°, and α=75°, are compared with each other, as shown in Table 3. It can be seen that theStrouhal Numberresults for the case of α=15° (α=30°) are almost equal to those for the case of α=75° (α=60°), which indicates that the above analysis is correct.Table 3
Comparison ofStrouhal Number with different attack angle.
U (m/s)
Strouhal Number
15°
75°
30°
60°
0.05049
0.1682
0.1682
0.1832
0.1831
0.0561
0.1714
0.1714
0.1873
0.1872
0.06278
0.1743
0.1742
0.1896
0.1896
0.07293
0.1746
0.1746
0.1922
0.1922To show the difference more clearly of the lock-in region of square column VIVPEH under different attack angles, we choose the dimensionless frequencyfsqu/fnsqu to indicate the frequency characteristics of square column VIVPEH, where fsqu is the vortex shedding frequency of the system, which can be obtained by Fast Fourier Transform (FFT) of the displacement time-history curve, shown in Figures 5–10; fnsqu is the natural frequency of the square column. The results in Figure 11 show the different lock-in regions of square column VIVPEH under different attack angles. For the case of α=0°, the range of values for the lock-in region is 6.5 to 7.0 (6.5≤Ursqu≤7.0); the corresponding bandwidth value is 0.5. For the case of α=15° and α=75°, the range of values for the lock-in region is 5.4 to 5.7 (5.4≤Ursqu≤5.7); the corresponding bandwidth value is just 0.3. For the case of α=30° and α=60°, the range of values for the lock-in region is 4.6 to 5.5 (4.6≤Ursqu≤5.5), the corresponding bandwidth value is 0.9, which is significantly larger than that of the above two cases. When the attack angle is α=45°, the range of values for the lock-in region is 4.2 to 5.4 (4.2≤Ursqu≤5.4); the corresponding bandwidth value increases to 1.2. In addition, the three different branch types of square column VIVPEH can be observed in Figure 11, such as “presynchronization,” “synchronization,” and “postsynchronization.”Figure 11
Nondimensional frequency of square column VIVPEH under different attack angles.Based on the above analysis results, we believe that the calculation parameters and response parameters for the case ofα=15° and α=30° are about the same as those for the case of α=75° and α=60°.


## 5.3.1.α=0°

The lock-in region is fromUrsqu=6.3 to the end of Ursqu=6.7, and the maximum amplitude is Ysqumax/Dsqu = 0.05, which appeared at Ursqu=6.5, as shown in Figure 4. In order to see more clearly about the vibration amplitude Ysqu/Dsqu, some displacement time-history curves and FFT analysis results for α=0° are given in Figure 5. It can be seen that when Ursqu is small, the amplitude of time-history curve presents a stable sinusoidal curve, and the maximum amplitude is so small that it can be ignored, which means that there is almost no vortex-induced vibration in the system. With the increases of Ursqu, the system enters “presynchronization” phase and the vibration amplitude Ysqu/Dsqu gradually increases. When the vortex shedding frequency approaches the natural frequency of square column VIVPEH, the oscillation amplitude of the square column VIVPEH increases significantly and the system enters synchronous phase, which means the phenomenon of “lock-in” occurs. Under synchronous phase, the vortex shedding frequency is kept constant and the amplitude of the system will remain at a higher value when the flow velocity Ursqu is increased from Ursqu=6.3 to the end of Ursqu=6.7.Figure 5
Displacement time-history and FFT analysis results forα=0°.
(a)
Nondimensional displacement forUrsqu=6.0
(b)
FFT analysis forUrsqu=6.0
(c)
Nondimensional displacement forUrsqu=6.5
(d)
FFT analysis forUrsqu=6.5


## 5.3.2.α=15°

It can be seen from Figure2(b) that there is a corner at the upper windward of the square column, which can obviously affect the flow field. As is shown in Figure 5, the maximum amplitude decreases and the lock-in region is narrow (from Ursqu=5.4 to Ursqu=5.7). The maximum amplitude is Ysqumax/Dsqu = 0.148, which appeared at Ursqu=5.7. The displacement time-history curves and FFT analysis results for α=15° are given in Figure 6. Similar with the results of α=0°, the displacement time-history curve can also be divided into the following four stages: “unsynchronization,” “presynchronization,” “synchronization,” and “postsynchronization.”Figure 6
Displacement time-history and FFT analysis results forα=15°.
(a)
Nondimensional displacement forUrsqu=5.307
(b)
FFT analysis forUrsqu=5.307
(c)
Nondimensional displacement forUrsqu=5.716
(d)
FFT analysis forUrsqu=5.716


## 5.3.3.α=30°

The results in Figure4 show that the maximum amplitude of the square column is obviously increased with a wider range of lock-in region (from Ursqu=4.5 to Ursqu=5.7); the maximum amplitude Ysqumax/Dsqu reaches 0.41. Similarly, some of the displacement time-history curves and FFT analysis results for α=30° are given in Figure 7. The results in Figure 7(a) show that the whole time-history curve appears as a parabolic shape and the growth rate decreased gradually to a stable level at late stage, which indicates the coupling of the vortex shedding frequency and natural frequency; that is, the lock-in phenomenon occurs. It can also be seen that there is a “beat” phenomenon in the amplitude curves, as shown in Figure 7(a). There are three peaks in the frequency spectrum curve, as shown in Figure 7(b), which indicates that three harmonics occur in the system. The explaining of the above phenomena is as follows. In the process of vortex-induced vibration of the square column, the vortex shedding frequency couples with the natural frequency when the flow velocity Ursqu reaches a certain value. The flow field near the square column surface is strongly disturbed because of the corner point of the square column, which results in the superposition of three vibration frequencies. This is because the unbalance of the upper and lower aerodynamic force is caused by the obvious asymmetry under attack angle α=30°. When the flow velocity is about Ursqu=5.49, as shown in Figures 7(c) and 7(d), the amplitude curve appears as a complete sine curve without noise. The maximum amplitude Ysqumax/Dsqu is increased to about 0.41 and the vibration frequency is stabilized at 6.98 Hz, which can be regarded as the best working condition of square column VIVPEH.Figure 7
Displacement time-history and FFT analysis results forα=30°.
(a)
Nondimensional displacement forUrsqu=4.758
(b)
FFT analysis forUrsqu=4.758
(c)
Nondimensional displacement forUrsqu=5.49
(d)
FFT analysis forUrsqu=5.49


## 5.3.4.α=45°

As is shown in Figure4, the displacement time-history curve can also be divided into the following four stages: “unsynchronization,” “presynchronization,” “synchronization,” and “postsynchronization.” The maximum amplitude is Ysqumax/Dsqu=0.28 and the lock-in region is from Ursqu=4.2 to the end of Ursqu=5.4, which appears earlier than that of other attack angles. The displacement time-history curves and FFT analysis results for α=45° are given in Figure 8. It can be seen that there is no severe aerodynamic disturbance, which shows that the vortex-induced vibration response of symmetric bluff body is relatively stable.Figure 8
Displacement time-history and FFT analysis results forα=45°.
(a)
Nondimensional displacement forUrsqu=3.678
(b)
FFT analysis forUrsqu=3.678
(c)
Nondimensional displacement forUrsqu=4.951
(d)
FFT analysis forUrsqu=4.951


## 5.3.5.α=60°

The displacement time-history curves and FFT analysis results forα=60° are given in Figure 9, which are similar in shape to the one for α=30° with only slight differences in values. According to the displacement time-history curves, as shown in Figures 7 and 9, it can be seen that the amplitude results are almost exactly the same as the one for α=30° in the lock-in region. The same conclusion can also be obtained from the spectral analysis results; that is, the phenomena of harmonic and noise are similar. This is because the spring force of the system and the gravity of the column provide a balance in the flow field. So it can be considered that the above two cases of α=30° and α=60° are symmetric.Figure 9
Displacement time-history and FFT analysis results forα=60°.
(a)
Nondimensional displacement forUrsqu=4.758
(b)
FFT analysis forUrsqu=4.758
(c)
Nondimensional displacement forUrsqu=5.49
(d)
FFT analysis forUrsqu=5.49


## 5.3.6.α=75°

The results in Figure10 show the displacement time-history curves and FFT analysis results for α=75°. It is easy to see that both the displacement time-history curves and the spectral analysis results are similar to those of α=15° case, which indicates that the cases of α=15° and α=75° are also symmetric.Figure 10
Displacement time-history and FFT analysis results forα=75°.
(a)
Nondimensional displacement forUrsqu=5.307
(b)
FFT analysis forUrsqu=5.307
(c)
Nondimensional displacement forUrsqu=5.716
(d)
FFT analysis forUrsqu=5.716In order to verify the above conclusions, theStrouhal Numbers for the following four cases, α=15°, α=30°, α=60°, and α=75°, are compared with each other, as shown in Table 3. It can be seen that theStrouhal Numberresults for the case of α=15° (α=30°) are almost equal to those for the case of α=75° (α=60°), which indicates that the above analysis is correct.Table 3
Comparison ofStrouhal Number with different attack angle.
U (m/s)
Strouhal Number
15°
75°
30°
60°
0.05049
0.1682
0.1682
0.1832
0.1831
0.0561
0.1714
0.1714
0.1873
0.1872
0.06278
0.1743
0.1742
0.1896
0.1896
0.07293
0.1746
0.1746
0.1922
0.1922To show the difference more clearly of the lock-in region of square column VIVPEH under different attack angles, we choose the dimensionless frequencyfsqu/fnsqu to indicate the frequency characteristics of square column VIVPEH, where fsqu is the vortex shedding frequency of the system, which can be obtained by Fast Fourier Transform (FFT) of the displacement time-history curve, shown in Figures 5–10; fnsqu is the natural frequency of the square column. The results in Figure 11 show the different lock-in regions of square column VIVPEH under different attack angles. For the case of α=0°, the range of values for the lock-in region is 6.5 to 7.0 (6.5≤Ursqu≤7.0); the corresponding bandwidth value is 0.5. For the case of α=15° and α=75°, the range of values for the lock-in region is 5.4 to 5.7 (5.4≤Ursqu≤5.7); the corresponding bandwidth value is just 0.3. For the case of α=30° and α=60°, the range of values for the lock-in region is 4.6 to 5.5 (4.6≤Ursqu≤5.5), the corresponding bandwidth value is 0.9, which is significantly larger than that of the above two cases. When the attack angle is α=45°, the range of values for the lock-in region is 4.2 to 5.4 (4.2≤Ursqu≤5.4); the corresponding bandwidth value increases to 1.2. In addition, the three different branch types of square column VIVPEH can be observed in Figure 11, such as “presynchronization,” “synchronization,” and “postsynchronization.”Figure 11
Nondimensional frequency of square column VIVPEH under different attack angles.Based on the above analysis results, we believe that the calculation parameters and response parameters for the case ofα=15° and α=30° are about the same as those for the case of α=75° and α=60°.


## 5.4. Phase Angle Analysis of Square Column VIVPEH under Different Attack Angles

The analysis results of the initial phase angle and the synchronous phase angle of square column VIVPEH under different attack angles are given in Figures12–15. It can be seen that the vibration is stable in the initial stage with a single amplitude and frequency, which is completely determined by the vortex shedding frequency. Therefore, the displacement time-history curve is almost synchronous with the lift coefficient curve without phase delay. This is because the amplitude is small, so that there is almost no stagnation when the vibration amplitude reaches its maximum or minimum value. Different from the initial vibration state, the vortex shedding frequency is locked again in the synchronous state, resulting in a multiple relationship between the vibration frequency and the natural frequency. That means, in synchronization region, there exists no phase difference between the vortex shedding frequency and vibration frequency. Accordingly, the displacement time-history curve is completely synchronized with the lift coefficient curve with some phase delay.Figure 12
Results of phase angle in initial and synchronization state forα=0°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure 13
Results of phase angle in initial and synchronization state forα=15°/α=75°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure 14
Results of phase angle in initial and synchronization state forα=30°/α=60°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure 15
Results of phase angle in initial and synchronization state forα=45°.
(a)
Initial phase angle
(b)
Synchronous phase angleFigure16 shows that the vortex-induced vibration of square column VIVPEH will stagnate for some time at the wave crest or the wave trough, because of the buffering effect of spring. As is shown in Figure 16(a), the vortices of V1 and V2 move forward to a distance of S1 when the wave crests of the two adjacent steps appear. At the moment, V1 obviously becomes thinner and longer, while the vibration of the square column is still in the wave crest. Similarly, the vortices of V3 and V4 move forward to a distance of S2 when the wave troughs of the two adjacent steps appear. It should be pointed out that, the vibration of the square column is always in the wave crest or the wave trough, which leads to the generation of the phase angle.Figure 16
Diagram of phase angle between displacement and lift force development history.
(a)
Y
m
a
x
+ (wave crest)
(b)
Y
m
a
x
- (wave trough)


## 5.5. Analysis of Near Wake Vortex Shedding of Square Column VIVPEH under Different Attack Angles

The shapes of wake vortices in synchronized state of square column VIVPEH under different attack angles are given in Figures17–20, in which T is the vibration period with subscript representing the value of an attack angle. The direction of the negative vorticity region is counterclockwise, which is expressed in blue; while the direction of the positive one is clockwise, which is expressed in red. It can be seen that, for the case of α=0°, the flow field around the square column is stable, the vibration amplitude is small, and the wake vortex structure presents a regular 2S shape. In other words, a positive and negative vortex pair sheds in a cycle. For α=15°/α=75°, the position of near wake vortices shedding of square column VIVPEH began moving forward to the near column wall with increasing vibration amplitudes. Correspondingly, the wake vortex shedding array gradually changes its shape from 2S shape to a circle. As for α=30°/α=60°, it can be observed that the wake vortex shedding mode changes obviously with an increasing width of the vortex array, due to the increase of vibration amplitude and the shape of the vortex pair is changed from a flat shape to a regular elliptical shape. When the attack angle is α=45°, the vibration amplitude decreases, which results in the wake vortex shedding mode changing back to the stable 2S mode.Figure 17
Vortex structures in synchronization region forα=0°.Figure 18
Vortex structures in synchronization region forα=15°/α=75°.Figure 19
Vortex structures in synchronization region forα=30°/α=60°.Figure 20
Vortex structures in synchronization region forα=45°.


## 5.6. Voltage Output and Power Output of Square Column VIVPEH under Different Attack Angles

Due to the influence of different attack angles, the vortex separation point of square column VIVPEH is different from that of cylinder, which results in the following analysis being more complicated. Considering the symmetric nature of the caseα=15°/α=75° and the case α=30°/α=60°, we take only the case of α=0°, α=15°, α=30°, and α=45° in the following analysis of voltage output. When the resistance load is R = 1 × 106 Ω, the output voltage of the system can be calculated by (14).The purpose of this section is to investigate the maximum value of the output voltage and the lock-in region to select the optimal attack angle of square column VIVPEH. The results of the maximum output voltage of square column VIVPEH under different attack angles and the effective working area of synchronization are given in Figure21. It can be seen that the maximum output voltage of square column VIVPEH appears at the case of α=45° and the corresponding value is 6.732 V, while the minimum value appears at the case of α=15°/α=75°, which indicates that these two cases are not suitable for energy harvesting. Similarly, the output voltage of the case of α=0° is slightly higher than that of the case of α=15°/α=75°; however, its value is still small. When the attack angle increases from α=15° to α=45°, the output voltage of square column VIVPEH system increases to its maximum value. The results of the working region of synchronization in Figure 21 show that the total bandwidth of the above six cases is (4.2–5.7, 5.3–5.7), and there are about 0.6 bandwidth of the region of nonsynchronization. The maximum bandwidth is (4.2–5.4) under attack angle α=45°, which is higher than that of all other attack angle cases. Accordingly, the output power of the system can be calculated by (15).Figure 21
Comparison of maximum value of voltage output and synchronization effective area between VIVPEH with different shapes.The results of the maximum output power of square column VIVPEH under different attack angles are given in Figure22. It can also be seen that the maximum value of the output power appears at the case of α=45° and the corresponding value is 4.5 × 10−5 W. Therefore, square column VIVPEH under attack angle α=45° is an ideal PEH, because of its higher voltage output value and larger working bandwidth.Figure 22
Comparison of maximum value of power output between VIVPEH with different shapes.


## 6. Conclusions

The energy harvesting features of square column VIVPEH under different attack angles are investigated in this paper with considering the vibration characteristics, phase characteristics, the near wake vortex shedding mode, and the output voltage and power of the system. The main conclusions are as follows:(1)
Within the range of reduced velocity studied in this paper, the vortex-induced vibration curves of square column VIVPEH under different attack angles can be obtained, which contains the phenomena of “presynchronization,” “synchronization,” and “postsynchronization.” The vortex shedding shape of square column VIVPEH is dominated by 2S mode.(2)
The attack angle has significant effect on the maximum value of vibration amplitude of square column and the lock-in region. Due to the influence of different attack angles, the boundary layer separation point does not move backwards like the cylinder. The maximum vibration amplitude and the lock-in vibration region of square column will also fluctuate. When the attack angle is equal to 45 degrees, the synchronization region can reach a value of 1.2 times of reduction velocity.(3)
The numerical analysis shows that the cases ofα=15°/α=30° and α=75°/α=60° for square column VIVPEH are symmetric, which indicates that the calculated results of the case of α=15°/α=30° are equal to the results of case of α=75°/α=60°, respectively.(4)
The maximum value of the output power is similar to that of the output voltage, which appears at the case ofα=45° and the corresponding value is 4.5 × 10−5 W and 6.732 V, respectively. Therefore, α=45° is a relatively ideal attack angle of square column VIVPEH.



---

*Source: 6439401-2017-08-10.xml*