Optimization of Galloping Piezoelectric Energy Harvester with V-Shaped Groove in Low Wind Speed

Abstract

A square cylinder with a V-shaped groove on the windward side in the piezoelectric cantilever flow-induced vibration energy harvester (FIVEH) is presented to improve the output power of the energy harvester and reduce the critical velocity of the system, aiming at the self-powered supply of low energy consumption devices in the natural environment with low wind speed. Seven groups of galloping piezoelectric energy harvesters (GPEHs) were designed and tested in a wind tunnel by gradually changing the angle of two symmetrical sharp angles of the V-groove. The GPEH with a sharp angle of 45° was selected as the optimal energy harvester. Its output power was 61% more than the GPEH without the V-shaped groove. The more accurate mathematical model was made by using the sparse identification method to calculate the empirical parameters of fluid based on the experimental data and the theoretical model. The critical velocity of the galloping system was calculated by analyzing the local Hopf bifurcation of the model. The minimum critical velocity was 2.53 m/s smaller than the maximum critical velocity at 4.69 m/s. These results make the GPEH with a V-shaped groove (GPEH-V) more suitable to harvest wind energy efficiently in a low wind speed environment.

1. Introduction

Energy harvesting is usually the process of converting waste or useless energy into electrical energy in the external environment, providing electrical energy for low-power electronic components. The ultimate goal is to develop self-powered sensors, actuators, and other electronic devices. The research on energy harvesters has always been one of the hot topics in the literature, and a large number of researchers have made extensive research in this field. The energy harvester based on vibration has attracted widespread attention because mechanical vibration can make energy in many places where thermal or light energy is not suitable [1]. For example, it can provide power for electromechanical systems and actuators, and can also be applied to environmental monitoring equipment, medical transplantation, and wireless sensors. In addition, due to the limited battery life and the cost and time required for maintenance, the vibration energy harvester (VEH) can also be employed to replace some small batteries. Generally speaking, the VEH mainly includes piezoelectric, electromagnetic, and electrostatic types, among which the piezoelectric energy harvester (PEH) has outstanding advantages such as high force electric coupling effect and energy density, no electromagnetic interference, and easy miniaturization process. These advantages make PEH ideal for low-power wireless sensor nodes [2].

PEH uses piezoelectric effect to convert mechanical energy into electrical energy. The external environment vibration causes mechanical deformation of the piezoelectric body, and polarization occurs inside it. Then positive charges and negatives appear on a pair of electrode faces of the piezoelectric wafer. The voltage can be obtained through an external resistance. The power of energy harvesting depends primarily on the form and amount of vibrational energy, as well as the performance of the harvester. In previous studies, VEH mainly focused on harvesting energy from foundation vibration, because it is a simple and common vibration mode [3]. A significant number of papers were issued to discuss the scheme of improving the level of energy harvester based on different mechanisms and techniques. [4]. Most researchers increase the efficiency of the PEH by using a nonlinear electronic interface circuit [5]. Wang et al. [6] proposed a model and experiment of a broadband piezoelectric vehicle with four-stable caused by the contact nonlinearity and magnetoelasticity interaction of the cantilever beam surface. Wang et al. [7] considered a tri-stable PEH with two external magnets.

However, with the development of energy acquisition, flow-induced vibration energy harvester (FIVEH) has attracted increasing attention over the last few years [8,9]. There are three main types of flow-induced vibrations (FIV) that cause periodic vibration of the energy harvester: vortex-induced vibration (VIV), flutter, and galloping. In the natural environment, VIV is a relatively common vibration mode in FIV. When the vortex falls off at a frequency close to one of the natural frequencies of the vibration system, the system will lock or synchronize, and then a large vibration will occur due to the resonance phenomenon. The vibration gets the characteristics of stability, self-limitation, and large amplitude, and based on the above characteristics, it is employed in the research of energy acquisition. A series of scientific research have been issued to probe the performance of different types of vortex-induced vibration energy harvesters (VIVEHs) [10,11,12,13]. Nevertheless, due to the existence of the frequency lock-in phenomenon, the energy harvester can only produce large vibration in a certain range of velocity [14]. The energy harvester of the wing-based plate based on the flutter principle was designed to work normally over a wide range of velocity. When the flow velocity increases from zero, such a fluid–solid coupling system may have a critical velocity, commonly referred to as the flutter velocity. At this flutter velocity, self-excited motion begins to occur because the absolute value of the structural damping is exactly equal to the absolute value of the damping caused by the aerodynamic effect [15]. As there is not any frequency lock-in phenomenon in flutter, the energy harvester can work normally as long as the flow velocity exceeds the critical velocity [16,17,18]. It is characterized by high frequency and small amplitude, because flutter belongs to multi-modal coupling vibration [19], which is not conducive to the transduction of energy by piezoelectric materials.

