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

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.


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][2][3][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][6][7][8][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 vortexinduced 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 vortexinduced 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" Geofluids 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 vortexinduced 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 vortexinduced 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.

Physical Model
The physical model of square column VIVPEH is shown in Figure 1, in which  stands for the inflow velocity,  stands for the angle between the square column and the inflow velocity,  is the mass of the column,  is the elastic coefficient of the system,  is the damping of the system, and  is the external resistance load.

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 the Navier-Stokes equation to describe the vortexinduced vibration, uses a linear second-order differential equation to describe the vortex-induced vibration of the single degree of freedom -- (mass spring damping) system, and finally uses coupled Gauss law and vibration equations to describe the electromechanical coupling system.

Fluid Solid Coupling Model.
The external flow field is calculated by the continuity equation and the Navier-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 timedependent viscous flow solutions can be obtained by numerical approximation of the incompressible unsteady Reynolds-Averaged Navier-Stokes (URANS) equations in conjunction with the one-equation Spalart-Allmaras (S-A) turbulence model [33], where a second-order Gauss integration scheme with a linear interpolation is used in the governing equations for the divergence, gradient, and Laplacian terms.For time integration, second-order backwards Euler method is employed.The numerical discretization scheme has secondorder 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 where  is pressure,  is fluid density, ] is dynamic viscosity,   is the mean flow velocity vector, and   is the strain rate tensor, 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       is the Reynolds stress tensor and can be modeled as       = 2    , where   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:

Geofluids
In S-A model, the turbulence eddy viscosity   can be obtained by in which  V1 =  3 /( 3 +  3 ]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  1 and  2 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 freedom -- system can be represented by the following linear second-order differential equation: The relationships between , , and  are as follows: where   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, Ẏ and Ÿ represent the firstor second-order derivative of the column vibration displacement, respectively,   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.

Electromechanical Coupling Model.
In order to describe the relationship between the amplitude and the voltage in the vortex-induced vibration circuit, Gauss law is adopted in this paper.Theoretical derivations are as follows: where  is the electromechanical coupling coefficient,   is the capacitance coefficient, and  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 twoorder nonhomogeneous ordinary differential equation (7).The homogeneous equation of ( 7) is as follows: Let  1 = ,  2 = Ẏ, and  3 = ; substitute (6) into ( 8) and ( 9); we can get The above equations can be expressed in the following matrix form: where The matrix () has three different eigenvalues of   , of which  = 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  1 and  2 , 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  3 is negative constant, we consider only the real part of  1 and  2 , when computing the trivial solution of the matrix ().

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: where  max 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.The corresponding output power can be obtained by the following equation: The above mathematical expressions include the solving process of fluid solid coupling and electromechanical coupling.
First of all, we can get the flow field pressures  and   by solving the URANS equations in OpenFOAM.Secondly, we can obtain the column vibration displacement , the damping , and natural circular frequency   , 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-vibrationelectromechanical coupling system.

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 Table 1, in which  squ is the side length of the square column,  Nor stands for the dimensionless characteristic length of the square column, and the relationship between  Nor and  squ is given in Based on the parameters of  2, 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 twodimensional, 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.
When the external resistance value is  = 1 × 10 6 Ω, the spring stiffness value of the square column VIVPEH system is  squ = 580 N/m, and the corresponding damping value is  squ = 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  squ under different attack angles, as shown in Table 2, in which  squ is the reduction velocity, defined as follows: where   =   /2.In 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 <   < 115).The main focus of this paper is to investigate the influence of external resistance load ( = 1 × 10 3 Ω, 1 × 10 4 Ω, 1 × 10 5 Ω, 1 × 10 6 Ω, 1 × 10 7 Ω) on the energy conversion of the square column VIVPEH system, including the effect on vibration amplitude, output voltage, and power.

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 Figure 3.
According to Figure 3, we can see that the total damping of the system is small when the resistance load is small.In the case of  < 1 × 10 5 Ω, the system total damping is increased with the resistance load increasing.Once the resistance load  reaches 1 × 10 5 Ω, 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,  > 1 × 10 5 Ω.For the natural circular frequency, when  < 3 × 10 4 Ω, the frequency value remains at 44 rad/s; when  > 2 × 10 6 Ω, 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.

