A Piezoelectric Energy Harvester with Bending–torsion Vibration in Low-speed Water

This paper presents a piezoelectric energy harvester using an eccentric cylinder undergoing bending–torsion vibration in low-speed water. It can harvest energy from water using vortex-induced vibration (VIV). A distributed parameter beam model with respect to the motion of the piezoelectric beam was established based on Euler–Bernoulli beam theory. The governing coupled equations of the harvester system were derived by Lagrange's equations. The optimal configurations and work conditions of harvesters were numerically analyzed according to the above mathematical models. Experiments were designed and performed to verify the numerical results. The numerical results were in good agreement with the experiment results, which verifies the validity of the mathematical models. The harvester with bending–torsion vibration generated an output power of 0.3978 mW, which is 1.99 times of that of the harvester with a solid-cylinder tip undergoing bending only.


Introduction
The demand for intermittent low-power, wireless, and small electronic devices is growing.More and more scholars have focused on energy harvesting from vibration since the electric energy can be converted from vibration using piezoelectric materials [1].Mitcheson et al. [2] suggested that a piezoelectric harvester with a non-resonant structure was the best choice for harvesting energy from low frequency ambient vibration.There are many investigations of piezoelectric harvester with respect to the modeling [3,4], structures [5,6] and piezoelectric materials [7,8].However, studies are mainly focused on the numerical analysis and experiments using harmonic base excitation.
Some scholars began to study a piezoelectric harvester by utilizing ambient vibrations.For example, Robbins et al. [9] studied a piezoelectric energy harvester which could harvest energy via the wind.Akaydin et al. [10] investigated a self-excited fluidic energy harvester, which was tested in a wind tunnel and 0.1 mW of electrical power was produced.Abdelkefi et al. [11] researched the phenomena of piezoelectric energy harvesting from freely oscillating cylinders undergoing vortex-induced vibration (VIV).Linear and nonlinear models were analyzed and it was found that the load resistance influences the onset of the synchronization region and the harvesters' characteristics.Energy harvesting by VIV has attracted considerable interest from researchers.The main research aspects are prototype design and fabrication [12,13], modeling, numerical and experimental analysis of the harvesters [14][15][16][17].As for the modeling of VIV-based energy harvesters, research on lift fluctuation and vortex interactions is necessary.Dai et al. [18][19][20] modeled the lift fluctuation with the existing mathematical representations [21][22][23].The model was applied to analyze the VIV of a cylinder linearly and nonlinearly.This was proven to be consistent with the practical situation by results of simulation and experiment.Goushcha et al. [24] investigated the fundamental mechanisms of vortex interactions with flexible structures, which were critical to fluidic electrical energy harvesting.A review of recent studies in the field of energy harvesting from aeroelastic vibrations was reported by Abdelkefi [25].Future recommendations in the field were discussed, including mathematical modeling, realistic loadings and small-sized power conditioning circuit optimization, and prototype fabrication of energy harvesters.Hence, scholars are going to continue to work on the above aspects.
The low-velocity water flow that can provide a low frequency vibration has been widely reported on.Taylor et al. and Allen et al. [26,27] researched an energy harvesting eel that utilized piezoelectric polymers to convert the flow energy in the oceans and rivers into electrical energy.Song et al. [28,29] studied an upright cylinder and a bicylinder vortex-induced piezoelectric energy harvester which harvest energy from water; the maximum output power of these was 84.49 µW and 21.86 µW, respectively.Shan et al. [30] researched the phenomena of piezoelectric energy harvesting using experiments with two tandem cylinders undergoing VIV in water, and the output power of 533 µW was obtained by the downstream harvester in specific conditions.However, the energy-harvesting ability of their single harvester was slightly weak.Hence, harvesters with high-efficiency in water have become an area of increasing research.
Therefore, this paper presents several piezoelectric energy harvesters with vortex-induced bending-torsion vibration.The piezoelectric beam undergoes bending-torsion with a tip eccentric cylinder, which creates a mass eccentric distance between its center of mass and the force target, and results in a coupling between the bending and torsion vibrations.Hence, there is a big effect on the harvesters' natural frequencies and mode shapes due to the bending-torsion coupling.In Section 2, governing coupled equations of the motion of the harvesting system are established.In Section 3, experiments are designed and performed.The work conditions and cylinder configurations of the harvester are numerically analyzed and experimentally investigated.Results discussions and conclusions are given in Section 4.

