Study on the Efficiency and Dynamic Characteristics of an Energy Harvester Based on Flexible Structure Galloping

Abstract

According to the engineering phenomenon of the galloping of ice-coated transmission lines at certain wind speeds, this paper proposes a novel type of energy harvester based on the galloping of a flexible structure. It uses the tension generated by the galloping structure to cause periodic strain on the piezoelectric cantilever beam, which is highly efficient for converting wind energy into electricity. On this basis, a physical model of fluid–structure interaction is established, and the Reynolds-averaged Navier–Stokes equation and SST K -ω turbulent model based on ANSYS Fluent are used to carry out a two-dimensional steady computational fluid dynamics (CFD) numerical simulation. First, the CFD technology under different grid densities and time steps is verified. CFD numerical simulation technology is used to simulate the physical model of the energy harvester, and the effect of wind speed on the lateral displacement and aerodynamic force of the flexible structure is analyzed. In addition, this paper also carries out a parameterized study on the influence of the harvester’s behavior, through the wind tunnel test, focusing on the voltage and electric power output efficiency. The harvester has a maximum output power of 119.7 μW/mm³ at the optimal resistance value of 200 KΩ at a wind speed of 10 m/s. The research results provide certain guidance for the design of a high-efficiency harvester with a square aerodynamic shape and a flexible bluff body.

1. Introduction

In recent years, the environmental problems caused by the burning of fossil fuels such as petroleum and coal have become more and more serious, and fossil energy is increasingly depleted. It is an important research focus to obtain clean and renewable energy from the environment [1]. Previous studies have shown that vortices may be generated alternately from the two side surfaces of a bluff body that is immersed in the flow, which results in the phenomenon of flow-induced motion (FIM) [2]. Although FIM may endanger the safety of structures, it can be potentially exploited for collecting energy from the environment [3]. To this end, various energy harvesters and technologies have been developed. It is also expected that such environmental energy harvesters can be utilized in practices to power micro-electromechanical systems (MEMS) and wireless sensor systems so that a more convenient realization of structural health monitoring, industry sense and detection, military track, and environmental monitoring [4] can be achieved.

In reference to wind energy harvesting techniques, the majority of related harvesters were developed using the principles of vortex-induced vibration (VIV), galloping, flutter, and buffeting, the main forms of wind-induced induced VIV and galloping. For the VIV energy harvester, the flow rate range of high-efficiency power generation requires that the vortex shedding frequency of the harvester is consistent with the natural frequency. Therefore, it is commonly used for fluid energy harvesting. Williamson’s team [5,6] conducted a lot of research on cylindrical VIV, identified various branches of the VIV amplitude, and summarized the vortex shedding mode into S and P modes. Ding et al. [7] conducted numerical simulation of VIV for different bluff bodies and studied their energy harvesting efficiency. An et al. [8] analyzed and studied the effect of CFD technology on the plate length of the VIV energy harvester on the flow velocity, dynamics, and performance of the wake structure. Zhang et al. [9] used the two-dimensional Reynolds number method to study the VIV of four staggered cylinders, and conducted a series of numerical simulations for energy harvesting. Akaydin et al. [10] developed a cylindrical bluff body VIV piezoelectric energy harvester. When the flow speed is 1.19 m/s, the harvester with a resonance frequency of 3.14 Hz can output a maximum power of 0.1 mW to an optimized load of 2.46 MΩ. Galloping is a typical self-excited vibration phenomenon caused by aeroelastic instability. It mostly occurs in rectangular, angular, and flexible structures, and is usually characterized by low-frequency and high-amplitude oscillations [11]. Due to the greater vibration and higher output power, this aerodynamic instability may be more suitable for energy harvesting than VIV [12]. Barrero-Gil et al. [13] were first to theoretically analyze the potential of using a single-degree-of-freedom (SDOF) system to harvest energy using lateral gallop. Javed et al. [14] used a distributed parameter pattern to study the influence of various aerodynamic force expressions on galloping. Zhao et al. [15] studied the influence of bluff wind exposure area, load resistance, mass of bluff, and piezoelectric sheet length on the output power of a galloping energy harvester. On the other hand, Hu et al. [16] examined the influence of aerodynamic configuration on wind harvesting performance, and found that the VIV of a cylinder could be transformed into galloping if the cylinder was treated via corner modification techniques. Additionally, Sirohi et al. [17] proposed a harvester based on triangular section rods attached to the cantilever beam.