Based on the deficiency of the above two kinds of FIVEHs, the method and concept of a galloping energy harvester (GEH) are proposed [20,21,22]. Therefore, in order to make full use of the FIVEH, a lot of researches were carried out, which was on the basis of the cylindrical FIVEH by changing the cross-sectional shape of the cylinder to stimulate the galloping motion of the system. During galloping, transverse self-excited vibration is generated by a structure in a stable flow field perpendicular to the direction of flow [23]. Some structures with non-streamlined or asymmetrical cross sections (rectangular, D-shaped, hexagonal, triangular, etc.) may undergo galloping at a certain inflow velocity and direction [24,25]. The galloping has a larger amplitude and a wider vibration frequency band than the VIV, and is generally considered to be a vibration of a single mode, which provides a more effective guarantee for stable acquisition of energy [26]. The structure of the circular cross section does not gallop, because for a circular cross-section, if no vortex shedding occurs, no lift will occur [27]. When the cylinder is vibrating perpendicular to the direction of the fluid, the force of the fluid is consistent with the direction of the relative flow velocity, and the component in the direction of the structural vibration is opposite to the direction of motion of the structure. This is a positive damping force and does not generate self-excitation vibration. Therefore, in order to make full use of the FIVEH, a lot of research has been carried out, which is on the basis of the columnar FIVEH changing the cross-sectional shape of the cylinder to stimulate the galloping motion of the system [28,29,30]. Barrero-Gil et al. [31] investigated the influence of section geometry on transverse galloping of a single degree of freedom system. Liu et al. [32] designed a galloping piezoelectric energy harvester (GPEH) with a Y-type three-blade bluff body. Amin et al. [33] studied the characteristics of the GPEH with a square cylinder under combined base excitations and galloping. Javed et al. [34,35] investigated the influence of different aerodynamic load expressions on galloping response in square cylinder GPEHs. Abdelmoula et al. [36] investigated the performance of a square cylinder GPEH when considering complex circuit systems.

To sum up, there have been numerous reports on the research of FIVEHs in academic circles, and many achievements have been achieved. Although most researchers only study the vibration characteristics of some basic shape section columns (such as circle, triangle, square, rectangle, and trapezoid). There is little research on the special section shape. In previous research, an acceleration sensor was installed on the square cylinder to collect the fluid force directly. This method is simple in principle and convenience in operation. However, in this study, due to the relatively small mass of the cylinder, the mass of the acceleration sensor may not be ignored, and the error of the results may not be predicted. Therefore, in this paper, the laser displacement sensor is used to evaluate the displacement signal of the cylinder. If the method of direct derivation of the data is used to obtain the fluid force, a large error will occur. Therefore, based on the system identification theory, combined with the singular value decomposition (SVD) and sparse identification algorithm (SINDy) [37], we get the empirical coefficient of fluid force. Based on the SVD theory, the noise reduction method has excellent robustness and the SINDy has the ability to make the governing equations sparse in a high-dimensional nonlinear function pace.

Based on the Euler–Lagrange equation and the principle of virtual work, combined with Galerkin discrete method, the control equation of GPEH-V is established, the main factors affecting the flow force are analyzed, and the structure of a V-shape on the windward side of the square cylinder is proposed to improve its amplitude in low velocity fluid. Through the wind tunnel experiment, the characteristics of different models to GPEH are analyzed.

2. Electromechanical Modeling of the Galloping

Inspired by the GPEH [31], we proposed the concept of the GPEH-V. Its schematic diagram shown in Figure 1. It is mainly composed of a piezoelectric cantilever beam and a square cylinder. The piezoelectric cantilever beam is composed of a 0.8 mm thick beryllium bronze support layer and its bonded macro fiber composite (MFC M8507 P2). At the free end of the cantilever beam, a square cylinder with a V-shaped groove on the windward side is connected. When the wind velocity exceeds the critical value, the cylinder will be subjected to the transverse wave force generated by the airflow. The blue part in Figure 1 represents the MFC piezoelectric sheet attached to the cantilever beam, which is connected in series with the resistance through two plane electrodes with negligible thickness. The geometric dimensions and physical parameters of the energy harvester are shown in

Table 1. Physical parameters of the GPEH-V.