Modeling of the Harvester System
Figure 1a shows the energy harvester using piezoelectric cantilever beam (PZT) and its working environment.The energy harvester is composed of a piezoelectric cantilever beam and an eccentric cylinder tip. Figure 1b shows the schematic of energy harvester.Table 1 lists the geometric and physical properties of the energy harvesting system.The symbolic variables e and d are the eccentricity distance and diameter of cylinder eccentric hole, respectively, which determine the value of the eccentric distance of mass l d .
the existing mathematical representations [21][22][23].The model was applied to analyze the VIV of a cylinder linearly and nonlinearly.This was proven to be consistent with the practical situation by results of simulation and experiment.Goushcha et al. [24] investigated the fundamental mechanisms of vortex interactions with flexible structures, which were critical to fluidic electrical energy harvesting.A review of recent studies in the field of energy harvesting from aeroelastic vibrations was reported by Abdelkefi [25].Future recommendations in the field were discussed, including mathematical modeling, realistic loadings and small-sized power conditioning circuit optimization, and prototype fabrication of energy harvesters.Hence, scholars are going to continue to work on the above aspects.
The low-velocity water flow that can provide a low frequency vibration has been widely reported on.Taylor et al. and Allen et al. [26,27] researched an energy harvesting eel that utilized piezoelectric polymers to convert the flow energy in the oceans and rivers into electrical energy.Song et al. [28,29] studied an upright cylinder and a bicylinder vortex-induced piezoelectric energy harvester which harvest energy from water; the maximum output power of these was 84.49 μW and 21.86 μW, respectively.Shan et al. [30] researched the phenomena of piezoelectric energy harvesting using experiments with two tandem cylinders undergoing VIV in water, and the output power of 533 μW was obtained by the downstream harvester in specific conditions.However, the energyharvesting ability of their single harvester was slightly weak.Hence, harvesters with high-efficiency in water have become an area of increasing research.
Therefore, this paper presents several piezoelectric energy harvesters with vortex-induced bending-torsion vibration.The piezoelectric beam undergoes bending-torsion with a tip eccentric cylinder, which creates a mass eccentric distance between its center of mass and the force target, and results in a coupling between the bending and torsion vibrations.Hence, there is a big effect on the harvesters' natural frequencies and mode shapes due to the bending-torsion coupling.In Section 2, governing coupled equations of the motion of the harvesting system are established.In Section 3, experiments are designed and performed.The work conditions and cylinder configurations of the harvester are numerically analyzed and experimentally investigated.Results discussions and conclusions are given in Section 4.

