A Piezoelectric Energy Harvester with Bending–Torsion Vibration in Low-Speed Water

Abstract

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.

1. 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.

2. 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. The geometric and physical properties of energy harvesting systems.

ρp	Density of the piezoelectric layer (kg/m3)	7800
ρs	Density of the substrate layer (aluminum; kg/m3)	2700
ρf	Fluid density (kg/m3)	1000
ρc	Density of the cylinder (acrylic; kg/m3)	1200
Ep	Young modulus of the piezoelectric layer (Gpa)	59.77
Es	Young modulus of the substrate layer (Gpa)	70
L	Length of the beam (mm)	85
b	Width of the beam (mm)	25
hp	Thickness of the piezoelectric layer (mm)	0.2
hs	Thickness of the substrate layer (mm)	0.2
e31	Piezoelectric constant (C/m2)	−16.6
ε11	Piezoelectric permittivity (nF/m)	41.78
Lc	Diameter of the cylinder (mm)	75
Lf	Length of the cylinder submerged in fluid (mm)	60

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.

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.

(3)

The beam is an Euler–Bernoulli beam.

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.

The local coordinate system I-ξηζ can be obtained by coordinating the rotation of the reference coordinate system, as shown in Figure 3. Firstly, O-X₁, O-Y₁ are obtained by a counter-clockwise rotation of O-X, O-Y with respect to the axis O-Z by the angle ψ. Secondly, I-ξ, O-Z₁ are obtained by a rotation of O-X₁, O-Z with respect to the axis O-Y₁ by the angle θ. Finally, I-η, I-ζ are obtained by a rotation of O-Y₁, O-Z₁ with respect to the axis O-ξ by the angle φ. As a result, the angular velocity of the three rotation processes can be derived by the coordinate transformation, as shown in Equation (1).

$$\begin{array}{l}\omega \left(s,t\right)=\dot{\psi }\stackrel{\rightharpoonup }{Z}+\dot{\theta }\overrightarrow{Y_1}+\dot{\varphi }\overrightarrow{\xi }\\ =(\dot{\varphi }-\dot{\psi }\sin \theta )\overrightarrow{\xi }+(\dot{\psi }\cos \theta \sin \varphi +\dot{\theta }\cos \varphi )\overrightarrow{\eta }+(\dot{\psi }\cos \theta \cos \varphi -\dot{\theta }\sin \varphi )\overrightarrow{\zeta }\\ ={\omega }_{\xi }\overrightarrow{\xi }+{\omega }_{\eta }\overrightarrow{\eta }+{\omega }_{\zeta }\overrightarrow{\zeta }\end{array}$$

(1)

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), ω_ξ = $\dot{\varphi }$, ω_η = ω_ζ = 0.

Next, the Lagrange’s equations are used to derive the governing equations of the motion of the harvesting system.

$$\frac{\partial }{\partial t}\left(\frac{\partial L}{\partial \dot{\mu }}\right)-\left(\frac{\partial L}{\partial \mu }\right)=\raisebox{1ex}{$\delta W$}\!\left/ \!\raisebox{-1ex}{$\delta \mu $}\right.$$

(2)

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).

$$T_b=\frac{1}{2}{\rho }_b\left[{\displaystyle {\int }_{V_b}\left({\dot{v}}^2+{\dot{w}}^2\right)d{V}_b}+{\displaystyle {\int }_0^L(I_{\xi }{{\omega }_{\xi }}^2)ds}\right]$$

(3)

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:

$$I_{\xi }=\raisebox{1ex}{$\left(b^2+h^2\right)$}\!\left/ \!\raisebox{-1ex}{$12$}\right.$$

where h is the thickness of the beam.

$$T_c=\frac{1}{2}{M}_c{\left[{\left(\dot{v}+{\dot{v}}^{\prime }\frac{L_c}{2}\right)}^2+{\left(\dot{w}+{\dot{w}}^{\prime }\frac{L_c}{2}\right)}^2\right]|}_{s=L}+\frac{1}{2}{J_{c\xi }{{\omega }_{\xi }}^2|}_{s=L}$$