The results from previous studies have shown that transmission lines can be covered with ice on cold days, and their cross section may change to a non-circular shape [18]. Under certain wind speed and wind attack angles, the incoming flow on both sides of the bluff body can produce vortices and shed backwards, and generate an alternate aerodynamic load, which results in the galloping of the transmission line [19]. Previous studies also showed that the tension of the wire can influence galloping significantly, and greater tensions tends to favor the occurrence of galloping. Meanwhile, many galloping energy harvesters were developed by using columns with a square section, as prisms with a square section have more complex cross-section geometric characteristics compared to cylinders. Since Den Hartog first studied and explained the galloping phenomenon, numerous studies have shown that galloping can be widely observed on bluff bodies with a square section [20,21]. Therefore, the square section is usually preferred for the study of galloping energy harvesters.

In this study, a square section energy harvester based on the galloping principle of an iced transmission line was developed. The performance of the harvester was examined via both experimental tests and CFD simulations. The CFD technique was utilized, since it provides a powerful tool to explore the characteristics of flow motions, and to further understand the working mechanism of the harvester. The remainder of the article is organized as follows: Section 2 introduces the design and modeling of the harvester, Section 3 and Section 4 detail the CFD method, and the experimental method, respectively. The specific results are presented and discussed in Section 5, and the main findings and conclusions are summarized in Section 6.

2. Design and Modeling of the Energy Harvester

2.1. Configuration

Figure 1 depicts the setup of the developed harvester. It consists of a flexible structure (i.e., the counterpart of transmission lines with a length of 0.42 m) that aims to receive wind energy from incoming flows, a supporting frame where the flexible structure is fixed, and a wind direction regulator that is mainly composed of a spindle and a deflector. The two ends of the flexible structure are respectively connected with the rigid body of the supporting frame, and the free end of a cantilever beam where a piezoelectric film is adhered. The cantilever beam is made of 0.3 mm × 150 mm stainless steel. The piezoelectric film, i.e., the PVDF (Polyvinylidene Fluoride) film (IPS-17020, ZHINK TECHNOLOGY) has dimensions of 30 mm (length) × 12.1 mm (width) × 0.28 μm (depth), and the density is 1.78 × 103 kg/m³. The spindle is fixed with the supporting frame and the deflector rigidly. Due to the guiding effects of the deflector, incoming wind flows tend to blow in a direction that is perpendicular to the windward surface of the flexible structure. When the wind speed exceeds a certain value, the flexible structure can vibrate up and down severely, i.e., galloping occurs. The vibrating flexible structure then drives the cantilever and therefore the PVDF film to sway, which fulfills the conversion of wind energy into electricity.

When the energy harvester faces the incoming wind, the flexible bluff body swings around the plane and causes the thin wire to generate tension. The tension at one end of the thin wire directly drives the cantilever beam to produce periodic alternating strain. At the same time, the piezoelectric film is deformed. A certain charge is generated, which effectively converts the flow energy of the air into electric energy. The energy harvester has a simple structure and low manufacturing costs, and is very suitable for the surrounding windy environment all year round. Compared with various forms of energy harvesters, this paper’s harvester is not affected by wind direction and is suitable for any wind direction angle. Since the designed energy harvester has excellent parameters, and the critical wind speed generated is small, the flexible structure is prone to galloping and improves the efficiency of energy harvesting.

2.2. CFD Mathematical Model

As shown in Figure 2, the harvester system can be mathematically simplified as a SDOF system. Taking the middle section of the flexible structure as the simulation research object, it is assumed that the wind field is incompressible and the effect of structural torsion and the inverse piezoelectric effect in the system can be reasonably neglected.

Thus, the dynamics of the harvester body over a unit length can be depicted by [22]:

$$M\ddot{y}+C\dot{y}+Ky=F_{\mathrm{y}}$$

(1)

where y denotes the lateral displacement of the square column, M is the mass per length, F_y is the force acting on the square structure in y direction, C is the mechanical damping coefficient, and K is the stiffness of the vibration model.

The details of the equation can be referred to in Barreiro-Gill et al. [23].

According to the quasi-steady theory [24], wind force can be expressed as function of wind velocity U, the characteristic length of the square section D (or the side length in this study) and a dimensionless coefficient, or the lift coefficient $C_L$ in this study:

$$F_y=\frac{1}{2}\rho U^{2}D{C}_L$$