Modeling of the Harvester System
Figure 1a shows the energy harvester using piezoelectric cantilever beam (PZT) and its working environment.The energy harvester is composed of a piezoelectric cantilever beam and an eccentric cylinder tip. Figure 1b shows the schematic of energy harvester.Table 1 lists the geometric and physical properties of the energy harvesting system.The symbolic variables e and d are the eccentricity distance and diameter of cylinder eccentric hole, respectively, which determine the value of the eccentric distance of mass ld.To model the harvester, the following assumptions are adopted: (1) The piezoelectric beam is assumed to be inextensible.
(2) The attachments of the tip cylinder, beam and fixed end are assumed to be rigid.
Figure 2 shows a schematic diagram of the cantilever beam.L and b are the length and width of the beam (Baoding Hengsheng Acoustics Electron Apparatus Co., Ltd ® , Baoding, China), respectively.The thickness of the PZT layer made of PZT 5H defined in Table 1 is h p , and that of the substrate layer is h s .Two coordinate systems are used to describe the coupled bending-torsion motions of the beam.The reference coordinate system (O-XYZ) is fixed to the clamped side of the beam, which is adopted to describe the bending motion.s is used as the arc-coordinate in the length direction of the beam.In the fixed end, s = 0, and s = L in the tip.s + u(s, t), v(s, t) and w(s, t) are the displacements of the infinitesimal element at three directions in the coordinate system O-XYZ.Meanwhile, a local coordinate system (I-ξηζ) is adopted to describe the torsion motion.The axis I-ξ, I-η and I-ζ are the central principal axes of the cross section in three directions.
Appl.Sci.2017, 7, 116 3 of 16 To model the harvester, the following assumptions are adopted: (1) The piezoelectric beam is assumed to be inextensible.
(2) The attachments of the tip cylinder, beam and fixed end are assumed to be rigid.
Figure 2 shows a schematic diagram of the cantilever beam.L and b are the length and width of the beam (Baoding Hengsheng Acoustics Electron Apparatus Co., Ltd ® , Baoding, China), respectively.The thickness of the PZT layer made of PZT 5H defined in Table 1 is hp, and that of the substrate layer is hs.Two coordinate systems are used to describe the coupled bending-torsion motions of the beam.The reference coordinate system (O-XYZ) is fixed to the clamped side of the beam, which is adopted to describe the bending motion.s is used as the arc-coordinate in the length direction of the beam.In the fixed end, s = 0, and s = L in the tip.s + u(s, t), v(s, t) and w(s, t) are the displacements of the infinitesimal element at three directions in the coordinate system O-XYZ.Meanwhile, a local coordinate system (I-ξηζ) is adopted to describe the torsion motion.The axis I-ξ, I-η and I-ζ are the central principal axes of the cross section in three directions.
where the over dot denotes the derivative with respect to time t, and the vector arrow denotes the vector along the axes.
where the over dot denotes the derivative with respect to time t, and the vector arrow denotes the vector along the axes.ωξ, ωη, ωζ are the angular velocity of the rotating beam with respect to the axes I-ξ, I-η and I-ζ.According to the practical work environment of the coupled bending-torsion motions of the beam, it is necessary to keep the mathematical model as simple as possible.Compared with other displacements, the longitudinal vibration displacement u(s, t), and rotation angular displacements ψ(s, t) and θ(s, t) are so small that they can be neglected.As a result, u(s, t) = 0, ψ(s, t) = θ(s, t) = 0, and from Equation (1), ωξ = φ  , ωη = ωζ = 0.
Next, the Lagrange's equations are used to derive the governing equations of the motion of the harvesting system.
The kinetic energy of the system is the sum of the kinetic energies of the beam (Tb), the eccentric cylinder (Tc) and the additional kinetic energy of fluid (Tf).
( ) where ρb and Vb are the density and volume of the beam, which concludes two parts: the PZT and substrate layer.Iξ is the rotary inertia per unit length with respect to the axis I-ξ, given by: ( ) where h is the thickness of the beam.
where D and d are the diameter of the cylinder and the eccentric hole, respectively.Mc and Jcξ are the mass and the rotary inertia with respect to the axis I-ξ of the eccentric cylinder, as follows: According to the practical work environment of the coupled bending-torsion motions of the beam, it is necessary to keep the mathematical model as simple as possible.Compared with other displacements, the longitudinal vibration displacement u(s, t), and rotation angular displacements ψ(s, t) and θ(s, t) are so small that they can be neglected.As a result, u(s, t) = 0, ψ(s, t) = θ(s, t) = 0, and from Equation (1), Next, the Lagrange's equations are used to derive the governing equations of the motion of the harvesting system.
∂ ∂t The kinetic energy of the system is the sum of the kinetic energies of the beam (T b ), the eccentric cylinder (T c ) and the additional kinetic energy of fluid (T f ).
where ρ b and V b are the density and volume of the beam, which concludes two parts: the PZT and substrate layer.I ξ is the rotary inertia per unit length with respect to the axis I-ξ, given by: where h is the thickness of the beam.
where D and d are the diameter of the cylinder and the eccentric hole, respectively.M c and J cξ are the mass and the rotary inertia with respect to the axis I-ξ of the eccentric cylinder, as follows: where M f and J fξ are the fluid-added mass and rotary inertia with respect to the axis I-ξ of the fluid due to the vibration of the eccentric cylinder, given by: where C M = 1 is the fluid-added mass coefficient.
The potential energies of the whole system include the potential energy (U i ), the electric potential energy (U V ) of the beam and the gravitational potential energy (U g ) of the system.
where the strains and stresses in the substrate and PZT layers are: Due to the poling of the PZT layer, the strain and electric displacement vector are governed by the constitutive equations, given by: where E s and E p are the Young's modulus at constant electric field, G s and G p are the shear modulus, E 2 is the electric field intensity in the poling direction, and V(t) is the voltage between the two piezoelectric electrodes.
where V p is the volume of the PZT layer.
where x(s, t) is the displacement of the beam in three axes direction of the reference coordinate system.The virtual work includes the virtual work of the external resistance δW R , lift force δW F , buoyancy of fluid δW b , mechanical damping δW c and fluid-added damping δW cf .
where Q(t) is the charge between the piezoelectric electrodes.
where l d is the eccentric distance of mass and given by: and F(t) is the lift force.It is generated when the eccentric cylinder is subjected to the incoming fluid [22], as: where U 0 is the speed of the incoming fluid.C L is the vortex lift coefficient.The fluid variable q is interpreted as a reduced vortex lift coefficient, which is governed by Van der Pol equation [31], and ω f is the vortex shedding frequency. .. where where C L0 = 0.3 [32].And f is the nondimensional force of vortex [22], given by: The virtual work of the buoyancy of fluid δW b is: The virtual work of fluid-added damping δW cf is: where where γ = 0.8 [22] is the fluid-added damping coefficient.
where c m is the mechanical damping coefficient measured by experiment.In order to solve the above equations, the bending vibration displacements v(s, t), w(s, t) and torsion angular displacement φ(s, t) are expressed by using Galerkin procedure [4] in the form of where φ(s) and r(t) are the mode shapes and the model coordinates, respectively.They can be expressed as: where the A n , B n , C n and D n are arbitrary constants and λ ni are determined by the natural frequency of the piezoelectric beam system, given by: where EI η and EI ζ are the bending stiffness of the beam with respect to axes I-η, I-ζ.GI p is the torsion stiffness of the beam with respect to axis I-ξ, as follows: where h a , h b and h c are the positions of the layers given with respect to the neutral axis which are defined in Figure 4.
as follows: where φ(s) and r(t) are the mode shapes and the model coordinates, respectively.They can be expressed as: ( ) ( ) ( ) where the An, Bn, Cn and Dn are arbitrary constants and λni are determined by the natural frequency of the piezoelectric beam system, given by: where EIη and EIζ are the bending stiffness of the beam with respect to axes I-η, I-ζ.GIp is the torsion stiffness of the beam with respect to axis I-ξ, as follows:

b E h h E h h EI b E h h E h h EI
where ha, hb and hc are the positions of the layers given with respect to the neutral axis which are defined in Figure 4.

p s p p s s h h E h h h E h E h
as follows: ( ) According to the work condition, their associated boundary conditions are used to solve the mode shapes above.
Appl.Sci.2017, 7, 116 where δ ij is the Kronecker delta, defined as 1 when i is equal to j and 0 otherwise.As a consequence, combining the boundary conditions ( 30)-( 35) and orthogonality conditions (36) and (37), as well as Equation ( 23), mode shapes and natural frequency of the beam can be obtained.As the energy harvesting system is designed to work in the low frequency environment, the first mode shape [33] is accurate enough to research the energy harvesting in this paper.

ds
The virtual work equation: where the coefficient H n are expressed as: Now, the Lagrange function which is the subtraction of kinetic energy Equation (38) and potential energy Equation (39) is established and expressed as: This occurs when substituting Equation (40) and Lagrange function (41) into Lagrange's Equation (2), and assigning variables r and V, respectively.Then, the following state equations are obtained: ..
where the coefficients of the state variables are expressed as: Substituting the Equations ( 19)-( 21) into Equation ( 15) the nondimensional force of vortex f is rewritten as: where the coefficient H 7 is given by: According to Equations ( 42)-( 44), a space state equation of the energy harvesting system is established, and expressed as: where X is the space state variable, written as: Equation (45) can be numerically solved by using the solver of ode45 in MATLAB (MathWorks, Inc ® , Natick, MA, USA).