Vibration Characteristics of Square Column VIVPEH
under Different Attack Angles.In the following numerical simulations, we take the resistance load  = 1 × 10 6 Ω 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   = 6.98 Hz, when the resistance load is  = 1 × 10 6 Ω.The dimensionless vibration amplitude  squ / squ is adopted to indicate the vibration response of the square column.The maximum value of  squ 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.
As we can see in Figure 4, 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 from  squ = 6.3 to the end of  squ = 6.7, and the maximum amplitude is  squmax /  squ = 0.05, which appeared at  squ = 6.5, as shown in Figure 4.In order to see more clearly about the vibration amplitude  squ / squ , some displacement time-history curves and FFT analysis results for  = 0 ∘ are given in Figure 5.It can be seen that when  squ 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  squ , the system enters "presynchronization" phase and the vibration amplitude  squ / squ 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  squ is increased from  squ = 6.3 to the end of  squ = 6.7.

𝛼 = 15
∘ .It can be seen from Figure 2(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 lockin region is narrow (from  squ = 5.4 to  squ = 5.7).The maximum amplitude is  squmax / squ = 0.148, which appeared at  squ = 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."

𝛼 = 30
∘ .The results in Figure 4 show that the maximum amplitude of the square column is obviously increased with a wider range of lock-in region (from  squ = 4.5 to  squ = 5.7); the maximum amplitude  squmax / squ 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 timehistory 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  squ 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  squ = 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  squmax / squ 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.

𝛼 = 45
∘ .As is shown in Figure 4, the displacement time-history curve can also be divided into the following four stages: "unsynchronization," "presynchronization," "synchronization," and "postsynchronization." The maximum amplitude is  squmax / squ = 0.28 and the lock-in region is from  squ = 4.2 to the end of  squ = 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.

𝛼 = 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 timehistory 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.

𝛼 = 75
∘ .The results in Figure 10 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.
In order to verify the above conclusions, the Strouhal 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 the Strouhal Number results 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.
To show the difference more clearly of the lock-in region of square column VIVPEH under different attack angles, we choose the dimensionless frequency  squ / squ to indicate the frequency characteristics of square column VIVPEH, where  squ 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;  squ 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 ≤  squ ≤ 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 ≤  squ ≤ 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 ≤  squ ≤ 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 ≤  squ ≤ 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." 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 ∘ .column VIVPEH under different attack angles are given in Figures 12-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.

Phase Angle Analysis of
Figure 16 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  1 and  2 move forward to a distance of  1 when the wave crests of the two adjacent steps appear.At the moment,  1 obviously becomes thinner and longer, while the vibration of the square column is still in the wave crest.Similarly, the vortices of  3 and  4 move forward to a distance of  2 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.

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 Figures 17-20, in which  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.

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  = 1 × 10 6 Ω, 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 Figure 21.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).

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.

Figure 1 :
Figure 1: Physical model of square column energy harvesting system.

Figure 2 :
Figure 2: Computational grid and boundary conditions.

Figure 3 :
Figure 3: Real and imaginary parts of the circuit conjugate solutions.

Figure 4 :
Figure 4: Dimensionless vibration amplitude of square column VIVPEH under different attack angles and different flow velocities.
Square Column VIVPEH under Different Attack Angles.The analysis results of the initial phase angle and the synchronous phase angle of square FFT analysis for  squ = 4.951

Figure 21 : 5  = 15 ∘Figure 22 :
Figure 21: Comparison of maximum value of voltage output and synchronization effective area between VIVPEH with different shapes.

Table 1 :
Calculation parameters for square column VIVPEH system.

Table 1
(16)formula(16), we can calculate the numerical values of  Nor 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.The computational grid and boundary conditions of square column vortex-induced vibration under different attack angles are given in Figure

Table 2 :
Run case of square column in OpenFOAM.

Table 3 :
Comparison of Strouhal Number with different attack angle.