Parameters	Description	Values
$L_b$ (m)	Length of beryllium bronze substrate	110 × 10−3
$w_b$ (m)	Width of beryllium bronze substrate	15 × 10−3
$t_b$ (m)	Thickness of beryllium bronze substrate	0.8 × 10−3
$L_p$ (m)	Length of MFC M8507 P2	100 × 10−3
$w_p$ (m)	Width of MFC M8507 P2	0.3 × 10−3
$t_p$ (m)	Thickness of MFC M8507 P2	100 × 10−3
$L_0$ (m)	Length of the tip mass	20 × 10−3
$D$ (m)	Width of the tip mass	8.3 × 10−3
${\rho }_b$ (kg/m3)	Density of beryllium bronze	8.3 × 103
${\rho }_p$ (kg/m3)	Density of MFC M8507 P2	5.4 × 103
${\rho }_0$ (kg/m3)	Density of the tip mass	4.5 × 102
$E^b$ (Pa)	Young’s modulus of beryllium bronze	128 × 109
$E^p$ (Pa)	Young’s modulus of MFC M8507 P2	30.336 × 109
$d_{31}$	Strain coefficient of MFC M8507 P2	−320 × 10−12
${\epsilon }_{33}$	Permittivity component at constant strain	13.28 × 10−9
${\rho }_a$ (kg/m3)	Air density	1.24
R (Ω)	Load resistance	9000

. The shear deformation and moment of inertia of the cantilever can be ignored, because the thickness of the cantilever is very small relative to its length.

Based on the Lagrangian equation and the virtual work principle, combined with the Galerkin discrete method, the governing equations of the GPEH are established. Therefore, expressions of kinetic energy T, total potential energy U, and the non-conservative work W due to fluid force and damping force should be established first. Considering the moment of inertia caused by the lateral movement of the cylinder, the expressions of kinetic energy and potential energy are given as:

$$\begin{array}{l}T=\frac{1}{2}\left[{\displaystyle \underset{V_b}{\int }\rho {\left[\frac{\partial w\left(x,t\right)}{\partial t}\right]}^{2}d{V}_b}+{\displaystyle \underset{V_b}{\int }{\rho }_p{\left[\frac{\partial w\left(x,t\right)}{\partial t}\right]}^{2}d{V}_p}\right]+\\ \begin{array}{cc}& \frac{1}{2}{m}_0{\left[{\frac{\partial w\left(x,t\right)}{\partial t}|}_{x=L_b}+{\frac{D}{2}\frac{{\partial }^{2}w\left(x,t\right)}{\partial x\partial t}|}_{x=L_b}\right]}^2+\frac{1}{2}J{\left[{\frac{{\partial }^{2}w\left(x,t\right)}{\partial x\partial t}|}_{x=L_b}\right]}^2\end{array}\end{array}$$

(1)

$$U=\frac{1}{2}\left[{\displaystyle \underset{V_b}{\int }{\sigma }_x^b{\epsilon }_x^{b}d{V}_b}+{\displaystyle \underset{V_p}{\int }{\sigma }_x^p{\epsilon }_x^{p}d{V}_p}-{\displaystyle \underset{V_b}{\int }{E}_3{D}_{3}d{V}_p}\right]$$

(2)

The strains and stresses in the piezoelectric layer and substrate are expressed as:

$${\epsilon }_x^b={\epsilon }_x^p={\epsilon }_x=-z\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}$$

(3)

$${\sigma }_x^b=E^b{\epsilon }_x$$

(4)

$${\sigma }_x^p=E^p\left({\epsilon }_x-d_{31}{E}_3\right)=E^p{\epsilon }_x-e_{31}{E}_3$$

(5)

where E^b and E^p represent Young’s modulus under a constant electric field, and E₃ is the electric field in the direction of polarization [14]. The relationship between the generated voltage V(t) and the electric field E₃ is E₃ = −V(t)/t_p [38]. Additionally, the expression for D₃ [14] is written as:

$$D_3=d_{31}{E}^p{\epsilon }_x+{\epsilon }_{33}{E}_3=e_{31}{\epsilon }_x+{\epsilon }_{33}{E}_3$$

(6)

The work done by non-conservative forces includes three parts: one is the power generated by the circuit, the other is the work done by fluid force F(t), and the third is the work done by the damping force of structural system. The variation expression of work done by the non-conservative force is as follows:

$$\delta W=-V\delta Q_R+F\left(t\right)\delta \overline{w}\left(t\right)-{\displaystyle \underset{0}{\overset{L_0}{\int }}c\frac{\partial w\left(x,t\right)}{\partial t}\delta w\left(x,t\right)dx}$$

(7)

where

$$\overline{w}\left(t\right)={w\left(x,t\right)|}_{x=L_b}+{\frac{D}{2}\frac{\partial w\left(x,t\right)}{\partial x}|}_{x=L_b}.$$

The external force is expressed as:

$$F\left(t\right)=\frac{1}{2}{\rho }_{0}D{L}_0{U}_0^2{\displaystyle \sum _{j=1}^7{a}_j{\alpha }^j=}\frac{1}{2}{\rho }_{0}D{L}_0{U}_0^2{\displaystyle \sum _{j=1}^n{a}_j{\left[\frac{1}{U_0}\frac{\partial w\left(x,t\right)}{\partial t}+\frac{D}{2U_0}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t\partial x}\right]}^j}$$

(8)

where the selection of a_i value is related to the shape of the cylinder [36].