(2)

The main parameters involved in this model of harvester are listed in

Table 1. Main physical model parameters of the harvester.

Item	Symbol	Value
Mass(g)	M	0.144
Diameter(mm)	D	30
Stiffness (N/m)	K	30.2
Damping ratio	$\zeta$	0.0018
Air density(Kg/m³)	$\rho$	1.225
Motion viscosity (m2·s−1)	v	1.41 × 10−6

, where $\zeta$ represents the damping ratio, which is defined as $\zeta =C/2\sqrt{MK}$. Note that the values of damping ratio and effective stiffness coefficient can be determined through experiment.

2.3. Piezoelectric Performance Equation

The energy harvester uses a d₃₁ piezoelectric sheet, i.e., the polarization direction is perpendicular to the applied stress direction, and the working mode of the piezoelectric material is the d₃₁ mode with the polarized direction parallel to the surface normal. To simplify the analysis, a load resistance (R_L) was attached to the conductive electrode of the PVDF piezoelectric ceramic film. The piezoelectric constitutive relationship can be directly expressed as [25]:

$$D_3=d_{31}{T}_1+{\epsilon }_{33}^T{E}_3$$

(3)

where d₃₁ is the piezoelectric charge constant, D₃ is the electric displacement component, T₁ is the stress component, E₃ is the electric field component, and ${\epsilon }_{33}^T$ is the dielectric constant.

3. CFD Simulation Method

3.1. Numerical Methods

The fluid–structure interaction solution method in conjunction with the weak coupling solution technique were used to numerically calculate the evolutions of wind field over the computational domain and to obtain the aerodynamic force acting on the harvester. The two-dimensional Reynolds number average Navier–Stokes (RANS) equation and the k-ω SST (shear stress transport) turbulent model [26] were used to numerically simulate the flow around a square cross-section. This turbulence model was chosen because, compared to other turbulence models, the k-ω model of SST is a hybrid of the k-ω and k-ε models. It activates the k-ω model near the wall and the k-ε model in free flow. It has better performance when predicting the boundary layer flow of the reverse pressure gradient. Therefore, this paper adopts the SST k-ω turbulence model in CFD simulation.

The interaction effects between the wind flows and the harvester structure were embedded into the CFD to implement the link with the main program. At the same time, the structure dynamic equation was solved and the structure domain calculated, and then the motion lift coefficient and displacement of the structure were obtained. At the same time, the overset grid technology was used to update the grid of the computational domain in real time.

User-defined functions (UDFs) were compiled (in C Language) and embedded into the CFD software to realize the simulation of fluid–structure interaction. In this paper, DEFINE_CG_MOTION was used as the motion function to define the motion mode. The overset grid model was used to define computational grid movement and data exchange. The Newmark-β method was used in UDF (refer to the Appendix A for details). This method is a numerical integration method for solving differential equations. It has been widely used to solve problems with oscillating systems [27].

Based on the basic assumptions of the central difference method, the explicit difference method of dynamic equations was adopted. The basic assumption of the Newmark-β method is that the acceleration changes linearly in each time increment, and the characteristics of the system remain constant during this interval. The equilibrium requirement of the force system acting on the mass at any time is as follows: the Newmark-β method approximates the speed and displacement of the system at time ($t+\Delta t$) through the following two assumptions:

$$\dot{y}(t+\Delta t)=\dot{y}(t)+(1-\gamma )\Delta t\ddot{y}(t)+\gamma \Delta t\ddot{y}(t+\Delta t)$$

(4)

$$y(t+\Delta t)=y(t)+\Delta t\dot{y}(t)+(\frac{1}{2}-\beta )\Delta t^2\ddot{y}(t)+\beta \Delta t^2\ddot{y}(t+\Delta t)$$

(5)

The two parameters $\gamma$ and $\beta$ in the formula can be selected as required, and different combinations correspond to different processing effects. When $\gamma =\frac{1}{2}$ and $\beta =\frac{1}{2}$, it is the central difference method [28]. Therefore, the acceleration and velocity can be obtained.

$$\ddot{y}(t+\Delta t)=\frac{1}{\beta \Delta t^2}(y(t+\Delta t)-y(t))-\frac{1}{\beta \Delta t}\dot{y}(t)-(\frac{1}{2\beta }-1)\ddot{y}(t)$$

(6)