(4)

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:

$$\begin{gathered} J_{c\xi }=\frac{1}{2}{M}_c\left(\frac{D^2}{4}+\frac{d^2}{4}\right)-\frac{1}{2}\pi {\rho }_c\frac{d^2}{4}{L}_c{e}^2 \\ {T}_f=\frac{1}{2}{M}_f\left[{\left({\dot{v}+{\dot{v}}^{\prime }\frac{2L_c-L_f}{2}|}_{s=L}\right)}^2+{\left({\dot{w}+{\dot{w}}^{\prime }\frac{2L_c-L_f}{2}|}_{s=L}\right)}^2\right]+\frac{1}{2}{J}_{f\xi }{{\omega }_{\xi }}^2{|}_{s=L} \end{gathered}$$

(5)

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:

$$M_f=\raisebox{1ex}{$C_M{\rho }_f\pi D^2{L}_f$}\!\left/ \!\raisebox{-1ex}{$4$}\right.$$

where C_M = 1 is the fluid-added mass coefficient.

$$J_{f\xi }=\raisebox{1ex}{$M_f{D}^2$}\!\left/ \!\raisebox{-1ex}{$8$}\right.$$

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.

$$U_i=\frac{1}{2}{\displaystyle {\int }_{V_s}({\sigma }^s{\epsilon }^s+{\tau }^s{\gamma }^s)}d{V}_s+\frac{1}{2}{\displaystyle {\int }_{V_p}\left[{\sigma }^p{\epsilon }^p+{\tau }^p{\gamma }^p\right]}d{V}_p$$

(6)

where the strains and stresses in the substrate and PZT layers are:

$$\begin{array}{lll}{\epsilon }_{11}=z\varphi {\eta }^{″}-y{\eta }^{″}& {\gamma }_{12}=2{\epsilon }_{12}=-z{\varphi }^{\prime }& {\gamma }_{13}=2{\epsilon }_{13}=y{\varphi }^{\prime }\\ {{\sigma }_{11}}^s=E_s{\epsilon }_{11}& {{\sigma }_{12}}^s=2G_s{\epsilon }_{12}& {{\sigma }_{13}}^s=2G_s{\epsilon }_{13}\\ {{\sigma }_{12}}^p=2G_p{\epsilon }_{12}& {{\sigma }_{13}}^p=2G_p{\epsilon }_{13}& \end{array}$$

(7)

Due to the poling of the PZT layer, the strain and electric displacement vector are governed by the constitutive equations, given by:

$$\begin{array}{l}{{\sigma }_{11}}^p=E_p{\epsilon }_{11}-e_{31}{E}_2\\ D_2=e_{31}{\epsilon }_{\eta }+{\in }_{11}{E}_2\end{array}$$

(8)

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₂ is the electric field intensity in the poling direction, and V(t) is the voltage between the two piezoelectric electrodes.

$$\begin{gathered} E_2(t)=V(t)/h_p \\ {U}_V=-\frac{1}{2}{\displaystyle {\int }_{V_p}{E}_2{D}_2}d{V}_p \end{gathered}$$

(9)

where V_p is the volume of the PZT layer.

$$U_g=M_{c}g{\frac{\partial x(s,t)}{\partial s}|}_{{}_{s=L}}\left[x(s,t)+\frac{L_c}{2}{\frac{\partial x(s,t)}{\partial s}|}_{{}_{s=L}}\right]+{\displaystyle {\int }_{V_b}{\rho }_b}gx(s,t)\frac{\partial x(s,t)}{\partial s}d{V}_b$$

(10)

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.

$$\delta W_R=-V(t)\delta Q(t)$$

(11)

where Q(t) is the charge between the piezoelectric electrodes.