Using the extended Hamilton principle, the governing equations of motion and boundary conditions of the galloping systems are expressed as:

$$\begin{array}{l}EI\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}+c\frac{\partial w\left(x,t\right)}{\partial t}+meq\frac{{\partial }^2\left(x,t\right)}{\partial t^2}+\theta V\left(t\right)\left(\frac{d\delta \left(x-x_1\right)}{dx}-\frac{d\delta \left(x-x_1-L_p\right)}{dx}\right)=\\ \frac{1}{2}{\rho }_a{D}_c{L}_c{U}_0^2\times {{\displaystyle \sum _{i=1}^3{a}_1\left[\frac{1}{U_0}\frac{\partial w\left(x,t\right)}{\partial t}+\frac{D_c}{U_0}\right]}}^i\left(\delta \left(x-L_b\right)+\frac{D}{2}\frac{d\delta \left(x-L_b\right)}{dx}\right).\end{array}$$

(9)

The associated boundary conditions are given by:

$$\begin{array}{l}w\left(0,t\right)=0,\frac{\partial w\left(0,t\right)}{\partial x}=0,\\ EI\frac{{\partial }^{3}w\left(L_b,t\right)}{\partial x^3}=m_c\frac{{\partial }^{2}w\left(L_b,t\right)}{\partial t^2}+m_c\frac{D_c}{2}\frac{{\partial }^{3}w\left(L_b,t\right)}{\partial x\partial t^2},\\ EI\frac{{\partial }^{2}w\left(L_b,t\right)}{\partial x^2}=-m_c\frac{D_c}{2}\frac{{\partial }^{2}w\left(L_b,t\right)}{\partial t^2}-\left(J+m_c{\left(\frac{D_c}{2}\right)}^2\right)\frac{{\partial }^{3}w\left(L_b,t\right)}{\partial x\partial t^2}.\end{array}$$

(10)

where the equivalent bending stiffness and mass can be expressed as:

$$EI=EIb\frac{x_1}{L_b}+EIbp\frac{L_p}{L_b}+EIb\frac{L_b-x_1-L_p}{L_b},$$

$$meq=\frac{w_b{\rho }_b{L}_b{t}_b+w_p{\rho }_p{L}_p{t}_p}{L_b},$$

where

$$\begin{array}{l}EIb=\frac{E_b{w}_b{t}_b^3}{12},\\ EIbp=\frac{1}{3}{E}_b{w}_b\left(y{y}_1^3-y{y}_0^3\right)+\frac{1}{3}{E}_p{w}_p\left(y{y}_2^3-y{y}_1^3\right),\\ yy=\frac{E_p{w}_p{t}_p^2+E_b{w}_b{t}_b^2+2E_p{w}_p{t}_p{t}_b}{2\left(E_p{w}_p{t}_p+E_b{w}_b{t}_b\right)},\\ y{y}_0=-yy,\\ y{y}_1=t_b-yy,\\ y{y}_2=\left(t_b+t_p\right)-yy.\end{array}$$

We express the displacement in the form:

$$w\left(x,t\right)={\displaystyle \sum _{i=1}^n{\phi }_i\left(x\right)r_i\left(t\right)}$$

(11)

Bibo et al. [33] proved that for galloping motion, the first mode in the Galerkin procedure is enough to obtain the accurate response of the system, so the first mode results of the system are taken into account in the subsequent calculation. Therefore, the control equations of the first order vibration model of GPEH-V based on galloping are as follows:

$$\begin{array}{l}{\ddot{r}}_1+2{\zeta }_1{\omega }_1{\dot{r}}_1+{\omega }_1^2{r}_1-{\theta }_{1}V\left(t\right)=\\ \frac{1}{2}{\rho }_{a}D{L}_0{U}_0{a}_1\left({\phi }_1\left(L_b\right)+\frac{D}{2}{{\phi }^{\prime }}_1\left(L_b\right)\right)\left[{\phi }_1\left(L_b\right)+\frac{D}{2}{{\phi }^{\prime }}_1\left(L_b\right)\right]{\dot{r}}_1+\\ \frac{1}{2U_0}{\rho }_{a}D{L}_0{a}_3\left({\phi }_1\left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right){\left[{\phi }_1\left(L_b\right)+\frac{D}{2}{{\phi }^{\prime }}_1\left(L_b\right)\right]}^3{\dot{r}}_1^3,\end{array}$$

(12)

$$C_p\dot{V}+{\theta }_1{\dot{r}}_1+\frac{V}{R}=0$$

(13)

where

$$\begin{array}{l}{C}_p=\frac{{\epsilon }_{33}{w}_p{L}_p}{2t_p},\\ {\theta }_1=\frac{e_{31}{w}_p\left(t_p+t_b\right)}{2}\left({{\phi }^{\prime }}_1\left(x_1+L_P\right)-{{\phi }^{\prime }}_1\left(x_1\right)\right).\end{array}$$