$$\dot{y}(t+\Delta t)=\frac{\gamma }{\beta \Delta t}(y(t+\Delta t)-y(t))+(1-\frac{\gamma }{\beta })\dot{y}(t)+(1-\frac{\gamma }{2\beta })\Delta t\ddot{y}(t)$$

(7)

The motion of the system at time $t+\Delta t$ can be expressed as:

$$M\ddot{y}(t+\Delta t)+C\dot{y}(t+\Delta t)+Ky(t+\Delta t)=F_y(t+\Delta t)$$

(8)

Integrating the above formulas leads to:

$$\tilde{K}y(t+\Delta t)={\tilde{F}}_y(t+\Delta t)$$

(9)

where

$$\tilde{K}=\frac{1}{\beta \Delta t^2}M+\frac{\gamma }{\beta \Delta t}C+K$$

(10)

$$\begin{array}{l}\tilde{F}(t+\Delta t)=F(t+\Delta t)+M\left[\frac{1}{\beta \Delta t^2}y(t)+\frac{1}{\beta \Delta t}\dot{y}(t)+(\frac{1}{2\beta }-1)\ddot{y}(t)\right]\\ +C\left[\frac{\gamma }{\beta \Delta t}y(t)+(\frac{\gamma }{\beta }-1)\dot{y}(t)+(\frac{\gamma }{2\beta }-1)\Delta t\ddot{y}(t)\right]\end{array}$$

(11)

The numerical simulation flowchart is shown in Figure 3.

3.2. Computational Fluid Domain and Grid

As shown in Figure 4a, the simulation used a rectangular computational fluid domain. The calculation domain had a width of 20D and a length of 30D. The middle position of the square was 10D away from the inlet. The boundary conditions used included the left speed inlet, the right pressure outlet, and the upper and lower sides’ sliding wall surface. The center of the square column was the origin of the coordinate system; x and y respectively represent the along-wind and crosswind directions.

The model establishment and grid division were all carried out in the commercial software Ansys ICEM. The grids under different working conditions were all divided by non-uniform structured grids to ensure the accuracy of calculation and save a lot of computing resources. In computational fluid dynamics, the dynamic overset grid method has efficient dynamic grid processing capability and can guarantee the quality of the grid, so it has been widely used in unsteady flow simulation. Therefore, this research used overset grid technology; many concepts related to nested grid computing used today can be traced back to the breakthrough idea of Joseph Steger [29,30]. The overset grid method was used in Fluent; when FIM occurred, the constructed grid oscillated with the square column to ensure that the grid was not updated, which effectively avoided the negative element grid. The square coordinate grid (called the component grid) matched the square cylinder, and the background grid adopted a unified Cartesian grid. Details of the overset grid are shown in Figure 4b.

In order to make the grid height Y+ of the first layer of the wall satisfy the wall function, the height value was determined according to the required Reynolds number [31]. Enough densification was performed close to the wall, and at least 20 layers of wall surface grids were set with a rate of change of 1.2. The Y+ of the wall grids were all around 1.

3.3. Solving Method

The turbulence, pressure, momentum, and other fluid control equations used the second-order upwind spatial discretization algorithm, and the gradient calculation used an algorithm based on a least squares unit [32]. The discretization was performed using the second-order upwind difference scheme, and the transient term was discretized using the second-order backward Euler scheme. The first-order algorithm was used for preliminary calculation and then the second-order algorithm was used to improve the calculation accuracy. All simulation conditions adopted the pressure-velocity (SIMPLEC) coupling method, and the convergence was determined by the size of the residual error of the governing equation [33]. The continuity equation, momentum equation, k equation, and ω of all scale residuals below 1 × 10⁻⁵ were used. The equation was used as the convergence criterion, and the corresponding calculation results were observed to be stable, which can be considered reliable results.

3.4. Validation of Numerical Model

To ensure reliability, and that the model used numerical simulation and obtained a convergent and accurate solution within a reasonable calculation time, this paper considers different grid densities and time steps.

Table 2. Verification of grid density independence (Re = 2.2 × 10⁴).

Grid Density	CD	CL	St
Coarse (32,414)	2.20	1.44	0.128
Medium (43,954)	2.25 (2.27%)	1.41 (2.08%)	0.127 (0.78%)
Fine (96,153)	2.26 (2.65%)	1.45 (0.90%)	0.127 (0.78%)

and

Table 3. Validation of time-step independence (Re = 2.2 × 10⁴).