$$\delta W_F=F\left(t\right)\delta \left[{v\left(s,t\right)|}_{s=L_s}+{v{\left(s,t\right)}^{\prime }|}_{s=L_s}({\mathrm{L}}_{\mathrm{c}}-\frac{L_{\mathrm{f}}}{2})\right]+F\left(t\right)l_d\delta \varphi$$

(12)

where l_d is the eccentric distance of mass and given by:

$$l_d=\frac{d^2}{D^2-d^2}e$$

and F(t) is the lift force. It is generated when the eccentric cylinder is subjected to the incoming fluid [22], as:

$$F\left(t\right)=\frac{C_L\left(t\right){\rho }_{f}D{U_0}^2{L}_f}{2}$$

(13)

where U₀ 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.

$$\ddot{q}+\epsilon {\omega }_f\left(q^2-1\right)\dot{q}+{\omega }_f{}^{2}q=f$$

(14)

where

$$q(t)=2C_L(t)/C_L{}_0$$

$${\omega }_f=\raisebox{1ex}{$2\pi S_t{U}_0$}\!\left/ \!\raisebox{-1ex}{$D$}\right.$$

where C_L0 = 0.3 [32]. And f is the nondimensional force of vortex [22], given by:

$$f=(\frac{A}{D})\frac{d^2\overline{x}}{d{t}^2}=(\frac{A}{D})\left[{\frac{{\partial }^{2}x\left(s,t\right)}{\partial t^2}|}_{s=L_s}+\left(L_c-\frac{L_f}{2}\right){\frac{{\partial }^{3}x\left(s,t\right)}{\partial s\partial t^2}|}_{s=L_s}+{\frac{{\partial }^2\varphi \left(s,t\right)}{\partial t^2}|}_{s=L_s}\right]$$

(15)

The virtual work of the buoyancy of fluid δW_b is:

$$\delta W_b={\displaystyle \underset{V_w}{\int }{\rho }_{f}g{\frac{\partial x(s,t)}{\partial s}|}_{{}_{s=L_s}}\delta \left[x(s,t)+\left(L_c-\frac{L_f}{2}\right){\frac{\partial x(s,t)}{\partial s}|}_{{}_{s=L_s}}\right]}d{V}_f$$

(16)

The virtual work of fluid-added damping δW_cf is:

$$\begin{array}{l}\delta W_{c_f}=c_f{\frac{\partial v\left(s,t\right)}{\partial t}|}_{s=L_s}\left(L_c-\frac{L_f}{2}\right)\delta \left[v\left(s,t\right){|}_{{}_{s=L_s}}\right]\\ +c_f{\frac{\partial w\left(s,t\right)}{\partial t}|}_{s=L_s}\left(L_c-\frac{L_f}{2}\right)\delta \left[w\left(s,t\right){|}_{{}_{s=L_s}}\right]+c_f{\omega }_{\xi }\frac{D}{2}\delta {\varphi \left(s,t\right)|}_{s=L_S}\end{array}$$

(17)

where

$$c_f=\gamma {\omega }_f{\rho }_f{D}^2$$

where γ = 0.8 [22] is the fluid-added damping coefficient.

$$\delta W_c=-{\displaystyle {\int }_0^L{c}_m\left[\dot{v}\delta v\left(s,t\right)+\dot{w}\delta w\left(s,t\right)\right]}ds$$

(18)

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

$$v\left(s,t\right)={\displaystyle \sum _{i=1}^n{\phi }_{vi}\left(s\right)}{r}_i\left(t\right)$$

(19)

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

(20)

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

(21)

where φ(s) and r(t) are the mode shapes and the model coordinates, respectively. They can be expressed as:

$$\begin{array}{l}{\phi }_v{}_i\left(s\right)=A_1\sin {\lambda }_1{}_{i}s+B_1\cos {\lambda }_{1i}s+C_1\sinh {\lambda }_{1i}s+D_1\cosh {\lambda }_{1i}s\\ {\phi }_w{}_i\left(s\right)=A_2\sin {\lambda }_2{}_{i}s+B_2\cos {\lambda }_2{}_{i}s+C_2\sinh {\lambda }_2{}_{i}s+D_2\cosh {\lambda }_2{}_{i}s\\ {\phi }_{\varphi }{}_i\left(s\right)=A_3\sin {\lambda }_3{}_{i}s+B_3\cos {\lambda }_3{}_{i}s\end{array}$$