The detailed derivation process of GPEH-V control equation using generalized Hamiltonian principle can refer to Abdelkefi et al. [39] and Varoto [40]. By solving Equation (12) with the Runge–Kutta method, the numerical solution of the displacement and voltage response of the piezoelectric coupling system can be obtained.

The expression of the total efficiency of the energy harvester or the electromechanical efficiency of the air is defined as:

$$\eta =\frac{P_{avg}}{P_f}$$

(14)

where P_avg is the average resonance power output of the energy harvester. P_f is the mechanical power output by the fluid, and its expression is proportional to the cube of the external velocity [31]:

$$P_f=\frac{D{L}_0}{2}{\rho }_a{U}_0^3$$

(15)

3. Galloping Characteristics Analysis

Galloping is a kind of self-excited vibration. The motion of the object causes the force along the direction of motion of the object. At this time, the fluid works on the object. When the force exceeds the energy dissipation of the system, dynamic instability occurs, so the vibration is induced. It develops into a large amplitude vibration. In other words, when the direction of action of the lift is consistent with the direction of motion of the structure, energy is continuously obtained from the fluid, and the amplitude is gradually increased. This process continues until the damping has a non-linear effect or other non-linear effects. The frequency of lift variation is the same as the natural frequency of the structure. This frequency is generally much lower than the frequency f_S of S_t, and its reduction speed V_r is more than 10.

3.1. The Influence Lateral Force Coefficient on the Critical Velocity

To calculate the critical velocity U_g of galloping wind velocity (i.e., the initial velocity of galloping), a linear analysis of the system is carried out. In Equation (12), the motion equation can be expressed in the state space by neglecting the nonlinear term of the fluid force. The correlation linear matrix of the state space formula is written as follows:

$$A=\left[\begin{array}{ccc}0& 1& 0\\ -{\omega }_1^2& -2{\zeta }_1{\omega }_1+\frac{{\rho }_{0}D{L}_0{U}_0}{2}{a}_1{\left[{\phi }_1\left(L_b\right)+\frac{D}{2}{{\phi }^{\prime }}_1\left(L_b\right)\right]}^2& {\theta }_1\\ 0& -\frac{{\theta }_1}{C_p}& -\frac{1}{R{C}_p}\end{array}\right]$$

(16)

According to Equation (12), when the linear term of the fluid force is equal to the damping force of the system, the critical velocity can be calculated as:

$$U_g=\frac{4{\zeta }_1{\omega }_1}{{\rho }_{0}D{L}_0{a}_1{\left[{\phi }_1\left(L_b\right)+\frac{D}{2}{{\phi }^{\prime }}_1\left(L_b\right)\right]}^2}$$

(17)

It can be seen from Equation (16) that linear matrix A has three eigenvalues, and the first two eigenvalues are complex conjugate term ${\lambda }_1={\overline{\lambda }}_2$; however for all wind velocity, the third eigenvalue is negative. The real part of the conjugate complex reflects the coupling damping of the system and also controls the dynamics of the system. As the wind velocity in the environment gradually increases, the linear term of the fluid force of the system gradually increases. When the linear term is greater than the electromechanical damping of the GPEH-V, the real part of the complex conjugate becomes positive, and the system begins to gallop. Any wind velocity greater than the critical value will make the real part of ${\lambda }_1$ and ${\lambda }_2$ positive, and make the system galloping. It is worth noting that the natural frequency ω₁ represents the first natural frequency of the system, because the first mode is enough to produce galloping oscillation. The experimental and theoretical results show that the value of ω₁ is between 13 Hz and 14 Hz, which is of great value for obtaining energy from piezoelectric materials.

It can be seen from Equation (17) that when the piezoelectric cantilever beam of the system is confirmed, the critical velocity is mainly related to the area of the windward side of the cylinder $D{L}_0$ and the coefficient of the linear term of the fluid force a₁. When the area of the windward side of the cylinder is constant, the value of a₁ determines the critical velocity. It can be seen from Equation (17) that the critical velocity U_g is inversely proportional to a₁, and when a₁ increases, U_g decreases. The value of a₁ is determined by the section shape of the cylinder. When the section of the cylinder is square, a₁ = 2.3 [34]. Therefore, a₁ is supposed to study near 2.3 in this paper. As showed in Figure 2, when the value of a₁ changes nearly 1–4, the corresponding critical velocity U_g changes between 1–7 m/s, and the critical velocity decreases with the diminish of mass of the cylinder.

3.2. Fluid Force Analysis

When the piezoelectric sheet is not attached, the system can be simplified into a single degree of freedom system:

$$\ddot{r}+2\zeta \omega \dot{r}+{\omega }^{2}r=F_y$$

(18)

where