Experimental Setup
Harvesters are fabricated with piezoelectric cantilever beams and tip eccentric cylinders.The beam has an aluminum layer and a PZT layer, which are bonded together by epoxy.The tip cylinder is made of acrylic with an eccentric hole.They are joined by a symmetry rectangular clamping fixture made of polyamide.
The experiment was carried out in a test rig system, as shown in Figure 5.This test system is composed of a water channel, pump (Shimge Co., Ltd ® ; Hangzhou, China), frequency converter (Shenzhen Junhui Electronics Co., Ltd ® , Shenzhen, China), data acquisition and processing system.The water channel is mainly composed of three sections: the setting chamber, (equipped with) a cellular device (Qingdao Tonglide plastic buzzer hives Co., Ltd ® , Qingdao, China), and damping meshes (Hangzhou Kuangshi Co., Ltd ® , Hangzhou, China), aiming to steady the flow.The contraction section returns the flow speed of water from the setting chamber.The experiments were performed in the working section.The pump makes the water cycle in the channel.The frequency converter is used to vary the work frequency of the pump, regulating the discharge of water, and then varying the flow speed of water.The data acquisition and processing systems include a DAQ of NI (National Instruments, Austin, TX, USA), and a PC (Lenovo Group Ltd ® , Beijing, China), which can measure and process the output voltage across the external resistance R in real time.
The water channel is mainly composed of three sections: the setting chamber, (equipped with) a cellular device (Qingdao Tonglide plastic buzzer hives Co., Ltd ® , Qingdao, China), and damping meshes (Hangzhou Kuangshi Co., Ltd ® , Hangzhou, China), aiming to steady the flow.The contraction section returns the flow speed of water from the setting chamber.The experiments were performed in the working section.The pump makes the water cycle in the channel.The frequency converter is used to vary the work frequency of the pump, regulating the discharge of water, and then varying the flow speed of water.The data acquisition and processing systems include a DAQ of NI (National Instruments, Austin, TX, USA), and a PC (Lenovo Group Ltd ® , Beijing, China), which can measure and process the output voltage across the external resistance R in real time.It is found that the power increases at first then decreases with the increase of the resistance value.There is an optimal resistance around 170 and 140 kΩ for the output power when the widths of the beams are 20 and 25 mm, respectively.From Figure 6, one can find that the optimal resistance is not related to the flow velocities, and the optimal resistance of the numerical analysis is consistent with that of the experiment.It is found that the power increases at first then decreases with the increase of the resistance value.There is an optimal resistance around 170 and 140 kΩ for the output power when the widths of the beams are 20 and 25 mm, respectively.From Figure 6, one can find that the optimal resistance is not related to the flow velocities, and the optimal resistance of the numerical analysis is consistent with that of the experiment.Table 2 shows the configurations of eccentric cylinder used in experiment.To document the measured results, it is necessary to number the cylinders.The sequence numbers represent the diameter of the cylinder (D) and the diameter of the eccentric hole (d).For example, the number Con 2006 represents a cylinder with D = 20 and d = 6 mm.In addition, the number Con 2000 represents a solid cylinder with a diameter of 20 mm.     Figure 7 illustrates the average power of harvesters with the configurations in Table 2, which are numerically calculated and experimentally measured at various flow velocity U 0 with the same resistance of 80 kΩ.It can be found that the output power P first increases to P max (maximum value of power) with U 0 , while U 0 continues to increase, P decreases.The optimal velocity is about 0.339-0.366m/s for the nine configurations.As the flow velocity U 0 rises up, the vortex shedding frequency ω f increases, which is gradually close to the natural frequency of the harvester.As a result, the lift force from the vortex drives the eccentric cylinder to resonance, and leads to a significant increase in amplitude of the cylinder rigidly connected with piezoelectric beam.At the same time, the optimal velocity becomes greater with the increase of d, which causes a lighter mass of the cylinder.As a consequence, the natural frequency of the harvester increases.Meanwhile, there should be a positive correlation between P max and l d dominated by the eccentricity distance e and the diameter of cylinder eccentric hole d.It can be concluded that P max increases with the increase of eccentricity of mass l d , as shown in Figure 7a and Table 3. P max is 0.3978 mW with a maximum value of l d in configurations 2006, 2008, and 2010 in Figure 7a.Besides, the P max of the solid cylinder is 0.1999 mW, and there is an improvement of 99%.The torque excited by vortex becomes bigger when l d increases.It leads to a bigger twist angle.Further, more energy is harvested.