(22)

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:

$$\begin{array}{l}{\lambda }_1{}_i=\sqrt[4]{\frac{\rho E{I}_{\eta }S}{\rho E{I}_{\zeta }S}}{\lambda }_2{}_i\\ {\lambda }_1{{}_i}^4\frac{E{I}_{\zeta }}{\rho S}={\lambda }_3{{}_i}^2\sqrt{\frac{G{I}_p}{J_{c\xi }}}\end{array}$$

(23)

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:

$$\begin{array}{l}E{I}_{\eta }=\frac{b\left(E_s\left(h_b{}^3-h_a{}^3\right)+E_p\left(h_c{}^3-h_b{}^3\right)\right)}{3}\\ E{I}_{\zeta }=\frac{b^3\left(E_s\left(h_b-h_a\right)+E_p\left(h_c-h_b\right)\right)}{12}\end{array}$$

$$G{I}_p=G_s\left[\frac{h_s{b}^3}{12}+\frac{b(h_b{}^3-h_a{}^3)}{3}\right]+G_p\left[\frac{h_p{b}^3}{12}+\frac{b(h_c{}^3-h_b{}^3)}{3}\right]$$

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.

$$h_0=\frac{\left(h_s+h_p\right)E_p{h}_p}{2\left(E_p{h}_p+E_s{h}_s\right)}+\raisebox{1ex}{$h_s$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

as follows:

$$h_a=-h_0{h}_b=h_s-h_0{h}_c=\left(h_s+h_p\right)-h_0$$

According to the work condition, their associated boundary conditions are used to solve the mode shapes above.

$$v(0,t)=0{\frac{\partial v}{\partial s}|}_{(0,t)}=0w(0,t)=0{\frac{\partial w}{\partial s}|}_{(0,t)}=0\varphi (0,t)=0$$

(24)

$$\begin{array}{l}E{I}_{\eta }{\frac{{\partial }^{2}v}{\partial s^2}|}_{s=L_s}+\frac{L_c}{2}{M}_c{\left(\frac{{\partial }^{2}v}{\partial t^2}+\frac{{\partial }^{3}v}{\partial s\partial t^2}\frac{L_c}{2}\right)|}_{s=L_s}\\ +\left(L_c-\frac{L_f}{2}\right)M_f{\left[\frac{{\partial }^{2}v}{\partial t^2}+\frac{{\partial }^{3}v}{\partial s\partial t^2}\left(L_c-\frac{L_f}{2}\right)\right]|}_{s=L_s}+M_c{l}_d\frac{L_c}{2}{\frac{{\partial }^2\varphi }{\partial t^2}|}_{s=L_s}=0\end{array}$$

(25)

$$\begin{array}{l}E{I}_{\eta }{\frac{{\partial }^{3}v}{\partial s^3}|}_{s=L_s}-M_c{\left(\frac{{\partial }^{2}v}{\partial t^2}+\frac{{\partial }^{3}v}{\partial s\partial t^2}\frac{L_c}{2}\right)|}_{s=L_s}\\ -M_f{\left[\frac{{\partial }^{2}v}{\partial t^2}+\frac{{\partial }^{3}v}{\partial s\partial t^2}\left(L_c-\frac{L_f}{2}\right)\right]|}_{s=L_s}-M_c{l}_d{\frac{{\partial }^2\varphi }{\partial t^2}|}_{s=L_s}=0\end{array}$$

(26)

$$\begin{array}{l}E{I}_{\zeta }{\frac{{\partial }^{2}w}{\partial s^2}|}_{s=L_s}+\frac{L_c}{2}{M}_c{\left(\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^{3}w}{\partial s\partial t^2}\frac{L_c}{2}\right)|}_{s=L_s}\\ +\left(L_c-\frac{L_f}{2}\right)M_f{\left[\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^{3}w}{\partial s\partial t^2}\left(L_c-\frac{L_f}{2}\right)\right]|}_{s=L_s}=0\end{array}$$