$$F_y=\frac{1}{2}{\rho }_{0}D{L}_0{U}_0^2\left[a_1{\left(\phi \left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right)}^2\frac{\dot{r}}{U_0}+a_3{\left(\phi \left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right)}^4{\left[\frac{\dot{r}}{U_0}\right]}^3\right]$$

(19)

where F_y is the instantaneous fluid force in the transverse direction to the incident flow.

In this study, we used quasi-static assumptions to express fluid forces. This is because the essential feature of the transverse galloping phenomenon is that the timescale of the system vibration ($\sim 2\pi /\omega$) is much larger than the characteristic timescale of the fluid flow ($\sim D/U_0$). For the convenience of theoretical calculations, the fluid force F_y can be approximated as a polynomial function of the cylinder velocity. For our purposes a cubic polynomial can be used to approximate the vertical fluid force [34]. Where a₁ and a₃ are empirical coefficients of fluid force. When the length of the cylinder is constant, their values are only related to the cross-sectional shape of the cylinder [31]. Since the cylinder vibrates with respect to the direction of the flow velocity, there is only odd harmonics in the representation of the fluid force. From the analysis in the previous section, it can be known that only when a₁ > 0, is the system likely to gallop. Empirical coefficient a₃ represents the nonlinear dependence of the fluid force, which is negative. When a₃ is positive, the fluid force will increase unrestricted as the cylinder moves faster, which is impossible in practice. It is worth noting that a₁ > 0 and a₃ < 0 are significant for subsequent studies.

3.3. Galloping Response and Conversion Factor

It can be known from Equation (18) that the system can be theoretically solved by asymptotic methods if the nonlinear term is small [31]. In the case that both aerodynamic and damping forces are small compared with inertia and stiffness forces, solutions to Equation (18) will tend to a limit cycle of quasi-harmonic oscillations.

In this thesis, the Krylov–Bogoliubov [31] method is used to solve Equation (18), and the first-order approximate analytical solution of the system (see Appendix A) can be obtained as Equation (20).

$$r=y\sin \left(\omega t\right)$$

(20)

where

$$y=\sqrt{\left(\frac{8\left(-2\zeta \omega +\frac{1}{2}{\rho }_{0}D{L}_0{U}_0{\left(\phi \left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right)}^2{a}_1\right)}{-3{\rho }_{0}D{L}_0{\omega }^2{\left(\phi \left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right)}^4\frac{a_3}{U_0}}\right)}$$

(21)

From Equation (21), one can find the magnitude of the system amplitude depending on the flow velocity U₀ and the cross-sectional shape of the cylinder (a₁ and a₃), when the mechanical properties of the system and the surrounding fluid are constant.

The corresponding analytical amplitude curves of the cylinder with different flow velocity U values are plotted in Figure 3 which are basically in agreement with the numerical solution. It is shown that when the flow velocity exceeds the critical value, the limit cycle oscillations take place owing to the presence of nonlinearities. This local bifurcation is named Hopf bifurcation.

Figure 4, Figure 5 and Figure 6 show the amplitude of the cylinder under the cross-sectional shape of the cylinder (a₁ and a₃). They show that when the wind speed is constant, the amplitude of the cylinder increases with the increase of a₁ and a₃. From Figure 4, we can find that a₁ also has a critical value.

The mechanical efficiency of the air:

$${\eta }_{f-m}=\frac{P_m}{P_f}$$

(22)

where

$$P_m=\frac{1}{T}{\displaystyle {\int }_0^T{F}_y\phi \left(L_b\right)\dot{r}dt}$$

(23)

Finally, consider Equations (15), (20)–(23), where one obtains:

$${\eta }_{f-m}=\frac{a_1}{2}\left[\phi \left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right]{\left(\frac{y}{U_0}\right)}^2+\frac{3a_3}{8}{\left[\phi \left(L_b\right)+\frac{D}{2}{\phi }^{\prime }\left(L_b\right)\right]}^3{\left(\frac{y}{U_0}\right)}^4$$

(24)

4. Experiment and Identification of Empirical Coefficient of Fluid

This part is a wind tunnel experiment on the galloping motion of V-shaped square cylinder with a single degree of freedom in the “low turbulence reflux wind tunnel” finished at the Department of Mechanics, Tianjin University, and the experimental results are compared with the theoretical and static numerical results. The experimental results are used to verify the correctness of the theory and simulation. The physical figure of the main body of the wind tunnel tested is shown in Figure 7 and Figure 8. The MFC piezoelectric sheet is connected in series with the rheostat. The voltage signals at both ends of the rheostat are output to the data collector. The laser displacement sensor is used to monitor the displacement signal of the cylinder and input them to the data collector. The computer controls the operation and storage of the data acquisition system.