Time Steps	CD	CL	St
0.001	2.25	1.45	0.127
0.0008	2.30 (2.22%)	1.56 (7.59%)	0.129 (1.57%)
0.0006	2.29 (1.78%)	1.66 (14.5%)	0.132 (3.93%)
LES [31]	2.21 (1.78%)	1.71 (17.9%)	0.128 (0.79%)
DNS [32]	2.18 (3.11%)	1.51 (4.14%)	0.132 (3.93%)

list the results of different grid parameters and time steps, in which C_D denotes the mean drag coefficient over the time domain, C_L stands for the root mean square (RMS) of the lift coefficient, and St represents the Strouhal number.

As shown in Table 2, the results under the three different grid densities were similar. The scheme with medium grid resolution was adopted in this study so that a balance between computational efficiency and cost could be achieved.

Meanwhile, the transient numerical method was adopted and time-step independent verification was carried out to evaluate the uncertainty of the value. The independence of the time steps was verified by using different time step sizes of 0.0006 s~0.001 s/step for numerical simulation. According to the calculation results for different time steps, a small time step cannot improve the simulation accuracy but takes a long time. Therefore, the time step of 0.0008 s per step was selected. The calculated results of the square at Re = 2.2 × 10⁴ were similar to the results presented by Cao [34] and Trias [35], indicating that the turbulence model adopted can obtain more reliable results.

4. Experimental Method

The experiment was carried out in a subsonic closed recirculation wind tunnel. The design range of the flow velocity at the working section of the tunnel was 1–22 m/s, with the turbulence intensity not exceeding 0.5%. The layout of the experiment is shown in Figure 5. When exploring the influence of the external resistance R_L, the external load resistance of 10 kΩ–1 MΩ was connected in sequence, and the wind tunnel control speed changed in a gradient of 3 m/s–10 m/s. The output voltage was collected and stored by a data acquisition instrument (NI 9234). The highest sampling frequency was 2840 Hz for acquisition, and the acquisition time exceeded 120 s. The voltage and displacement were inputted into labview for visual observation and analysis. The displacement of the transmission-line-like flexible structure was monitored by a laser displacement sensor that was installed on the side wall of the wind tunnel. The output voltage of the harvester was recorded when the vibration of the flexible structure became stable.

5. Results and Discussion

5.1. Performance Analysis of the Harvester

On obtaining the voltage time-history data through wind tunnel experiments, the root mean square (RMS) value of voltage V_RMS can be calculated:

$$V_{RMS}=\sqrt{\frac{1}{T_2-T_1}{\displaystyle {\int }_{T_1}^{T_2}u(t)}dt}$$

(12)

Accordingly, the mean power P can be obtained as:

$$P=\frac{1}{T_2-T_1}{\displaystyle {\int }_{T_1}^{T_2}\frac{u{(t)}^2}{R_L}}dt$$

(13)

in which T₂ and T₁ are the start and end time during the acquisition period, respectively, $R_L$ is the external load resistance, and u(t) is the instantaneous voltage.

Figure 6 displays the results of the RMS voltage and mean power output by connecting various load resistances from 10 kΩ to 1 MΩ in the PVDF film in the wind speed range of 3–10 m/s. It can be seen in Figure 6a that under the same load resistance, the wind speed continues to increase, and the output RMS voltage has an obvious upward trend. The minimum RMS voltage recorded under a load resistance of 10 kΩ is 350 mV. In addition, the output RMS voltage of this resistance model is relatively small. In the wind speed range of 3–10 m/s, the output RMS voltage steadily increased from 350 mv to 0.14 V. The highest RMS voltage recorded under the load resistance of 1 MΩ was 2.96 V. It can be found from Figure 6b that in the wind speed range of 3–10 m/s, the rising trend of output power was similar to the rising trend of RMS output voltage. When the wind speed was 10 m/s, the maximum load resistance of 200 kΩ in this experiment was 12.1 μW. The results indicate that when the load resistance is constant, as the input air speed increases, the voltage and power between the load resistances increase almost in advance.