(27)

$$\begin{array}{l}E{I}_{\zeta }{\frac{{\partial }^{3}w}{\partial s^3}|}_{s=L_s}-M_c{\left(\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^{3}w}{\partial s\partial t^2}\frac{L_c}{2}\right)|}_{s=L_s}\\ -M_f{\left[\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^{3}w}{\partial s\partial t^2}\left(L_c-\frac{L_f}{2}\right)\right]|}_{s=L_s}=0\end{array}$$

(28)

$$G{I}_p{\frac{\partial \varphi }{\partial s}|}_{s=L_s}+J_{c\xi }{\frac{{\partial }^2\varphi }{\partial t^2}|}_{s=L_s}+M_c{l}_d{\left(\frac{{\partial }^{2}v}{\partial t^2}+\frac{{\partial }^{3}v}{\partial s\partial t^2}\frac{L_c}{2}\right)|}_{s=L_s}=0$$

(29)

By substituting Equations (19)–(21) into Equations (24)–(29), the equations of boundary conditions expressed by mode shapes can be written as:

$${\phi }_{vi}\left(0\right)=0{\phi }_{vi}^{\prime }\left(0\right)=0{\phi }_{wi}\left(0\right)=0{\phi }_{wi}^{\prime }\left(0\right)=0{\phi }_{\varphi i}\left(0\right)=0$$

(30)

$$\begin{array}{l}E{I}_{\eta }{\phi }_{vi}^{″}\left(L_s\right)-\frac{L_c}{2}{M}_c{{\omega }_i}^2\left({\phi }_{vi}\left(L_s\right)+{\phi }_{vi}^{\prime }\left(L_s\right)\frac{L_c}{2}\right)\\ -\left(L_c-\frac{L_f}{2}\right)M_f{\omega }_i{}^2\left[{\phi }_{vi}\left(L_s\right)+{\phi }_{vi}^{\prime }\left(L_s\right)\left(L_c-\frac{L_f}{2}\right)\right]-M_c{\omega }_i{}^2{l}_d\frac{L_c}{2}{\phi }_{\varphi i}\left(L_s\right)=0\end{array}$$

(31)

$$\begin{array}{l}E{I}_{\eta }{\phi }_{vi}{}^{\prime \prime \prime }\left(L_s\right)+M_c{\omega }_i{}^2\left({\phi }_{vi}\left(L_s\right)+{\phi }_{vi}^{\prime }\left(L_s\right)\frac{L_c}{2}\right)\\ +M_f{\omega }_i{}^2\left[{\phi }_{vi}\left(L_s\right)+{\phi }_{vi}^{\prime }\left(L_s\right)\left(L_c-\frac{L_f}{2}\right)\right]+M_c{\omega }_i{}^2{l}_d{\phi }_{\varphi i}\left(L_s\right)=0\end{array}$$

(32)

$$\begin{array}{l}E{I}_{\zeta }{\phi }_{wi}{}^{″}\left(L_s\right)-\frac{L_c}{2}{M}_c{\omega }_i{}^2\left({\phi }_{wi}\left(L_s\right)+{\phi }_{wi}^{\prime }\left(L_s\right)\frac{L_c}{2}\right)\\ -\left(L_c-\frac{L_f}{2}\right)M_f{\omega }_i{}^2\left[{\phi }_{wi}\left(L_s\right)+{\phi }_{wi}^{\prime }\left(L_s\right)\left(L_c-\frac{L_f}{2}\right)\right]=0\end{array}$$

(33)