Piezoelectric material and its structure are the key points of the piezoelectric energy harvester. The selection of piezoelectric material should consider the working environment, material performance, and other factors. In the existing literature of wind-induced vibration wind energy harvester, the majority of lead zirconate titanate piezoelectric ceramics (PZT) materials are used, which is due to that piezoelectric ceramics have large piezoelectric constant, dielectric constant, and electromechanical coupling coefficient, and since their performance is excellent, a small deformation can obtain a large voltage. Due to the irregular vibration and high frequency of the structure in the wind, and friableness of the piezoelectric ceramic material, it is easy to fatigue fracture under long-term vibration, and since the stiffness is large, it is difficult to deform under the action of general wind. Therefore, another common piezoelectric material, piezoelectric composite material, is selected in this paper, which has small mass, large flexibility, and does not easily incur fracture fatigue. According to the characteristics of wind-induced vibration, the structure that can harvest wind energy is designed and an experimental platform is built, as shown in Figure 9.

In order to study the influence of cross-section shape on fluid force, seven groups of comparative experiments were carried out, with the cross-section shape of the cylindrical structure at the end of the cantilever beam changed, as shown in

Table 2. The cross-section shape of a cylinder.

Model Number	1	2	3	4	5	6	7
Section shape
$\alpha$	90	13	20	27	34	45	64

. Model 1 is a cylinder with a square cross-section, which mainly verifies the theoretical and numerical results, and can also be compared with another six models. The cross-section shape of models 2–7 is controlled by the angle $\alpha$. The seven groups of models are made by 3D printing technology.

In order to accurately analyze the force acting on the cylinder and calculate the flow force coefficient, it is necessary to simplify the experimental model and eliminate the unnecessary interference factors. Therefore, the wind tunnel experiment without sticking to the MFC piezoelectric sheet is carried out first.

When the wind speed is 5 m/s, the relationship between the time and the displacement of the cylinder of model 1 is shown in Figure 10. It can be seen from Figure 10 that there are some unavoidable errors in the acquisition of experimental data. In order to reduce the errors in the subsequent calculation, the singular value decomposition (SVD) method is used to reduce the noise of the data. Secondly, the SINDy method [37] is used to identify the fluid force coefficient a_i in the formula. Assume that a₁ and a₃ are sufficient to represent the fluid forces. The identification result is a₁ = 2.2961 and a₃ = −18.7142, which is in line with Abdelkefi [34] and Paulo [40]. The coefficients of fluid force used by others is almost the same. Figure 11 shows the relationship between the fluid velocity and the amplitude of the cylinder, where the solid line is the result of theoretical calculation using the identified fluid force coefficient, and the circle represents the experimental results. By comparison, it is found that the theory is almost consistent with the experiment.

The experimental data do not confirm the theoretical or numerical solution when the flow velocity is between about 1.5 and 3 m/s in Figure 11. This is because when the velocity is near the critical velocity, Hopf bifurcation occurs. The amplitude may suddenly increase or decrease with a slight disturbance. The trend of experimental data is basically the same as that of theoretical and numerical results, when the velocity is far away from the critical velocity and it is more according to the numerical solution.

Depending on the above analysis method, the fluid coefficients of the remaining six models are identified in turn, and seven groups of energy harvesters are calculated numerically.

Table 3. Summary of seven groups of energy harvester characteristics at wind velocity of 10 m/s.

Section Shape
Model number	1	2	3	4	5	6	7
$\alpha$	90	13	20	27	34	45	65
$a_1$	2.2961	1.2538	1.7823	1.6411	1.7067	1.9334	2.3276
$a_3$	−18.7142	−6.2630	−14.1578	−9.1829	−9.6315	−8.2705	−25.0288
$U_g$	2.56	4.69	3.30	3.58	3.45	3.04	2.53
Power of numerical solution (W)	5.68×10–4	4.66×10–4	4.67×10–4	6.12×10–4	6.70×10–4	9.26×10–4	4.16×10–4
Power of experimental solution (W)	5.94×10–4	4.84×10–4	5.57×10–4	5.84×10–4	6.90×10–4	9.30×10–4	3.82×10–4
Total efficiency of numerical solution (%)	0.044	0.036	0.036	0.048	0.051	0.071	0.033
Total efficiency of experimental solution (%)	0.044	0.036	0.036	0.047	0.052	0.072	0.032

shows the detail characteristics of seven groups of energy harvesters. It is clear that the power of the numerical solutions are basically consistent with the experimental solutions; so are the total efficiency solutions. These demonstrate that the empirical coefficients of fluid obtained by the sparse recognition method is basically in line with the actual situation.

From Figure 12, we can conclude that model 7 has the minimum critical velocity, while its output power is the least at a flow velocity of 10 m/s from Figure 13. Whereas, model 6 does not have the minimum critical velocity. Nevertheless, its output power is the largest. It is extremely important to consider the characteristics of different models according to the actual situation.