According to the RMS voltage and output power calculated by Formulas (12) and (13), the power density and voltage density per scanning volume, which are defined as the power value and RMS voltage value divided by the volume of piezoelectric film, were calculated. In the analysis that followed, the load resistance first had to be matched to evaluate the harvesting performance. Figure 7 displays the change of output voltage density and power density with load resistance under different wind speeds. Figure 7a shows that the overall output trend of the measured voltage density increases with the increase of resistance within a certain wind speed range, and the maximum voltage density is 29.1 V/mm³. Figure 7b shows that within a certain range, the overall output power trend of the measured power density at all wind speeds increases as the resistance increases. The optimal value of the load resistance for different wind speeds is between 100 kΩ and 430 kΩ. Outside this load resistance range, the average power generated is significantly reduced. After intensive testing, it was found that when the load resistance is about 200 kΩ, it provides the highest output average power. Therefore, R₀ = R_L = 200 kΩ is the optimal electrical load resistance, and this result will be used for the efficiency of the energy harvester in the subsequent experimental discussion when the load resistance is 200 KΩ, and the power density is 119.7 μW/mm³. Previous studies have pointed out that there is an optimal resistance R₀ for the best collection efficiency. Neglecting the effects of damping and dielectric loss, the optimization resistance can be determined as [36]:

$$R_0=\frac{1}{2\pi fC}$$

(14)

where C is the harvester capacitance and f is the vibration frequency.

To further explore the working performance of the harvester, it was necessary to analyze the piezoelectric film dynamics characteristics. Figure 8 displays the time history of the output voltage with a load resistance of 300 kΩ at a wind speed of 3–7 m/s. As demonstrated, the voltage signal fluctuates periodically for each case associated with a fixed wind speed, and the amplitude of the voltage increases distinctly as the wind speed increases (from 0.14 V to 1.12 V). Apparently, the periodic variation of output voltage is attributed to the periodic evolution of the strain around the fix end of the cantilever beam. Because the galloping amplitude increases consistently with the increase of wind speed, the strain becomes larger for the cases with stronger incoming wind, which further leads to higher output voltage.

5.2. Dynamic Characteristics Analysis

By calculating the damping coefficient, mass ratio, and Reynolds number selected in this study, it can be found that the value of U_g/U_v is much less than 1, which means that the square column FIM has a strong interaction between VIV and galloping [37]. It can be considered that galloping continuously affects the amplitude and locking area of VIV, and VIV also changes the form of galloping. In order to better understand FIM and flow structure, this paper introduces four kinds representative simulation examples under various Reynolds numbers (Re = 2054/6161/16,430/20,537). This includes four FIM regions: VIV initial branch, VIV upper branch, VIV-galloping transition, and galloping. These vortex patterns, at different Reynolds numbers, are analyzed. In addition, the effects of the amplitude and frequency responses for square columns are discussed as well.

This paper presents the instantaneous vorticity diagram, velocity cloud diagram, and velocity streamline of the square column at different wind speeds (Reynolds number), including the vortex shedding at the highest and lowest points of the structure. In those figures of the vortex pattern, the red vorticity is positive, which means counterclockwise rotation; the blue vorticity is negative, which means clockwise rotation. There are two vortex shedding modes: 2S mode (in one cycle, two vortices shed from opposite sides of the cylinder) and 2P mode (in each cycle, two pairs of vortices shed from each side of the cylinder). As the wind speed changes, the vortex shape can change, and the branch transitions affect each other [27]. In this paper, the vortex shedding mode is studied by observing the instantaneous vorticity contours around the square column oscillator, where the vorticity is defined as $\omega =\partial v/\partial x-\partial u/\partial y$.

When the wind speed is 1 m/s (Re = 2054), the harvester is in the initial branch of VIV. As shown in Figure 9, the 2S model of the vortex shedding mode in the wake can be clearly observed; that is, the positive vortex and the negative vortex separate out during the vibration period. This mode is the classical von Kármán Street. In this range, as Re increases, the size of the vortex is larger than when Re is low, which is consistent with the numerical results given by Ding et al. [27].

As the wind speed increases to 3 m/s (Re = 6161), the number of vortices shed in each oscillation cycle increases. In each oscillation period, 6 vortices shed from the square column and shed in the 2P + 2S mode of two pairs of vortices and one single vortex, as shown in Figure 10. This vortex pattern is known as Quasi-2P, which means two pairs of vortices shedding per cycle.

When the wind speed reaches 8 m/s (Re = 16,430), the vortex pattern in the VIV-galloping transition zone is formed. As shown in Figure 11, ten vortices shed from the square column in an oscillation cycle. According to previous research by Ding [38], a similar near-wake vortex structure was captured as 4P + 2S. There are still many controversies about the identification method of multiple vortices shedding from square columns. At the highest and lowest points of vibration, a single vortex appears first. In the next downward or upward process, a total of two pairs of vortices appear.