Numerical and Experimental Analysis
In Figure 7, it also can be summarized that the cylinder diameter D has a negative influence on P max .P max decreases from 0.3978 to 0.2482 mW, while D increases from 20 to 30 mm with a similar l d in Con 2010, Con 2510, and Con 3010.In addition, the same conclusion can be obtained from Table 3.The lift force reduces with the increase of diameter D, which weakens cylinder vibration response.By comparing the results of different experiments, the best configuration and working environment for the harvester is summarized.The optimal load resistance is 170 kΩ when b = 20 mm, and 140 kΩ when b = 25 mm.The greater the l d is, the greater P max will be.Therefore, l d should be designed as large as possible if the structure processing permits.Similarly, a harvester with a smaller cylinder diameter D has a greater P max , and D is the minimum value of 20 mm in this paper.
The best configurations and working environments for the harvester are concluded by experiments.The validity of the model derived is verified.As a result, the model can be used to numerically analysis the harvesting characteristics of harvesters in varieties of situations and configurations.Figure 8 shows the numerical analysis results of the harvester with D = 20 mm and b = 20 mm at the optimal load resistance.Compared with the harvester undergoing bending only with an output power of 0.2188 mW shown in the bottom of Figure 8, the output power of the harvester with configurations d = 10 and e = 5 mm undergoing bending-torsion is 0.4854 mW.This result implies that the present work is meaningful for the design of harvesters.

Figure 1 .
Figure 1.Energy harvester: (a) Harvester and its working environment; (b) Schematic of harvester.Figure 1.Energy harvester: (a) Harvester and its working environment; (b) Schematic of harvester.

Figure 1 .
Figure 1.Energy harvester: (a) Harvester and its working environment; (b) Schematic of harvester.Figure 1.Energy harvester: (a) Harvester and its working environment; (b) Schematic of harvester.

Figure 6
Figure 6 illustrates the numerical and experimental values of average power versus external load resistance R. The configurations of the harvester are D = 20, d = 10 and e = 5 mm.The powers are measured at flow velocities of 0.355 and 0.404 m/s, respectively.It is found that the power increases at first then decreases with the increase of the resistance value.There is an optimal resistance around 170 and 140 kΩ for the output power when the widths of the beams are 20 and 25 mm, respectively.From Figure6, one can find that the optimal resistance is not related to the flow velocities, and the optimal resistance of the numerical analysis is consistent with that of the experiment.

Figure 6
Figure 6 illustrates the numerical and experimental values of average power versus external load resistance R. The configurations of the harvester are D = 20, d = 10 and e = 5 mm.The powers are measured at flow velocities of 0.355 and 0.404 m/s, respectively.It is found that the power increases at first then decreases with the increase of the resistance value.There is an optimal resistance around 170 and 140 kΩ for the output power when the widths of the beams are 20 and 25 mm, respectively.From Figure6, one can find that the optimal resistance is not related to the flow velocities, and the optimal resistance of the numerical analysis is consistent with that of the experiment.

Figure 7 .
Figure 7. Values of average power numerically calculated and experimentally measured at various flow velocities with different cylinder configurations of energy harvesters with a beam width of 25 mm (a) D = 20 mm; (b) D = 25 mm; (c) D = 30 mm.

Figure 7 .
Figure 7. Values of average power numerically calculated and experimentally measured at various flow velocities with different cylinder configurations of energy harvesters with a beam width of 25 mm (a) D = 20 mm; (b) D = 25 mm; (c) D = 30 mm.

Table 1 .
The geometric and physical properties of energy harvesting systems.

Table 1 .
The geometric and physical properties of energy harvesting systems.ρpDensity of the piezoelectric layer (kg/m 3 ) 7800 ρs Density of the substrate layer (aluminum; kg/m 3 ) 2700 ρf Fluid density (kg/m 3 )

Table 2 .
Configurations of the eccentric cylinder.

Table 2 .
Configurations of the eccentric cylinder.

Table 3 .
P max of different configurations of the eccentric cylinder when b = 20 mm.