5. Energy Harvesting Characteristics of Seven Experimental Models

To design an efficient and optimized GPEH-V, an accurate prediction of the levels of output voltage, output power, and total efficiency of the energy harvester is of the utmost importance. Through the analysis of Section 4, we can accurately obtain the fluid coefficients of the system. The output voltage of the energy harvester under different flow velocity can be calculated by substituting the fluid coefficients into Equations (12) and (13) of Section 2. In the previous section, we know that the numerical solution can analyze the energy harvester characteristics of each model more accurately. Therefore, we use the Runge–Kutta method and MATLAB software to solve Equations (12) and (13) numerically to get the voltage of the system under different flow velocity. According to Equations (14) and (15), the output power and the total energy conversion efficiency of the energy harvest can be calculated.

The flow velocity and output voltage of the energy harvesters are shown in Figure 14. The asterisk represents the experimental results and the solid line in the figure represents the numerical results. It can be seen from the figure that when the wind velocity of model 2 is 6–7 m/s, the experimental results are quite different from the numerical results. This is because the critical velocity of model 2 is 4.7 m/s. According to the previous analysis, Hopf bifurcation occurs near the critical velocity. The experimental results may suddenly increase or decrease with a slight disturbance. In the low wind speed environment, the output voltage of model 6 is the largest, which is because the appropriate sharp angle is conducive to the falling off of the fluid vortex, making the cylinder produce greater pressure difference in the direction of motion. The lift of the cylinder is increased, which leads to an increase in the oscillation amplitude and the output voltage of the energy harvest.

The flow velocity and output power of the energy harvesters are shown in Figure 15. It shows that when the wind speed is about 4 m/s, model 1 and model 7 are the first to output electric energy. With the increase of wind speed, the output energy of model 6 is gradually higher than the other six models. Model 2 starts to output electric energy when the wind speed is greater than 6 m/s, and its output power is the lowest among the seven models when the wind speed is between 6 m/s and 9 m/s. After the wind speed is greater than 9 m/s, the output power gradually exceeds model 7.

The flow velocity and the total efficiency of the energy harvesters are shown in Figure 16. When the wind speed is about between 4 m/s and 5 m/s, the total conversion efficiency of model 1 is the highest. When the wind speed is about between 5 m/s and 10 m/s, the conversion efficiency of model 6 is far higher than the other six models, while model 2 is the lowest.

Overall, model 6 can generate higher voltage and output higher electric energy in the low wind speed environment of 5–10 m/s, and more effectively convert wind energy into electric energy. Figure 17 and Figure 18 show the voltage and output power of model 6 under different resistance loads in different flow velocities. The solid line in the figure represents the numerical results and the asterisk represents the experimental results. It can be obviously seen from the Figure 17 and Figure 18 that when the resistance value changes, the output voltage and power under different resistance values can be accurately predicted by the identified fluid force coefficient. The two figures also show when the resistance of the load increases, the output voltage and power increases.

6. Conclusions

In this paper, a novel architecture of GPEH-V is proposed to harvest the vibration energy caused by air flow. The nonlinear distributed parameter model is established for the large deflection transverse vibration based on the Euler–Lagrange method. The one-order vibration equation is obtained via the Galerkin method. Based on the theoretical model, cross-section geometry, critical velocity, energy efficiency (or conversion factor), and mechanical efficiency of the air was established. The galloping response was solved by the Krylov–Bogoliubov method, and the first-order approximate analytical solution of the system was obtained. Hopf bifurcation was analyzed theoretically and observed in the experiment. Seven groups of comparative experiments were carried out in order to study the influence of V-shaped cross-section on fluid force. The SVD and SINDy methods were used to identify the fluid force coefficient. Finally, combined with the experimental results and theoretical analysis, the critical velocity and energy efficiency of GPEH-V with different V-shaped cross-sections were studied. The conclusions are as follows:

• The V-shaped groove on the windward of the square cylinder can affect not only the critical velocity of galloping but also the lateral force acting on it. The optimal geometry of the V-groove makes the GPEH-V produce the maximum output power in the low wind speed environment of 5–10 m/s.
• Based on the coefficient, the bifurcation mode of the GPEH-V system can be predicted accurately by the Krylov–Bogoliubov method.
• It may accurately identify the flow force coefficients using the SVD and SINDy methods. It also shows that the identified coefficients can accurately predict the energy acquisition characteristics under different wind speeds and different resistance loads.
• The deformation of the piezoelectric cantilever can be accurately predicted by the numerical calculation method to prevent the irreversible damage of MFC due to excessive deformation.

In the follow-up research, we will use the finite element method to analyze fluid–solid coupling vibration of the optimized GPEH-V, reveal the influence of the shape of the windward face on the vibration characteristics theoretically, and optimize the circuit part of the GPEH-V to comprehensively improve the energy harvester efficiency.