$$\begin{array}{l}E{I}_{\zeta }{\phi }_{wi}{}^{\prime \text{}\prime \text{}\prime }\left(L_s\right)+M_c{\omega }_i{}^2\left({\phi }_{wi}\left(L_s\right)+{\phi }_{wi}^{\prime }\left(L_s\right)\frac{L_c}{2}\right)\\ +M_f{\omega }_i{}^2\left[{\phi }_{wi}\left(L_s\right)+{\phi }_{wi}^{\prime }\left(L_s\right)\left(L_c-\frac{L_f}{2}\right)\right]=0\end{array}$$

(34)

$$G{I}_p{\phi }_{\varphi i}{}^{\prime }\left(L_s\right)-J_{c\xi }{\omega }_i{}^2{\phi }_{\varphi i}\left(L_s\right)-M_c{l}_d{\omega }_i{}^2\left({\phi }_{vi}\left(L_s\right)+{\phi }_{vi}^{\prime }\left(L_s\right)\frac{L_c}{2}\right)=0$$

(35)

Normalizing the eigenfunctions using the following orthogonality conditions of mass and stiffness, the relations between different mode shapes are obtained.

$$\begin{array}{l}{\rho }_b{\displaystyle {\int }_{V_b}({\phi }_{vi}\left(s\right){\phi }_{vj}\left(s\right)+{\phi }_{wi}\left(s\right){\phi }_{wj}\left(s\right))d{V}_b}\\ +M_c\left({\phi }_{vi}\left(L_s\right)+\frac{L_c}{2}{\phi }_{vi}{}^{\prime }\left(L_s\right)\right)\left({\phi }_{vj}\left(L_s\right)+\frac{L_c}{2}{\phi }_{vj}{}^{\prime }\left(L_s\right)\right)\\ +M_c\left({\phi }_{wi}\left(L_s\right)+\frac{L_c}{2}{\phi }_{wi}{}^{\prime }\left(L_s\right)\right)\left({\phi }_{wj}\left(L_s\right)+\frac{L_c}{2}{\phi }_{wj}{}^{\prime }\left(L_s\right)\right)\\ +M_f\left({\phi }_{vi}\left(L_s\right)+\frac{\left(2L_c-L_f\right)}{2}{\phi }_{vi}{}^{\prime }\left(L_s\right)\right)\left({\phi }_{vj}\left(L_s\right)+\frac{\left(2L_c-L_f\right)}{2}{\phi }_{vj}{}^{\prime }\left(L_s\right)\right)\\ +M_f\left({\phi }_{wi}\left(L_s\right)+\frac{\left(2L_c-L_f\right)}{2}{\phi }_{wi}{}^{\prime }\left(L_s\right)\right)\left({\phi }_{wj}\left(L_s\right)+\frac{\left(2L_c-L_f\right)}{2}{\phi }_{wj}{}^{\prime }\left(L_s\right)\right)\\ +\left(J_{c\xi }+J_{f\xi }\right)\left[{\phi }_{\varphi i}\left(L_s\right)+{\phi }_{wi}{}^{\prime }\left(L_s\right){\phi }_{vi}{}^{\prime }\left(L_s\right)\right]\left[{\phi }_{\varphi j}\left(L_s\right)+{\phi }_{wj}{}^{\prime }\left(L_s\right){\phi }_{vj}{}^{\prime }\left(L_s\right)\right]\\ ={\delta }_{ij}\end{array}$$

(36)

$${\displaystyle {\int }_0^L{\phi }_{vi}{}^{″}E{I}_{\eta }{\phi }_{vj}{}^{″}}ds+{\displaystyle {\int }_0^L{\phi }_{\varphi i}{}^{\prime }G{I}_p{\phi }_{\varphi j}{}^{\prime }}ds={\delta }_{ij}{\omega }_i{}^2$$

(37)

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.

By substituting Equations (19)–(21) and the orthogonality conditions (36) and (37) into the kinetic energy Equations (3)–(5), potential energy Equations (6), (9) and (10), and the virtual work Equations (11), (12), (16)–(18), the energy equations and virtual work equations are rewritten as:

The kinetic energy equation T:

$$T=\raisebox{1ex}{${\dot{r}}^2\left(t\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

(38)

The potential energy equation U:

$$U=H_{1}Vr\left(t\right)+\frac{1}{2}{\omega }_i{}^2{r}^2\left(t\right)-\frac{1}{2}{H}_1\mathrm{V}r\left(t\right)+H_2{\mathrm{V}}^2+{\mathrm{H}}_3{r}^2(t)$$

(39)

where H_n is the coefficient of the variables, given by:

$$H_1=-e_{31}\frac{b\left(h_b{}^2-h_a{}^2\right)}{2h_p}{\displaystyle {\int }_0^{L_p}{\phi }_v{}^{″}\left(s\right)}ds$$

$$H_2=-\frac{{\epsilon }_{11}b{L}_p}{2h_p}$$

$$\begin{array}{l}{\mathrm{H}}_3=M_{c}g\left[{\phi }_v\left(L_s\right){\phi }_v{}^{\prime }\left(L_s\right)+\frac{L_c}{2}{\phi }_v{{}^{\prime }}^2\left(L_s\right)+{\phi }_w\left(L_s\right){\phi }_w{}^{\prime }\left(L_s\right)+\frac{L_c}{2}{\phi }_w{{}^{\prime }}^2\left(L_s\right)\right]\\ +{\rho }_b{S}_{b}g{\displaystyle {\int }_0^{Ls}\left({\phi }_v{\phi }_v{}^{\prime }+{\phi }_w{\phi }_w{}^{\prime }\right)ds}\end{array}$$

The virtual work equation:

$$\frac{\delta W}{\delta r(t)}=H_{4}F\left(t\right)+H_{5}r(t)+H_6\dot{r}(t)$$

(40)

where the coefficient H_n are expressed as:

$$H_4=\left[{\phi }_v\left(L_s\right)+{\phi }_v^{\prime }\left(L_s\right)\left(L_c-\frac{L_{\mathrm{f}}}{2}\right)\right]+l_d{\phi }_{\varphi }\left(L_s\right)$$

$$\begin{array}{l}{H}_5={\rho }_{w}g\frac{\pi D^2}{4}{\displaystyle {\int }_0^{L_{\mathrm{f}}}\left[{\phi }_v\left(s\right){\phi }_v{}^{\prime }\left(L_s\right)+\left(L_c-\frac{L_{\mathrm{f}}}{2}\right){\phi }_v{{}^{\prime }}^2\left(s\right)\right]ds}\\ +{\rho }_{w}g\frac{\pi D^2}{4}{\displaystyle {\int }_0^{L_{\mathrm{f}}}\left[{\phi }_w\left(L_s\right){\phi }_w{}^{\prime }\left(s\right)+\left(L_c-\frac{L_{\mathrm{f}}}{2}\right){\phi }_w{{}^{\prime }}^2\left(s\right)\right]ds}\end{array}$$

$$H_6=-c_f{\phi }_v{}^2\left(L_s\right)-c_f{\phi }_w{}^2\left(L_s\right)-c_f\frac{D}{2}{\phi }_{\varphi }{}^2\left(L_s\right)-{\displaystyle {\int }_0^L{c}_m\left[{\phi }_v{}^2\left(s\right)+{\phi }_w{}^2\left(s\right)\right]}ds$$

Now, the Lagrange function which is the subtraction of kinetic energy Equation (38) and potential energy Equation (39) is established and expressed as:

$$\begin{array}{l}L=T-U\\ =\frac{1}{2}{\dot{r}}^2\left(t\right)-\left(\frac{1}{2}{\omega }_i{}^2+{\mathrm{H}}_3\right)r^2\left(t\right)-\frac{1}{2}{H}_{1}Vr\left(t\right)-H_2{\mathrm{V}}^2\end{array}$$

(41)

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:

$$\ddot{r}\left(t\right)+\upsilon \dot{r}\left(t\right)+\lambda r(t)+\theta V=\kappa F\left(t\right)$$

(42)

$$\frac{V\left(t\right)}{R}+C_p\frac{dV\left(t\right)}{dt}-\theta \dot{r}\left(t\right)=0$$

(43)

where the coefficients of the state variables are expressed as:

$$\upsilon =-H_6,\lambda ={\omega }_i{}^2+2H_3-H_5,\kappa =H_4,C_p=-2H_2,\mathrm{and}\theta =\raisebox{1ex}{$H_1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

Substituting the Equations (19)–(21) into Equation (15) the nondimensional force of vortex f is rewritten as:

$$f=H_7\ddot{r}(t)$$

(44)

where the coefficient H₇ is given by:

$$H_7=(\frac{A}{D})\left[{\phi }_v\left(L_s\right)+(L_c-\frac{L_{\mathrm{f}}}{2}){\phi }_v^{\prime }\left(L_s\right)+{\phi }_w\left(L_s\right)+(L_c-\frac{L_{\mathrm{f}}}{2}){\phi }_w^{\prime }\left(L_s\right)+{\phi }_{\varphi }\left(L_s\right)\right]$$

According to Equations (42)–(44), a space state equation of the energy harvesting system is established, and expressed as:

$$\dot{X}=\left[\begin{array}{c}{X}_2\\ \frac{\kappa C_{L0}{\rho }_{w}D{U}^2{L}_c}{4}{X}_4-\upsilon X_2-\lambda X_1-\theta X_3\\ \frac{\theta }{C_p}{X}_2-\frac{X_3}{C_{p}R}\\ X_5\\ \ddot{q}\end{array}\right]$$

(45)

where X is the space state variable, written as:

$$X=\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\\ X_5\end{array}\right]=\left[\begin{array}{c}r\\ \dot{r}\\ V\\ q\\ \dot{q}\end{array}\right]$$

Equation (45) can be numerically solved by using the solver of ode45 in MATLAB (MathWorks, Inc^®, Natick, MA, USA).

3. Numerical Analysis and Experimental Validation

3.1. 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.

3.2. Numerical and Experimental Analysis

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 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. Configurations of the eccentric cylinder.

Configurations	2006	2008	2010	2508	2510	2512	3010	3012	3014
D (mm)	20	20	20	25	25	25	30	30	30
d (mm)	6	8	10	8	10	12	10	12	14
e (mm)	5	4	3	6.25	5	4	7.5	6	5
ld (mm)	0.495	0.762	0.952	0.713	0.925	1.198	0.938	1.142	1.392

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₀ 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₀, while U₀ continues to increase, P decreases. The optimal velocity is about 0.339–0.366 m/s for the nine configurations. As the flow velocity U₀ 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 of different configurations of the eccentric cylinder when b = 20 mm.

Configurations	2000	2006	2008	2010	2500	2508
Pmax (mW)	0.1831	0.2663	0.281	0.3267	0.2139	0.2763
Configurations	2510	2512	3000	3010	3012	3014
Pmax (mW)	0.3254	0.3273	0.1764	0.249	0.2624	0.2866

. 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.

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.

4. Conclusions

In this paper, a vortex-induced vibration piezoelectric energy harvester undergoing bending–torsion vibrations by the tip eccentric cylinder was designed and fabricated. The water-to-electricity conversion characteristics in low-speed water flow was studied. A distributed parameter beam model was established based on Euler–Bernoulli beam theory. The mode shapes and model coordinates were used as basic functions to represent the displacements of harvester system in a Galerkin procedure. Energy-governing coupled equations of the harvesting system were derived by using the Lagrange’s equations. Numerical models were used to analyze and predict the best configurations and optimal work environments for the harvesters by calculating the output voltage. Experiments have been designed and conducted. The correctness of the numerical models was verified. A harvester with D = 20 mm, d = 10 mm, e = 3 mm has a maximum power of 0.3978 mW. The eccentric hole setting of the harvester increases its harvesting power as much as 99%. This paper provides a practical significance in an efficient harvester’s design and manufacture.