When the wind speed is 10 m/s (Re = 20,537), the FIM of the square column is in the fully developed galloping state zone. As shown in Figure 12, there is a fundamental difference between galloping vibration and VIV. When entering the galloping state, a total of 14 vortices will shed in the manner of 6P + 2S in each oscillation period, two of which are separate vortices. According to the results of previous studies [39], it can be found that the vortex of the square column is at high Reynolds number. The type is np + 2S, where n represents the number of p and increases with the increase of Reynolds number or flow rate. The 2S vortex appears at the highest point and the lowest point of the displacement and falls off respectively.

Within a certain wind speed range (as shown in Figure 13, Figure 14, Figure 15 and Figure 16), the lift coefficient decreases with the increase of wind speed. When the wind speed is relatively small, the lift coefficient becomes a regular oscillation under low Reynolds number [40]. In the initial stage of VIV, comparing the lift coefficient and the displacement graph, the lift coefficient and the displacement performance are highly consistent. The displacement curve and lift coefficient curve are transformed by FFT to obtain the frequency spectrum. In Figure 15b, it can be seen that both the displacement spectrum and lift spectrum have multiple peak frequencies and that they maintain good consistency and their main frequency is close to the natural frequency. When galloping occurs, the dominant frequency of the displacement spectrum decreases, and the lift spectrum frequency is more obviously controlled by the vortex deflation frequency. In Figure 16b, when the square column enters the state of full galloping, it can be found that the proportion of low frequency in the displacement spectrum increases. The increase in its proportion is also an important signal for the full development of galloping. When galloping occurs, the frequency response of the square column is obviously different; its vibration is hardly controlled by the vortex breakaway frequency and the main frequency is always lower than the natural frequency, which also leads to the occurrence of large vibrations.

Through the numerical simulation results, as the wind speed increases (the Reynolds number also increases), although the lift coefficient decreases, the lift applied to the two sides of the bluff body continues to increase. The maximum displacement of the system also increases, which causes an increase in the instantaneous stress of the piezoelectric thin film. Combined with the results of the output power of the energy harvester, it can be found that when the Reynolds number is low, the bluff body’s vibration amplitude and frequency are low, so the harvesting efficiency is at a low level. However, when the incoming wind speed increases and begins to enter the galloping state, the coupling effect of VIV and galloping can accelerate the onset of galloping and increase the displacement response of the bluff body at low and medium wind speeds. After entering full galloping, the amplitude of the entire system will increase significantly, and the output power will be increased. In addition, compared with VIV, the vibration of the square column in the galloping process is unstable, and its amplitude will increase with the increase of the incoming flow velocity. This means that the kinetic energy of the fluid captured by the square column FIM harvesting system will continue to increase with the increase of the incoming wind, and there is no limit [41]. Therefore, it can be seen that the square cylinder has obvious advantages in the application of FIM energy conversion, especially under high Reynolds number. The form of wind-induced vibration is closely related to the damping, elastic stiffness, and mass of energy harvester, which is used as the power source of energy collection. The harvesting efficiency of the flexible structure can be further improved only through analysis of its dynamic characteristics and reasonable design.

6. Conclusions

In this study, a novel type of piezoelectric energy harvester was developed based on the principle of wind-induced galloping. The dynamic model of the fluid–structure interaction system was established, and CFD techniques in conjunction with the Newmark-β method were adopted to obtain corresponding numerical solutions. Wind tunnel tests were also conducted to detail the working performance of the harvester. The main findings and conclusions are summarized below:

• The flexible structure experienced four FIM stages. The wind speed increased to 8 m/s in the transition zone between vortex-induced resonance and galloping. When the wind speed reached 10 m/s, it started to gallop completely, and the harvester system experienced low-frequency and high-amplitude vibration.
• Within the studied wind speed range of 3–10 m/s, it was found that as the wind speed increased, the vortex-induced shedding frequency of the flexible structure first increased and then decreased, but the associated displacement increased consistently.
• When the resistance was constant, the output voltage and power of the harvester increased with the increase of wind speed. When the wind speed was constant, the output voltage increased consistently with increasing resistance. By contrast, the output electric power first increased and then decreased with the resistance, and there was optimal resistance at around 100 kΩ. The maximum output power of the energy harvester was 119.7 μW/mm³.