Analysis of clamped-clamped piezoelectric energy harvester under vortex induced vibration considering the stretching effect

Abstract

This paper investigates the middle layer stretching effect on the dynamic behavior of an energy harvester clamped-clamped beam. Two foam cylinders that are attached together as a dumbbell are mounted on the middle part of the beam for vortex-induced vibrations. The piezoelectric patch collects the electrical energy produced by vortex induced vibration. In this analysis the coupled differential equations governing on the structure oscillation, harvested voltage and fluid lift force are established applying the Hamilton’s principle, Gauss law and wake oscillator model. The obtained differential equations are discretized using Galerkin method, and solved by both numerical and analytical perturbation methods. The results demonstrate that middle layer stretching effect has a significant effect on the lock-in domain, and the output electric voltage must be evaluated by considering this effect. It has been shown that for large values of cylinder diameter the difference between numerical and theoretical results increases due to increasing the middle layer stretching effect.

Keywords

Energy harvesting non-linear vibration piezoelectric vortex-induced-vibration middle layer stretching effect

1. Introduction

For two decades, researchers have focused on the electrical energy harvesting from vibrating systems to replace electrical batteries. One of the ways to harvest energy from vibrations is the use of piezoelectric effect. This effect says that when the mechanical strain is applied to a piezoelectric material, an electric potential difference is created in it (Covaci and Gontean, 2020). Consequently, when a layer of piezoelectric material is deposited or attached on a vibrating beam, it will act as an energy collector. Compared with other transducer of vibrations to electric energy such as electromagnetic and electrostatic transducers, it should be said that piezoelectric transducer has inherent electromechanical coupling and high-power density. Therefore, they have received more attention in the scope of electrical energy production (Dahiya and Valle, 2013; Firouzi et al., 2023).

When the air flow passes over a bluff body attached to the tip of a cantilever beam, the beam will be actuated due to the aerodynamic forces. The amplitude and frequency of the actuating aerodynamic forces applied on the beam depend on the flow pattern behind of the bluff body (Sarpkaya, 1979). In a range of flow speed, the periodic vortexes are created behind of the cylinder. In this domain a vortex is generated from the top of the cylinder, and then another vortex starts from its bottom. The periodic natural of the vortexes causes a periodic lift force on the system. It causes the periodic oscillation of system that is called vortex induced vibration (VIV) (Blevins, 1977). The frequency of the created periodic vortex behind of the cylinder depends on the interaction between the structure and air flow. When the shedding of periodic vortexes is started, at first the frequency of periodic vortex shedding is independent of the structure, and its frequency increases by increasing of the value of flow speed. In the following by increasing the value of flow speed and nearing the frequency of the vortex shedding to the natural frequency of the system, a considerable interaction is presented between flow and structure. In this domain of the flow speed the frequency of vortex periodic shedding will be locked on the natural frequency and causes a large amplitude for the oscillations that is called lock-in domain (Williamson and Govardhan, 2004, 2008). After this domain that is, by more increases of the value of flow speed the amplitude of the oscillation decreases. Unlike resonance, lock-in or synchronization is a nonlinear phenomenon that occurs in a band of frequency that is, in a band of flow speed (Williamson and Govardhan, 2004).

It is clear that an energy harvester piezoelectric beam under vortex induced vibration has a complex aero-elastic and electro-mechanic coupling. It means that it is difficult to present a complete theoretical model for it. Therefore, many researchers focused on experimental methods. However, the experimental works require spending money, so a major part of works have been presented using semi-empirical models proposed for interaction between structure and fluid. The most important of them is the Facchinetti et al. oscillator model, which is based on the Van der Pol equation. In this modeling the Van der Pol equation including motion acceleration term at its right-hand side is coupled with the equation of motion. The dependent variable in Van der Pol equation is in direct relation to the lift coefficient (Facchinetti et al., 2004; Qu and Metrikine, 2020). Before Facchinetti model the coupling term has been considered as expressions of displacement and velocity (Skop and Balasubramanian, 1997; Skop and Griffin, 1973; Skop and Luo, 2001). The important work based on the experimental and semi-empirical method are reviewed in follow.

Mehmood et al. (2013) modeled the VIV piezoelectric energy harvester as a mass-spring-damper system and considered the effect of the piezoelectric layer as a load resistance. The aerodynamic forces due to vortex shedding have been obtained by using CFD simulations where the incompressible continuity and unsteady Navier–Stokes equations have been solved using an acceleration reference frame. The results showed that the load resistance has a significant effect on the lock-in area in such a way that with increasing it the lock-in area become wider. Abdelkefi et al. (2012) theoretically investigated the effect of various parameters on the mechanical and electrical responses of a VIV harvester modeled as a mass-spring-damper system with an electrical resistance. They found that electrical load resistance of the piezoelectric layer shifts the onset of the lock-in range to higher fluid velocities. The nonlinearity of electric load resistance causes hardening behavior in displacement domains and output power.

Akaydin et al. (2010) studied a cantilever piezoelectric beam excited by vortex behind a cylinder. The cylinder has been considered in a distance from the free end of cantilever beam. They used Navier-Stokes equations for evaluating the fluid force, and showed that numeric results have a good agreement with the experimental results. Jia et al. (2018) presented a horizontal piezoelectric energy harvester Euler Bernoulli beam where a cylinder has been attached to it parallel its length. They used the semi empirical van der Pol wake oscillator model based on Facchinetti model for simulating the aerodynamic force induced by the vortex. Gao et al. (2013) showed that turbulent airflows generates higher power and electrical voltage rather than the laminar flow in the configuration studied in Jia et al. (2018). Dai et al. (2016) compared the behavior of cantilever piezoelectric energy harvester when the cylindrical bluff body are installed at the tip of cantilever beam with different orientations (bottom, top, horizontal, and vertical). Their experiments results showed that the synchronization regions of the bottom, top, and horizontal configurations are almost the same at low wind speeds, and the vertical configuration has the highest wind speed for synchronization with the largest harvested power. The effect of hybrid bluff bodies on VIV energy harvesting performance has been investigated in Li et al. (2023). They considered the hybrid bluff body as a combination of two kinds of bluff elements, that is, O-shaped cylinder and D-shaped prism. The vortex induced vibration of two bimorph piezoelectric cantilevers attached to a single cylinder in their tips has been studied in other research (Zhang and Wang, 2016).

Many works have been done to consider the piezoelectric energy harvester when a cubic bluff body is attached to the tip of the beam. The aerodynamic force in this category of works is due to galloping phenomena (Ewere and Wang, 2014; Rezaei and Talebitooti, 2019; Seyed-Aghazadeh et al., 2020; Wang et al., 2020; Zhao et al., 2012, 2016, 2019). Unlike VIV that exhibits large amplitudes only when the vortex-shedding frequency is near the structure’s natural frequency (lock-in or synchronization), the galloping phenomenon exhibits large amplitudes after a critical speed. In this phenomenon a series of steady limit cycles will be generated at each flow speed value larger than critical speed. The combination of both galloping and VIV phenomena for energy harvesting has been proposed recently in some works (Sun et al., 2019; Yang and He, 2019). Li et al. considered the interaction between vortex-induced vibration (VIV) and galloping for a cantilever piezoelectric energy harvester considering the effects of geometrical nonlinearity (Li et al., 2022a).

Recently magnetic-coupling piezoelectric energy harvester (MCPEH) based on the vortex-induced vibration has been proposed (Hou et al., 2020; Li et al., 2022b). Ma et al. indicated that the working wind speed range can be widened and adjusted using a magnetic-coupled tri-stable configuration (Ma et al., 2023). The harvester consists of two working units, a piezoelectric unit, and an electromagnetic unit under combined vortex-induced and base excitations has been investigated in Hou et al. (2022). It has been shown that combination of the excitation improves the performances of the harvester. Magnetically coupling bending-torsion piezoelectric energy harvester based on vortex-induced vibration has been investigated in Sui et al. (2022). The literature review shows that the main element of piezoelectric energy harvesters are often a cantilever beam with a bluff body in its tip. Therefore, the nonlinear effect due to middle layer stretching effect has been neglected in the previous works (Abdelkefi et al., 2012; Akaydin et al., 2010; Dai et al., 2016; Ewere and Wang, 2014; Gao et al., 2013; Hou et al., 2020, 2022; Jia et al., 2018; Li et al., 2022a, 2022b, 2023; Ma et al., 2023; Mehmood et al., 2013; Rezaei and Talebitooti, 2019; Seyed-Aghazadeh et al., 2020; Skop and Griffin, 1973; Sui et al., 2022; Sun et al., 2019; Wang et al., 2020; Yang and He, 2019; Zhang et al., 2022; Zhang and Wang, 2016; Zhao et al., 2012, 2016, 2019).

Recently a configuration comprised of a horizontal main beam and two vertical side beams has been proposed by Su and Lin in order to collect wind energy in two horizontal and vertical directions (Su and Lin, 2020; Su and Wang, 2021). They assumed that the piezoelectric patches are bonded to the main beam, and two foam cylinders are attached to the center of the main beam as a dumbbell for vortex-induced vibrations. If the side beams stiffness of the configuration proposed in their work is assumed to be high then, the main beam will be similar to a clamped-clamped beam. One important phenomenon affecting on the dynamic behavior of clamped–clamped beam is the middle layer stretching effect (Rezaei and Zamanian, 2017; Younis and Nayfeh, 2003; Zamanian and Khadem, 2010a, 2010b). This phenomenon has been neglected in previous work (Su and Lin, 2020; Su and Wang, 2021). Therefore, this paper investigates the effect of stretching effect on the dynamic behavior of an energy harvester piezoelectric beam, where two foam cylinders are attached to the center of the main beam as a dumbbell for vortex-induced vibrations

It has been shown in Li et al. (2024); Mackowski and Williamson (2013) that nonlinear stiffness in a spring mass system under the vortex induced vibration causes a different pattern for lock-in domain. It increases the necessary of considering the mid-plane stretching effect on the amount of harvested energy when the configuration is a dumbbell shape cylinder mounted on the middle part of the clamped-clamped beam.

In this analysis the effect of vortex flow behind of the cylinders are considered based on wake oscillator model proposed by Facchinetti. The equation of motion coupled with wake oscillator model is obtained using the Hamilton’s principle. The Gauss law is used for obtaining the differential equation governing on the output voltage coupled with equation of motion. Firstly, free vibration mode shape of system are obtained, and then using them as comparison functions in Galerkin method, the coupled differential equations governing on the system are discretized. The discretized equations are solved by both numerical and analytical perturbation method. It can be said that in addition to considering the middle layer stretching effect as an innovation, this paper is an effort to present an analytical solution for the considered configuration which has not been presented in previous studies.

2. Modeling and formulation

The considered system is a fixed-fixed beam with a length of $l_{b}$ , which is shown in Figures 1 and 2. A pair of cylinders with diameter $D$ and length l_s in the form of a dumbbell is mounted on the middle part of the beam. The fluid flow parallel to the length of the beam is passed on the cylinders. A piezoelectric layer is mounted on each sides of the beam, which is located at distance $l_{1}$ from both ends of the beam. The width of the beam and piezoelectric layer are assumed to be equal. The interaction between the fluid and structure occurs due to passing flow over the cylinders, and system oscillates due to producing the vortex behind of the cylinders. The beam is a slender beam, and so the Euler-Bernoulli beam theory assumptions are used in derivation of the motion equations.

To describe the displacement of the structure, the fixed global coordinate system $syz$ shown in Figure 2 and the local coordinate system $\bar{s} \bar{y} \bar{z}$ shown in Figure 3(a) are defined. The displacement of the beam element with respect to the fixed coordinate systems is assumed to be $W$ and $U$ . It is assumed that displacement in “y” direction is negligible, and the rotation of the cylinder about s axis in negligible.

Considering the effect of middle-layer stretching, the strain of the beam and the piezoelectric layer in the cross-section can be expressed as (Nayfeh and Pai, 2004)

ε (\bar{z}) = e - κ \bar{z}

(1)

Where $e$ is the stretching strain of neutral axis of cross-section and $κ$ is the curvature around $\bar{y} axis$ , respectively. Considering Figure 3(b) and using Taylor series expansions, the description of $e$ and $κ$ may be written as follow (Nayfeh and Pai, 2004).

\begin{matrix} e = \sqrt{{(1 + U')}^{2} + {W'}^{2}} - 1 = U' + 1 / 2 {W'}^{2} - 1 / 2 U' {W'}^{2} + \dots, \\ κ = (\tan^{- 1} \frac{W'}{1 + U'})' = (W' - U' W' + {U'}^{2} W' - 1 / 3 {W'}^{3})' \end{matrix}

(2)

There is an electromechanical coupling in the piezoelectric materials. Here, the electrical displacement is assumed one-dimensional; thus, the relation between stress and strain in the beam and piezoelectric layer may be written as follow (Zamanian and Khadem, 2010a).

\begin{matrix} σ_{b} = E_{b} ε_{b}, \\ σ_{p} = E_{p} ε_{p} - E_{p} d_{31} \frac{V (t)}{h_{p}} \end{matrix}

(3)

Where $σ_{b}$ and $σ_{p}$ are the stress in cross section of the beam and piezoelectric layer, respectively, $d_{31}$ denotes the constant of the piezoelectric layer, and $V (t)$ is the output voltage harvested by the piezoelectric layer. Also $E_{b}$ and $E_{p}$ are the elasticity modulus of the beam and piezoelectric layer, respectively. As mentioned before the strain in the beam and piezoelectric element is related to $\bar{z}$ that is, the distance between neutral axis and the position of the element. As shown in Figure 4, the distance between the neutral axis and midline of the beam cross section is shown by ${\bar{z}}_{n}$ , and is evaluated setting the normal force in the cross section of system equal to zero (Mahmoodi and Jalili, 2007). After this operation one obtains:

{\bar{z}}_{n} = \frac{E_{p} h_{p} (h_{p} + h_{b})}{2 (E_{b} h_{b} + E_{p} h_{p})}

(4)

Where $h_{b}$ and $h_{p}$ are the thickness of the beam and piezoelectric layer, respectively. The motion equations of the system are derived using Hamilton’s principle that is,

δ \int_{0}^{l_{b}} (T - V_{c} + W_{f}) dt = 0

(5)

In this relation $T$ is the kinetic energy, $V_{c}$ is the potential energy of the system, and $W_{f}$ is the work of the external force on system. The kinetic energy $, T,$ may be written as follows assuming the mass of cylinder in the middle of the beam as a concentrated mass and representing it by $m_{T}$

T = 1 / 2 \int m (s) ({\dot{U}}^{2} + {\dot{W}}^{2}) d s + 1 / 2 \int m_{T} δ (s - l_{b} / 2) ({\dot{U}}^{2} + {\dot{W}}^{2}) d s

(6)

Where $U$ and $W$ are the longitudinal and lateral displacement of the beam element in the fixed coordinate systems. The dot and prime signs show the differential with respect to $s$ and $t$ in the above equation. The expression of $δ (s - \frac{l_{b}}{2})$ denotes the Dirac function at point of $s = \frac{l_{b}}{2}$ . The expression of $m (s)$ shows the mass of the system per unit length, and it is may be written as follow

m (s) = b_{b} (ρ_{b} h_{b} + (H_{l_{1}} - H_{l_{1} + l_{p}}) ρ_{p} h_{p})

(7)

Where $ρ_{b}$ and $ρ_{p}$ are the density of the beam and piezoelectric layer. The expression of $H_{l_{i}}$ is Heaviside function which is considered as follow

H_{li} = Heaviside function (s - l_{i}) = {\begin{matrix} 1 s \geq l_{i} \\ 0 s < l_{i} \end{matrix}

(8)

Also, the potential energy of the system may be written as follows

V_{c} = \int \int_{0}^{ε_{b}} σ_{b} d ε_{b} d v_{b} + \int \int_{0}^{ε_{p}} σ_{p} d ε_{p} d v_{p}

(9)

Where $d v_{b}$ and $d v_{p}$ are the volume element for the beam and piezoelectric layer. By substituting equations (1) and (3) into equation (9)

\begin{matrix} V_{c} = \frac{1}{2} \int_{0}^{l_{b}} c_{ξ} (s) e^{2} ds + \frac{1}{2} \int_{0}^{l_{b}} c_{η} (s) k^{2} ds - \int_{0}^{l_{b}} c_{α} (s) V (t) eds \\ + \int_{0}^{l_{b}} c_{γ} (s) V (t) kds \end{matrix}

(10)

In which

\begin{matrix} C_{ξ} (s) = (1 - H_{l_{1}}) E_{b} A_{b} + (H_{l_{1}} - H_{lp + l_{1}}) E_{p} A_{p} \\ + (H_{l_{1}} - H_{lp + l_{1}}) E_{b} A_{b} + (H_{lp + l_{1}}) E_{b} A_{b}, \\ C_{η} (s) = (1 - H_{l_{1}}) E_{b} I_{1}^{b} + (H_{l_{1}} - H_{lp + l_{1}}) E_{p} I_{2}^{p} \\ + (H_{l_{1}} - H_{lp + l_{1}}) E_{b} I_{2}^{b} + (H_{lp + l_{1}}) E_{b} I_{1}^{b}, \\ C_{α} (s) = (H_{l_{1}} - H_{lp + l_{1}}) \frac{E_{p} d_{31} A_{p}}{h_{p}}, \\ C_{γ} (s) = (H_{l_{1}} - H_{lp + l_{1}}) \frac{E_{p} d_{31} I_{1}^{p}}{h_{p}} \end{matrix}

(11)

and

\begin{matrix} A_{p} = b_{p} h_{p}, \\ I_{1}^{b} = \frac{1}{12} b_{b} h_{b}^{3}, \\ I_{2}^{b} = \frac{1}{12} b_{b} h_{b}^{3} + h_{b} b_{b} {\bar{z}}_{n}^{2}, \\ I_{2}^{p} = b_{p} [\frac{1}{3} (h_{p}^{3} + \frac{3}{2} h_{b} h_{p}^{2} + \frac{3}{4} h_{p} h_{b}^{2} + h_{p} {\bar{z}}_{n}^{2} - (h_{p}^{2} + h_{b} h_{p}) {\bar{z}}_{n}], \\ I_{1}^{p} = \frac{b_{p}}{2} (h_{p}^{2} + h_{b} h_{p} - 2 h_{p} {\bar{z}}_{n}) \end{matrix}

(12)

Where $b_{b}$ and $b_{p}$ are the width of the beam and piezoelectric layer. The work by external force may be written as below:

δ W_{f} = - \int_{0}^{l_{b}} \hat{C} \dot{W} d s δ W + \int_{0}^{l_{b}} F_{L} (t) δ (s - l_{b} / 2) d s δ W

(13)

Where $\hat{C}$ is the viscous damping coefficient and $F_{L} (t)$ is the fluid lift force applied to the cylinder. By considering equation (6) and substituting equation (2) into equation (10), neglecting the inertia effect along the longitudinal direction, the equation of motion in longitudinal directions will be as follow

(U' C_{ξ} (s) + 1 / 2 {W'}^{2} C_{ξ} (s))' - C_{α} (s)' V (t) = 0

(14)

Also, the equation of motion in transverse direction will be as follow

\begin{array}{l} m (s) \overset{\cdot\cdot}{W} + m_{T} (s) δ (s - l_{b} / 2) \overset{\cdot\cdot}{W} + \hat{C} δ (s - l_{b} / 2) \dot{W} + (C_{η} (s) W^{″})^{″} \\ = C_{ξ} (s) (W^{'} U^{″} + W^{″} U^{'} - 2 U^{'} U^{″} W^{'} - {U^{'}}^{2} W^{″} + \frac{3}{2} {W^{'}}^{2} W^{″}) \\ + {C^{'}}_{ξ} (s) (- W^{'} {U^{'}}^{2} + W^{'} U^{'} + \frac{1}{2} {W^{'}}^{3}) \\ + C_{α} V (t) (- W^{″} + W^{'} U^{″} + W^{″} U^{'}) \\ + C_{α}^{‣} V (t) (- W^{'} + W^{'} U^{'}) + C_{γ}^{‣} V (t) U^{″} \\ - {C^{″}}_{γ} V (t) + F_{L} δ (s - l_{b} / 2) \end{array}

(15)

It must be mentioned that the inertia effect in the longitudinal direction of the beam has been neglected since the natural frequency of longitudinal vibrations of the beam is higher than the natural frequency of transverse vibrations. Also, the nonlinear terms including bending stiffness coefficient have been neglected in comparison to the nonlinear terms including axial stiffness coefficient because the beam is considered to be thin. Equations (14) and (15) for the beam with a piezoelectric layer can be verified by comparing with the ones obtained by Zamanian and Khadem (Zamanian and Khadem, 2010a). According to equation (14)

U' = - 1 / 2 {W'}^{2} + \frac{C_{α} (s)}{C_{ξ} (s)} V (t) + \frac{C_{1} (t)}{C_{ξ} (s)}

(16)

Where $C_{1} (t)$ is the constant of integration. Integration from equation (14)

U (l_{b}) - U (0) = - \frac{1}{2} \int_{0}^{l_{b}} {W'}^{2} ds + \int_{0}^{l_{b}} \frac{C_{α} (s)}{C_{ξ} (s)} V (t) ds + \int_{0}^{l_{b}} \frac{C_{1} (t)}{C_{ξ} (s)} ds

(17)

By substituting boundary conditions $U (0) = 0$ and $U (l_{b}) = 0$ into equation (17), it result that

C_{1} (t) = \frac{\bar{ε} + \frac{1}{2} \int_{0}^{l_{b}} {W'}^{2} ds - V (t) \frac{E_{p} d_{31} A_{p}}{h_{p} (E_{b} A_{b} + E_{p} A_{p})} l_{p}}{\frac{l_{p}}{E_{b} A_{b} + E_{p} A_{p}} + \frac{(l_{b} - l_{p})}{E_{b} A_{b}}}

(18)

The motion equation with one depended variable are obtained by substituting $U^{'}$ from equation (16) into equation (15) as

\begin{array}{l} m (s) \overset{\cdot\cdot}{W} + m_{T} (s) δ (s - l_{b} / 2) \overset{\cdot\cdot}{W} + \hat{C} δ (s - l_{b} / 2) \dot{W} \\ + (C_{η} (s) W^{″})^{″} - C_{1} (t) W^{″} \\ = - {C_{y}}^{''} V (t) + F_{L} δ (s - l_{b} / 2) \end{array}

(19)

As mentioned before $F_{L} (t)$ is the fluid lift force applied to the cylinder, which is expressed as (Jia et al., 2018).

F_{L} = \frac{l_{s} C_{L} (t) ρ_{f} D {U_{f}}^{2}}{2}

(20)

Where $ρ_{f}$ is the mass density of the flow, and $U_{f}$ is the flow velocity. Also, $C_{L} (t)$ is the lift coefficient due to the vortex behind of the cylinder, which is expressed as

C_{L} (t) = \frac{q (t) C_{L 0}}{2}

(21)

Where $C_{L 0}$ is the reference lift coefficient on an unfixed cylinder undergoing vortex shedding. Variable $q (t)$ depends on the vortex flow behind the cylinder, and it is coupled with cylinder acceleration based on wake oscillator model proposed by Facchinetti et al. (2004).

\frac{d^{2} q}{d t^{2}} + ω_{f} η (q^{2} - 1) \frac{d q}{d t} + {ω_{f}}^{2} q = (\frac{K}{D}) \frac{\partial^{2} W}{\partial t^{2}} δ (s - l_{b} / 2)

(22)

Where $η$ and $K$ are the damping and coupling coefficient in the wake oscillator model. Also $ω_{f}$ is the vortex shedding frequency, which may be described in terms of Strouhal number, $S_{t},$ as:

ω_{f} = \frac{2 π S_{t} U_{f}}{D}

(23)

Now the differential equation governing on the harvested electric voltage is obtained. The electric displacements for the piezoelectric layer along direction 3, $D_{3},$ which is here correspond to direction $z$ can be written as below (Erturk and Inman, 2011)

D_{3} = d_{31} ε_{p} + ϵ_{33}^{T} E_{3}

(24)

Where $ϵ_{33}^{T}$ is the electrical permittivity factor (dielectric constant), and $d_{31}$ is the strain constant in the piezoelectric layer, also $E_{3}$ is electric field along the $z$ axis, which is expressed as

E_{3} = - \frac{V (t)}{h_{p}}

(25)

The strain is calculated in terms of the beam deflection as follow

ε_{p} = - \bar{z} \frac{\partial^{2} W (s, t)}{\partial s^{2}}

(26)

According to piezoelectric relations

d_{31} = \frac{e_{31}}{E_{p}}

(27)

Where $e_{31}$ is the piezoelectric coefficient. Based on the Gauss’ law, the total electric charge inside a closed volume is related to the integral of electric displacement on the surface of closed volume as below

Q = \int \int_{A} D_{3} . n dA

(28)

Where $Q$ is electric charge, $dA$ is the area element of the piezoelectric layer surface and $n$ is the unit normal vector of the surface. By derivation from equation (28) with respect to time variable and using Ohm’s law, it can be written as

\frac{\partial}{\partial t} \int \int D_{3} . n dA - \frac{V (t)}{R} = 0

(29)

Where $R$ is the electrical resistance. By substituting equations (25), (26) and (27) into equation (24) and then substituting the outcome into equation (29) one obtains

\int_{l_{1}}^{l_{1} + l_{p}} b_{p} e_{31} h_{p} \frac{\partial}{\partial t} (\frac{\partial^{2} W (s, t)}{\partial s^{2}}) ds + \frac{l_{p} b_{p} ϵ_{33}^{T}}{h_{p}} \frac{\partial V (t)}{\partial t} + \frac{V (t)}{R} = 0

(30)

Now, the differential equation governed on the motion equation, wake oscillator model and output electric voltage are written in dimensionless form by using the following variables changes

w = \frac{W}{D} & x = \frac{s}{l_{b}} & τ = \frac{t}{T_{b}} & T_{b} = \sqrt{\frac{ρ_{b} b_{b} h_{b} l_{b}^{4}}{E_{b} I_{1}^{b}}}

(31)

Thus, by substituting equation (31) into equation (19), the dimensionless form of motion equation will be as follow

\begin{matrix} \frac{\partial^{2}}{\partial x^{2}} (H (x) \frac{\partial^{2} w}{\partial x^{2}}) + m (x) \frac{\partial^{2} w}{\partial τ^{2}} + m_{0} δ (x - 1 / 2) \frac{\partial^{2} w}{\partial τ^{2}} \\ + c δ (x - 1 / 2) \frac{\partial w}{\partial τ} - (N + α_{1} \int_{0}^{l_{b}} \frac{\partial w}{\partial x} \frac{\partial w}{\partial x} dx - α_{3} V (t)) \frac{\partial^{2} w}{\partial x^{2}} \\ = - α_{4} V (t) (\frac{d^{2}}{d x^{2}} H_{\frac{l_{1}}{l_{b}}} - \frac{d^{2}}{d x^{2}} H_{\frac{l_{1} + l_{p}}{l_{b}}}) + F_{L} {Ω_{f}}^{2} q (τ) δ (x - 1 / 2) \end{matrix}

(32)

In which

m (x) = 1 + \frac{ρ_{p} h_{p} b_{p}}{ρ_{b} h_{b} b_{b}} (H_{\frac{l 1}{lb}} - H_{\frac{lp + l 1}{lb}}),

m_{0} = \frac{m_{T}}{ρ_{b} b_{b} h_{b}},

\begin{matrix} H (x) = (1 - H_{l 1 / lp}) + (H_{l 1 / lp} - H_{\frac{lp + l 1}{lb}}) \frac{E_{b} I_{2}^{b}}{E_{b} I_{1}^{b}} \\ + (H_{l 1 / lp} - H_{\frac{lp + l 1}{lb}}) \frac{E_{p} I_{2}^{p}}{E_{b} I_{1}^{b}} + (H_{\frac{lp + l 1}{lb}}), \end{matrix}

c = \hat{C} \frac{l_{b}^{4}}{E_{b} I_{1}^{b} T_{b}},

α_{1} = 6 (\frac{D}{h_{b}})^{2} \frac{1 + \frac{E_{p} h_{p}}{E_{b} h_{b}}}{(\frac{E_{p} h_{p}}{E_{b} h_{b}}) \times (1 - \frac{l_{p}}{l_{b}}) + 1},

N = N_{b} \frac{1 + \frac{E_{p} h_{p}}{E_{b} h_{b}}}{(\frac{E_{p} h_{p}}{E_{b} h_{b}}) (1 - \frac{l_{p}}{l_{b}}) + 1}, N_{b} = \frac{\bar{ε} A_{b} l_{b}}{I_{1}^{b}},

α_{3} = 6 α_{5} l_{p} \frac{\frac{E_{p}}{E_{b}}}{(\frac{E_{p} h_{p}}{E_{b} h_{b}}) \times (1 - \frac{l_{p}}{l_{b}}) + 1}, α_{5} = 2 \frac{l_{b}^{2} d_{31}}{h_{b}^{3}},

α_{4} = 3 α_{5} \frac{E_{p} \times (h_{p}^{2} + h_{b} h_{p} - 2 h_{p} {\tilde{z}}_{n})}{E_{b} D h_{p}},

F_{L} = \frac{l_{s} C_{L 0} ρ_{f} D^{2}}{16 π^{2} {S_{t}}^{2}} \frac{1}{ρ_{b} b_{b} h_{b}}

(33)

Also, by substituting equation (31) into equation (22), the dimensionless form of the equation governed wake oscillator will be as

\begin{matrix} \frac{\partial^{2} q (τ)}{\partial τ^{2}} + η Ω_{f} [q {(τ)}^{2} - 1] \frac{\partial q (τ)}{\partial τ} + {Ω_{f}}^{2} q (τ) \\ = K \frac{\partial^{2} w}{\partial τ^{2}} δ (x - 1 / 2), \\ Ω_{f} = T_{b} ω_{f} \end{matrix}

(34)

Where $Ω_{f}$ is the vortex shedding frequency in dimensionless form, and finally by substituting equation (31) into equation (30) the dimensionless form of the differential equation governed on the output electric voltage will be as

\int_{\frac{l_{1}}{l_{b}}}^{\frac{l_{1} + l_{p}}{l_{b}}} \frac{e_{31} b_{p} h_{p} D}{l_{b} T_{b}} \frac{\partial}{\partial τ} (\frac{\partial^{2} w (x, τ)}{\partial x^{2}}) dx + \frac{ϵ_{33} b_{p} l_{p}}{h_{p} T_{b}} \frac{\partial V (τ)}{\partial τ} + \frac{V (τ)}{R} = 0

(35)

3. Mode shape and natural frequency

Neglecting the nonlinear and damping terms in equation (32), the equation of motion will be written as follow

\begin{matrix} \frac{\partial^{2}}{\partial x^{2}} (H (x) \frac{\partial^{2} w (x, τ)}{\partial x^{2}}) + m (x) \frac{\partial^{2} w (x, τ)}{\partial τ^{2}} \\ + m_{0} δ (x - \frac{1}{2}) \frac{\partial^{2} w (x, τ)}{\partial τ^{2}} = 0 \end{matrix}

(36)

Now, the free vibration response of the system is assumed as

w (x, τ) = Φ_{n} (x) e^{i ω τ}

(37)

Where $ω$ is the natural frequency of the system, and $Φ_{n} (x)$ is the mode shape of the system correspond to natural frequency $ω$ . By substituting equation (37) into equation (36), the differential equation governing the mode shape of the system is written as follow

\begin{matrix} \frac{d^{2}}{d x^{2}} (H (x) \frac{d^{2} Φ_{n} (x)}{d x^{2}}) - ω^{2} m (x) Φ_{n} (x) \\ - ω^{2} m_{T} δ (x - \frac{1}{2}) Φ_{n} (x) = 0 \end{matrix}

(38)

Equation (38) is solved by applying Galerkin method where the mode shape of a clamped-clamped beam without piezoelectric layer, $φ_{i} (x),$ are used as comparison functions (Nayfeh and Pai, 2004). Therefore, it is assumed that

Φ_{n} (x) = \sum_{i = 1}^{n} C_{i} φ_{i} (x)

(39)

Where $C_{i}, i = 1 . . n$ , are the constant coefficients that must be obtained after applying Galerkin approach. So, substituting equation (39) into (38), and multiplying the outcome by $φ_{j} (x)$ , and integrating it at the interval $x = 0 . . 1$ , it results that

\begin{matrix} \int_{0}^{1} \sum_{i = 1}^{n} C_{i} \frac{d^{2}}{d x^{2}} (H (x) \frac{d^{2} φ_{i} (x)}{d x^{2}}) φ_{j} (x) dx - ω^{2} \\ \int_{0}^{1} \sum_{i = 1}^{n} C_{i} φ_{i} (x) m (x) φ_{j} (x) dx - ω^{2} \\ \int_{0}^{1} \sum_{i = 1}^{n} C_{i} φ_{i} (x) m_{T} δ (x - \frac{1}{2}) φ_{j} (x) dx = 0 \end{matrix}

(40)

The natural frequency of system is obtained setting the determinants of constant coefficients $C_{i}$ equal to zero in equation (40). After that substituting the natural frequency in the algebraic system, the coefficients of $C_{i}$ are obtained.

4. Discretization of motion equation

In previous section the natural mode shapes of system are obtained. Now these mode shapes are used as comparison functions to discretize the partial differential equation governed on the motion equation, wake oscillator model and output electric voltage, therefore it is assumed that:

w (x, τ) = Φ_{n} (x) P (τ)

(41)

Where $P (τ)$ is time coordinate of the response. By substituting equation (41) into equation (32) and multiplying the outcome by $Φ_{n} (x)$ , and integrating it at the interval $x = 0 . . 1$ , one obtains

\begin{array}{l} \frac{\partial^{2} P (τ)}{\partial τ^{2}} + ω_{n}^{2} P (τ) + C_{e q} \frac{\partial P (τ)}{\partial τ} - K_{c u b i c} P {(τ)}^{3} \\ + α_{3 e q} V (t) P (τ) - N_{e q} P (τ) = α_{4 e q} V (t) + F_{e q} Ω f^{2} q (τ) \end{array}

(42)

Where the terms including coefficients of $K_{cubic}$ , $N_{eq}$ , $C_{eq}$ and $F_{eq}$ are presented due to considering the middle-layer stretching, pre axial load, viscosity damping and lift force effects, respectively. Moreover, the terms including coefficient $α_{3 eq}$ and $α_{4 eq}$ are presented due to the piezoelectric effect. These coefficients are described as below:

K_{cubic} = \frac{α_{1} \int_{0}^{1} Φ_{n} (\int_{0}^{1} {(\frac{d Φ_{n}}{dx})}^{2} dx) \frac{d^{2} Φ_{n} (x)}{d x^{2}} dx}{m_{eq}},

m_{eq} = \int_{0}^{1} m (x) {Φ_{n}}^{2} dx + m_{0} {Φ_{n}}^{2} |_{\begin{matrix} x = 0.5 \end{matrix}},

C_{eq} = \frac{c {Φ_{n}}^{2} |_{\begin{matrix} \begin{matrix} x = 0.5 \end{matrix} \end{matrix}}}{m_{eq}},

α_{3 eq} = \frac{α_{3} \int_{\frac{l_{1}}{l_{b}}}^{\frac{l_{1} + l_{p}}{l_{b}}} Φ_{n} \frac{d^{2} Φ_{n}}{d x^{2}} dx}{m_{eq}},

N_{eq} = \frac{N \int_{0}^{1} Φ_{n} \frac{d^{2} Φ_{n}}{d x^{2}} dx}{m_{eq}},

\begin{matrix} α_{4 eq} = \frac{α_{4} (\frac{d Φ_{n}}{dx} |_{\begin{matrix} x = \frac{l_{1}}{l_{b}} \end{matrix}} - \frac{d Φ_{n}}{dx} |_{\begin{matrix} x = \frac{l_{1} + l_{p}}{l_{b}} \end{matrix}})}{m_{eq}}, \end{matrix}

F_{eq} = \frac{F_{L} {Φ_{n}}^{2} |_{\begin{matrix} \begin{matrix} x = 0.5 \end{matrix} \end{matrix}}}{m_{eq}}

(43)

By substituting equation (41) into equation (34)

\frac{d^{2} q (τ)}{d τ^{2}} + η Ω_{f} [q {(τ)}^{2} - 1] \frac{dq (τ)}{d τ} + {Ω_{f}}^{2} q (τ) = G \frac{d^{2} P (τ)}{d τ^{2}}

(44)

In which

G = K Φ_{n} |_{\begin{matrix} x = 0.5 \end{matrix}}

(45)

Also by substituting equation (41) into equation (35)

χ \frac{dP (τ)}{d τ} + C_{p} \frac{dV (τ)}{d τ} + \frac{V (τ)}{R} = 0

(46)

In which

\begin{matrix} χ = \frac{e_{31} b_{p} h_{p} D}{l_{b} T_{b}} (\frac{d Φ_{n}}{dx} |_{\begin{matrix} x = \frac{l_{1 + l_{p}}}{l_{b}} \end{matrix}} - \frac{d Φ_{n}}{dx} |_{\begin{matrix} x = \frac{l_{1}}{l_{b}} \end{matrix}}), \\ C_{p} = \frac{ϵ_{33} b_{p}}{h_{p} T_{b}} l_{p} \end{matrix}

(47)

As seen, the discretized equation governing on the motion, wake oscillator and the output electric voltage are coupled to gather.

5. Stability and modal analysis

Neglecting the nonlinear and damping terms, the motion equations governing on the oscillations and wake oscillator model will be written as follow

\begin{matrix} \frac{d^{2} P (τ)}{d τ^{2}} + ω_{n}^{2} P (τ) = F_{eq} {Ω_{f}}^{2} q (τ), \\ \frac{d^{2} q (τ)}{d τ^{2}} + {Ω_{f}}^{2} q (τ) = G \frac{d^{2} P (t)}{d τ^{2}} \end{matrix},

(48)

The response of the equation (48) may be considered as follows

P (τ) = \bar{p} e^{υ τ} & q (τ) = \bar{q} e^{υ τ}

(49)

By substituting equation (49) into equation (48)

\begin{matrix} (ν^{2} + ω_{n}^{2}) \bar{p} - F_{eq} {Ω_{f}}^{2} \bar{q} = 0, \\ - G ν^{2} \bar{p} + (ν^{2} + {Ω_{f}}^{2}) \bar{q} = 0 \end{matrix}

(50)

The nonzero solution is obtained for system when the determinant of the coefficient $\bar{p}$ and $\bar{q}$ be equal to zero. It results that

ν^{4} + (ω_{n}^{2} + {Ω_{f}}^{2} - F_{eq} G {Ω_{f}}^{2}) ν^{2} + ω_{n}^{2} {Ω_{f}}^{2} = 0

(51)

By solving the algebra equation (51)

\begin{matrix} 2 (ν_{\pm})^{2} = - (ω_{n}^{2} + {Ω_{f}}^{2} - F_{eq} G {Ω_{f}}^{2}) \\ \pm \sqrt{{(ω_{n}^{2} + {Ω_{f}}^{2} - F_{eq} G {Ω_{f}}^{2})}^{2} - 4 ω_{n}^{2} {Ω_{f}}^{2}} \end{matrix}

(52)

When the value of expression inside the radical is positive then the eigenvalues will be pure imaginary, and response for linear system will be periodic, and when the value of expression inside the radical is zero, two eigenvalues will be matched to each other. Therefore, by setting the expression under the radical equal to zero, one obtains:

{(ω_{n}^{2} + (1 - F_{eq} G) {Ω_{f}}^{2})}^{2} - 4 ω_{n}^{2} {Ω_{f}}^{2} = 0

(53)

By solving equation (53), two values of $Ω_{f}$ that causes the eigenvalues are matched will be as

{Ω_{f}}_{\pm}^{2} = \frac{ω_{n}}{1 \pm \sqrt{F_{eq} G}}

(54)

In the range of ${Ω_{f}}_{+} < Ω_{f} < {Ω_{f}}_{-}$ , the undamped linear system will be unstable. By substituting vortex shedding frequency in dimensionless form into equation (54), the wind speed boundary value for unstable solution will be as

U = \frac{D}{2 π S_{T} T_{b}} \frac{ω_{n}}{1 \pm \sqrt{F_{eq} G}}

(55)

According to the equation (50), the mode shape of linear system vibration considering the wake oscillator effect will be as:

\frac{\bar{p}}{\bar{q}} = \frac{F_{eq} {Ω_{f}}^{2}}{(ν^{2} + ω_{n}^{2})} = \frac{(ν^{2} + {Ω_{f}}^{2})}{G ν^{2}}

(56)

6. Perturbation analysis

In previous analysis the effect of nonlinear term has been neglected. It causes that system be stable or unstable depend on the flow velocity. In practice, system will experience the periodic oscillation due to the nonlinear effect which is not predicted by linear solution. Therefore, in this case, the perturbation technique is used for predicting the nonlinear solution. Firstly, the booking parameter $ε = 1$ , which represents the order of terms is used as follow

\begin{matrix} \frac{\partial^{2} P (τ)}{\partial τ^{2}} + ω_{n}^{2} P (τ) + ε C_{eq} \frac{\partial P (τ)}{\partial τ} - ε K_{cubic} P {(τ)}^{3} \\ + ε α_{3 eq} V (t) P (τ) = ε α_{4 eq} V (t) + ε F_{eq} {Ω_{f}}^{2} q (τ) \end{matrix}

\frac{\partial^{2} q (τ)}{\partial τ^{2}} + ε η Ω_{f} [q {(τ)}^{2} - 1] \frac{\partial q (τ)}{\partial τ} + {Ω_{f}}^{2} q (τ) = ε G \frac{\partial^{2} P (τ)}{\partial τ^{2}},

χ \frac{\partial P (τ)}{\partial τ} + C_{p} \frac{dV (τ)}{d (τ)} + \frac{V (τ)}{R} = 0

(57)

Now the response of equation (57) is assumed as follows

P (τ) = P_{0} (T_{0}, T_{1}) + ε P_{1} (T_{0}, T_{1}),

q (τ) = q_{0} (T_{0}, T_{1}) + ε q_{1} (T_{0}, T_{1}),

V (τ) = V_{0} (T_{0}, T_{1}) + ε V_{1} (T_{0}, T_{1})

(58)

Where $T_{0} = τ$ and $T_{1} = ε τ$ are time scales. By substituting equation (58) into equation (57) and separating terms with equal order of $ε$ , it results that

Order (1)

\frac{\partial^{2}}{\partial {T_{0}}^{2}} P_{0} (T_{0}, T_{1}) + {ω_{n}}^{2} P_{0} (T_{0}, T_{1}) = 0,

\frac{\partial^{2}}{\partial {T_{0}}^{2}} q_{0} (T_{0}, T_{1}) + {Ω_{f}}^{2} q_{0} (T_{0}, T_{1}) = 0,

C_{p} \frac{\partial}{\partial T_{0}} V_{0} (T_{0}, T_{1}) + \frac{1}{R} V_{0} (T_{0}, T_{1}) = - χ \frac{\partial}{\partial T_{0}} P_{0} (T_{0}, T_{1})

(59)

Order ( $ε$ )

\begin{matrix} \frac{\partial^{2}}{\partial {T_{0}}^{2}} P_{1} (T_{0}, T_{1}) + {ω_{n}}^{2} P_{1} (T_{0}, T_{1}) = \\ - 2 \frac{\partial^{2}}{\partial T_{0} \partial T_{1}} P_{0} (T_{0}, T_{1}) - C_{eq} \frac{\partial}{\partial T_{0}} P_{0} (T_{0}, T_{1}) \\ + K_{cubic} P_{0} {(T_{0}, T_{1})}^{3} - α_{3 eq} (P_{0} (T_{0}, T_{1}) . V_{0} (T_{0}, T_{1})) \\ + α_{4 eq} V_{0} (T_{0}, T_{1}) + F_{eq} {Ω_{f}}^{2} q_{0} (T_{0}, T_{1}), \end{matrix}

\begin{matrix} \frac{\partial^{2}}{\partial {T_{0}}^{2}} q_{1} (T_{0}, T_{1}) + {Ω_{f}}^{2} q_{1} (T_{0}, T_{1}) = \\ - 2 \frac{\partial^{2}}{\partial T_{0} \partial T_{1}} q_{0} (T_{0}, T_{1}) - η Ω_{f} [q_{0} {(T_{0}, T_{1})}^{2} - 1] \frac{\partial}{\partial T_{0}} q_{0} (T_{0}, T_{1}) \\ + G \frac{\partial^{2}}{\partial {T_{0}}^{2}} P_{0} (T_{0}, T_{1}), \end{matrix}

\begin{matrix} C_{p} \frac{\partial}{\partial T_{0}} V_{1} (T_{0}, T_{1}) + \frac{1}{R} V_{1} (T_{0}, T_{1}) = - C_{p} \frac{\partial}{\partial T_{1}} V_{0} (T_{0}, T_{1}) \\ - χ (\frac{\partial}{\partial T_{1}} P_{0} (T_{0}, T_{1}) + \frac{\partial}{\partial T_{0}} P_{1} (T_{0}, T_{1})) \end{matrix}

(60)

The solution of equation (59) will be obtained as follow

P_{0} = A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}},

q_{0} = B (T_{1}) e^{i Ω_{f} T_{0}} + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}},

V_{0} = \frac{- i ω_{n} χ \bar{A} (T_{1}) R e^{- i ω_{n} T_{0}}}{iR C_{p} ω_{n} - 1} - \frac{i ω_{n} χ A (T_{1}) R e^{i ω_{n} T_{0}}}{iR C_{p} ω_{n} + 1} + \bar{E} e^{- \frac{T_{0}}{R C_{p}}}

(61)

Where $A (T_{1})$ and $B (T_{1})$ are the complex coefficient and $\bar{A} (T_{1})$ and $\bar{B} (T_{1})$ are their complex conjugates. It is clear that the term of $\bar{E} e^{- \frac{T_{0}}{R C_{p}}}$ will be equal to zero in steady state response. It means that obtaining the coefficient of $A (T_{1})$ and $B (T_{1})$ is sufficient for steady state response of system. Therefore, by substituting equation (61) into equation (60)

\begin{matrix} \frac{\partial^{2}}{\partial {T_{0}}^{2}} P_{1} (T_{0}, T_{1}) + {ω_{n}}^{2} P_{1} (T_{0}, T_{1}) = - 2 \frac{\partial^{2}}{\partial T_{0} \partial T_{1}} \\ (A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}}) \\ - C_{eq} \frac{\partial}{\partial T_{0}} (A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}}) \\ + K_{cubic} {(A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}})}^{3} \\ - α_{3 eq} (A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}}) \\ \times (\frac{- i ω_{n} χ \bar{A} R e^{- i ω_{n} T_{0}}}{iR C_{p} ω_{n} - 1} - \frac{i ω_{n} χ \bar{A} R e^{i ω_{n} T_{0}}}{iR C_{p} ω_{n} + 1} + \bar{E} e^{- \frac{T_{0}}{R C_{p}}}) \\ + α_{4 eq} (\frac{- i ω_{n} χ \bar{A} R e^{- i ω_{n} T_{0}}}{iR C_{p} ω_{n} - 1} - \frac{i ω_{n} χ \bar{A} R e^{i ω_{n} T_{0}}}{iR C_{p} ω_{n} + 1} + \bar{E} e^{- \frac{T_{0}}{R C_{p}}}) \\ + F_{eq} {Ω_{f}}^{2} (B (T_{1}) e^{i Ω_{f} T_{0}} + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}}) \end{matrix}

\begin{matrix} \frac{\partial^{2}}{\partial {T_{0}}^{2}} q_{1} (T_{0}, T_{1}) + {Ω_{f}}^{2} q_{1} (T_{0}, T_{1}) = \\ - 2 i Ω_{f} (\frac{\partial}{\partial T_{1}} B (T_{1}) e^{i Ω_{f} T_{0}} - \frac{\partial}{\partial T_{1}} \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}}) \\ - η i {Ω_{f}}^{2} [{(B (T_{1}) e^{i Ω_{f} T_{0}} + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}})}^{2} - 1] (B (T_{1}) e^{i Ω_{f} T_{0}} \\ + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}}) + G \frac{\partial^{2}}{\partial {T_{0}}^{2}} (A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}}) \end{matrix}

(62)

The unknown’s coefficients $A (T_{1})$ and $B (T_{1})$ are obtained by setting the secular terms equal to zero in the above equation. Now considering $Ω_{f}$ as follows

Ω_{f} = ω_{n} + σ ε

(63)

Where $σ$ is the detuning parameter, by substituting equation (63) into equation (62)

\begin{array}{l} \frac{\partial^{2}}{\partial {T_{0}}^{2}} P_{1} (T_{0}, T_{1}) + {ω_{n}}^{2} P_{1} (T_{0}, T_{1}) = \\ [- 2 i ω_{n} \frac{\partial}{\partial T_{1}} A (T_{1}) - C_{e q} i ω_{n} A (T_{1}) + K_{c u b i c} 3 A {(T_{1})}^{2} \bar{A} (T_{1}) \\ - α_{4 e q} (\frac{i ω_{n} χ A R}{i R C_{p} ω_{n} + 1}) - α_{3 e q} A (T_{1}) \bar{E} e^{- \frac{T_{0}}{R C_{p}}} \\ + F_{e q} {Ω_{f}}^{2} B (T_{1}) e^{i (σ ε) T_{0}}] e^{i ω_{n} T_{0}} \\ - α_{4 e q} \bar{E} e^{- \frac{T_{0}}{R C_{p}}} + α_{3 e q} \frac{i ω_{n} χ A \bar{A} R}{i R C_{p} ω_{n} + 1} + C C + N S T, \end{array}

\begin{array}{l} \frac{\partial^{2}}{\partial {T_{0}}^{2}} q_{1} (T_{0}, T_{1}) + {Ω_{f}}^{2} q_{1} (T_{0}, T_{1}) \\ = [- 2 i Ω_{f} \frac{\partial}{\partial T_{1}} B (T_{1}) - η i {Ω_{f}}^{2} B {(T_{1})}^{2} \bar{B} (T_{1}) + η i {Ω_{f}}^{2} B (T_{1}) \\ - G {ω_{n}}^{2} A (T_{1}) e^{- i σ ε T_{0}}] e^{i Ω_{f} T_{0}} \\ + C c + N S T \end{array}

(64)

Where $Cc$ indicates conjugate of complex term, and $NST$ demonstrates all terms which are not secular. To obtain the solvability condition the coefficients of $e^{i ω_{n} T_{0}}$ and $e^{i Ω_{f} T_{0}}$ must be equal to zero, therefore

\begin{matrix} - 2 i ω_{n} \frac{\partial}{\partial T_{1}} A (T_{1}) - C_{eq} i ω_{n} A (T_{1}) + K_{cubic} 3 A {(T_{1})}^{2} \bar{A} (T_{1}) \\ - α_{4 eq} (\frac{i ω_{n} χ AR}{iR C_{p} ω_{n} + 1}) \\ - α_{3 eq} A (T_{1}) \bar{E} e^{- \frac{T_{0}}{R C_{p}}} + F_{eq} {Ω_{f}}^{2} B (T_{1}) e^{i (σ ε) T_{0}} = 0, \end{matrix}

\begin{matrix} - 2 i Ω_{f} \frac{\partial}{\partial T_{1}} B (T_{1}) - η i {Ω_{f}}^{2} B {(T_{1})}^{2} \bar{B} (T_{1}) + η i {Ω_{f}}^{2} B (T_{1}) \\ - G {ω_{n}}^{2} A (T_{1}) e^{- i σ ε T_{0}} = 0 \end{matrix}

(65)

Now $A (T_{1})$ and $B (T_{1})$ are considered as polar form:

A (T_{1}) = \frac{1}{2} a (T_{1}) e^{i α (T_{1})} & B (T_{1}) = \frac{1}{2} b (T_{1}) e^{i β (T_{1})}

(66)

Noting that $ε T_{0} = T_{1}$ , and substituting equation (66) into equation (65) and multiplying $e^{- i α}$ at the first equation and multiply $e^{- i β}$ at the second equation

\begin{matrix} - i ω_{n} \overset{\cdot}{a} + ω_{n} a \overset{\cdot}{α} - \frac{C_{eq} i ω_{n} a}{2} + K_{cubic} \frac{3}{8} a^{3} \\ - \frac{1}{2} α_{4 eq} \frac{i ω_{n} χ aR}{iR C_{p} ω_{n} + 1} - \frac{1}{2} α_{3 eq} a \bar{E} e^{- \frac{T_{0}}{R C_{p}}} \\ + \frac{1}{2} F_{eq} {Ω_{f}}^{2} b e^{i (β - α - σ T_{1})} = 0, \end{matrix}

- i Ω_{f} \overset{\cdot}{b} + Ω_{f} b \overset{\cdot}{β} - \frac{i {Ω_{f}}^{2} η b^{3}}{8} + \frac{i {Ω_{f}}^{2} η b}{2} - \frac{G {ω_{n}}^{2} a}{2} e^{i (α - β - σ T_{1})} = 0

(67)

By separating the real part and imaginary part of equation (67) and equating the results equal to zero

\begin{matrix} \overset{\cdot}{α} = - \frac{3 α_{1 eq} a^{2}}{8 ω_{n}} - \frac{α_{3 eq} \bar{E} e^{- \frac{T_{0}}{R C_{p}}}}{2 ω_{n}} + \frac{α_{4 eq} ω_{n} R^{2} χ C_{p}}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} \\ - \frac{F_{eq} {Ω_{f}}^{2} b}{2 a ω_{n}} Cos (β - α - σ T_{1}), \end{matrix}

\overset{\cdot}{β} = \frac{G {ω_{n}}^{2} a}{2 b Ω_{f}} Cos (α - β - σ T_{1}),

\overset{\cdot}{a} = - \frac{C_{eq, n} a}{2} - \frac{α_{4 eq} a χ R}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} + \frac{F_{eq} {Ω_{f}}^{2} b}{2 ω_{n}} Sin (β - α - σ T_{1}),

\dot{b} = - \frac{η Ω_{f} b}{8} (b^{2} - 4) - \frac{G {ω_{n}}^{2} a}{2 Ω_{f}} S i n (α - β - σ T_{1})

(68)

Now by considering $ψ = β - α + σ T_{1}$ , it results that

\begin{array}{l} \dot{ψ} = - \frac{3 K_{c u b i c} a^{2}}{8 ω_{n}} + \frac{α_{3 e q} \bar{E} e^{- \frac{T_{0}}{R C_{p}}}}{2 ω_{n}} - \frac{α_{4 e q} ω_{n} R^{2} χ C_{p}}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} \\ + (\frac{G {ω_{s}}^{2} a}{2 b Ω_{f}} + \frac{F_{e q} {Ω_{f}}^{2} b}{2 a ω_{n}}) C o s (ψ) + σ, \end{array}

\overset{\cdot}{a} = - \frac{C_{eq} a}{2} - \frac{α_{4 eq} a χ R}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} + \frac{F_{eq} {Ω_{f}}^{2} b}{2 ω_{n}} Sin (β - α - σ T_{1}),

\overset{\cdot}{b} = - η Ω_{f} b (\frac{1}{8} b^{2} - \frac{1}{2}) + \frac{G {ω_{s}}^{2} a}{2 b Ω_{f}} Sin (ψ)

(69)

The equation governing equilibrium solution amplitude will be obtained by placing derivative of $a$ , $β$ and $ψ$ equal to zero. By substituting equation (66) into equation (61)

P_{0} = \frac{1}{2} a (T_{1}) e^{i α (T_{1})} e^{i ω_{n} T_{0}} + \frac{1}{2} a (T_{1}) e^{- i α (T_{1})} e^{- i ω_{n} T_{0}},

q_{0} = \frac{1}{2} b (T_{1}) e^{i β (T_{1})} e^{i Ω_{f} T_{0}} + \frac{1}{2} b (T_{1}) e^{- i β (T_{1})} e^{- i Ω_{f} T_{0}}

(70)

According to the obtained results, when the effect of damping and nonlinearity are considered, the frequency presented in wake oscillation is $Ω_{f} + \overset{\cdot}{β}$ and the frequency presented in structure oscillation is $ω_{n} + \overset{\cdot}{α}$ .

7. Numerical solutions

Here the discretized equation obtained through Galerkin approach are solved using Runge–Kutta–Fehlberg numerical method, which is called RKF45. Calculations are done using Maple software. Firstly, the initial conditions of the system are considered equal to zero, and frequency of shedding vortex is tuned much less than the natural frequcy of linear system. Then applying the numerical algorithm, time history of the structure oscillations, wake oscillations and electric output voltage are obtained. After that, the amplitude of the steady-state is recorded. The value of steady state solution is considered as initial condition for the next step. Moreover the vortex shedding frequency is increased by an increase of the value of flow speed. The equations are again solved and the amplitudes in steady state response are recorded. This process is followed until the Lock-in area is observed and terminated.

8. Results and discussion

The response of the system is obtained according to the characteristics stated in Table 1.

Table 1.

Geometric and material properties of the beam and piezoelectric layer.

$h_{b} (μ m)$	$l_{b} (mm)$	$E_{b} (Gpa)$	$b_{b} (mm)$	$ρ_{b} (\frac{Kg}{m^{3}})$	$l_{s} (mm)$	$D (mm)$
71	130	200	7.3	7800	400	20
$h_{p} (μ m)$	$l_{p} (mm)$	$E_{p} (Gpa)$	$b_{p} (mm)$	$ρ_{p} (\frac{Kg}{m^{3}})$	$d_{31} (mm) (\frac{pC}{N})$	$ϵ_{33} (\frac{nF}{m})$
300	60	30.336	7.3	5440	−2.1	3.15

The first, second and third free vibrations mode shapes of the system and their corresponding natural frequencies have been shown in Figure 5 and Table 2, respectively.

Table 2.

First, second and third natural frequency of the structure.

$ω_{1}$	$ω_{2}$	$ω_{3}$
$36.03587$	$173.6000$	$304.09128$

Variations of the eigenvalues of the system with respect to the variations of flow speed considering the effect of wake oscillator are shown in Figure 6 using equation (52). The first eigenvalue belongs to the structure, and for low values of flow speed its amount is approximately constant, and equals to the natural frequency of the free vibration of the system. The second eigenvalue belongs to the wake oscillator, and increases by increase of the value of flow speed until it reaches approximately to first eigenvalue. After this amount, the eigenvalues will be imaginary with positive real part, and the unstable domain will start. This domain continues by increasing the value of flow speed until two pure imaginary eigenvalues be observed again.

The stability analysis shown in Figure 6 demonstrates that if one neglects the nonlinear terms, then the amplitude of vibration increases as continuous by increase of time which does not occur in practice. Therefore, in follow the variations of amplitude of the system oscillations and harvested voltage is evaluated considering the nonlinear effect.

Equation (42) states that the effect of middle layer stretching appears as a cubic stiffness term in discretized equation of motion. Therefore, setting this coefficient to zero, this effect may be removed. The variations of the amplitude of the system oscillation and harvested voltage with respect to the variations of flow speed with and without considering the middle layer stretching effect for different values of the cylinder diameter is shown in Figure 7. It denotes that when the diameter of the attached cylinder is low, the dynamic behavior of systems with and without considering the middle layer stretching effect is approximately similar. It states that by an increase of the value of cylinder diameter, the difference between the dynamic behavior of system with and without the middle layer stretching effect will be considerable. The dynamic behavior of system without middle layer stretching effect shows that in a region of flow speed which is called Lock-in region, the amplitude of the oscillations first increases with a large slope from zero to a considerable value, then it remains almost constant, and then it decreases with a large slope to zero. Figure 7 shows that when the diameter of the cylinder is large and middle layer stretching is considered, then the slope of increasing amplitude at beginning of the lock-in domain is lower than the state without considering the stretching effect, but the lengths of lock in domain is larger. It shows that the starting of the lock in domain are the same for both systems with and without stretching effect. It is expected because the amplitude at the start of the lock in is approximately equal to zero and so the stretching effect has not been appeared. After starting the lock-in domain the dynamic behavior is due to the competition between lift force, linear stiffness and nonlinear stiffness from stretching. According to Table 3, as the diameter increases, the dimensionless values of lift force coefficient increases and the linear spring coefficient decreases, thus the value of oscillation amplitude increases. As the amplitude increases, the effect of the middle layer stretching on the behavior of the system will be significant. For small diameters the amplitude is small since the lift force is small, therefore the nonlinear stretching effect will be also small. Therefore, a considerable difference is not appeared for system with and without stretching effect. For system with large diameters of the cylinder, the amplitude will be large due to the effect of big lift force. So, if one considers the stretching effect, its effect will be considerable.

Table 3.

Dimensionless coefficients of lift force, linear stiffness term and non-linear stiffness term obtained from numerical solution.

Cylinder diameter	Lift force	Non-linear spring	Linear spring
10	0.00021	4,920,000	1670
12	0.00029	6,790,000	1600
14	0.00038	8,810,000	1530
16	0.00047	10,900,000	1450
18	0.00057	13,100,000	1370
20	0.00066	15,200,000	1300

As mentioned before, at the beginning of the lock in domain, there is no difference between the dynamic behavior of system with and without the stretching effect, but by increase of the value of wind speed, and increasing the oscillation amplitude, the stiffness due to the middle layer stretching effect will be added to the system. Considering equation (53) the final speed for closing the instability domain which corresponds approximately to finalizing the lock in domain depends on the natural frequency of system. It is clear that the nonlinear cubic stiffness due to stretching effect $, K_{cubic} P {(τ)}^{3},$ is similar to a linear stiffness as $K_{EQ} P (τ)$ where $K_{EQ} = k_{cubic} P {(τ)}^{2}$ . It means that the natural frequency of the equivalent system increase by increase of the value of $P (τ) .$ therefore the final speed for enclosing the unstable domain as well as the final speed of enclosing the lock in domain increases.

In order to better study the effect of the cylinder diameter on the behavior of the system, the diagram of amplitude variations with and without considering the middle layer stretching effect is shown separately in Figure 8(a) and (b). Figure 8(a) states that in the case without considering the middle layer stretching effect, the required speed to reach the Lock-in area increases by increasing the diameter of the cylinder.

Figure 8(b) shows that the same behavior is also observed in the system with middle layer stretching effect. It is due to the fact that before starting the lock in phenomena, the vortex shedding frequency is $ω_{f} = 2 π S_{t} U / D$ . In other words, the vortex shedding frequency decreases by increasing the diameter of the cylinder. The natural frequency of the structure is also decreased due to the decrease of the equivalent linear spring shown in Table 3. So, the competition between these parameters determine the flow speed for staring the lock in phenomena. Since the vortex shedding frequency is in inverse relation to the cylinder diameter, the rate of decreasing of the vortex shedding frequency will be much lower than decreasing the natural frequency. Considering that at the beginning of the Lock-in area, the vortex shedding frequency must be almost equal to the linear natural frequency, therefore by an increase of the cylinder diameter, the lock-in will be started at the larger wind speed. This condition is maintained for the case with considering middle layer stretching effect since the stretching effect is not appeared at the start of lock-in. In the system without middle layer stretching effect, the maximum amplitude for different values of the cylinder diameter is also increased by an increase of the diameter. As shown in Table 3, the coefficient of lift force is also increased by an increase of the cylinder diameter. Noting that the coefficient of linear stiffness is also decreased, therefore the maximum amplitude increases. The important and significant point is that the same behavior is also observed in the system with middle layer stretching effect but the amount of increasing of the maximum amplitude is lower due to affecting the middle layer stretching effect. Considering Figure 8(b), although for larger diameters, the maximum amplitude seems equal, but in actual the oscillation amplitude has been increased since the vertical axis show a dimensionless value W/D. It means that if the value be the same, the displacement W is system with larger diameters is larger.

The variations of the harvested voltages with considering and without considering the stretching effect has been shown in Figure 9. As seen before, the maximum amplitude of the oscillations increases by increase of the value of cylinder diameter. By increasing the amplitude of the oscillation, the strain of the piezoelectric layer increases and following that the output voltage increases. In addition, Figure 9 demonstrates that when the middle layer stretching effect is appeared then the energy may be collected in a wider range of wind speeds compare to the case that this effect is not appeared.

The natural frequency of the system decreases by an increase of the value of cylinder length because it causes that the mass of system increases. In addition, according to equation (20), the coefficient of lift force increases by increasing the length of the cylinder but this parameter has no effect on the vortex shedding frequency before starting the lock in phenomena. In other words, as shown in Figure 10(a) the flow speed for beginning of the Lock-in decreases by an increase of the value of cylinder length, and since the amplitude of the lift force increases directly therefore the oscillation amplitude is also increased compared to a system with smaller cylinder length.

Comparing Figure 10(a) with Figure 10(b) shows that the variations pattern of output voltage with respect to the variations of flow speed for different values of cylinder length is similar to the variation pattern of structure oscillation amplitude. It is due to the fact that output voltage is in direct relation to the strain of piezoelectric layer which is in direct relation to the structure oscillation amplitude.

Now, the variation effect of length of the beam is considered on the amplitude oscillation and harvested voltage. It is clear that the stiffness of the system decreases by an increase of the value of beam length. It causes that the natural frequency of the system decreases. It must be noted that the frequency of vortex shedding before lock-in phenomena is approximately independent of natural frequency of the structure. It means that system reaches to the synchronization region at lower speeds by increasing of beam lengths, see Figure 11(a). As mentioned before the output voltage depends on the strain on the piezoelectric layer which is in direct relation to the structure oscillation. It means that as shown in Figure 11(b), the variations pattern of output voltage with respect to the variations of beam length must follow the variation pattern of the structure oscillation. We know that the effect of decreasing the beam thickness on the natural frequency of system is similar to the effect of increasing the beam length. In addition, beam thickness variations have no direct effect on the frequency of vortex shedding before starting the lock-in phenomena. Therefore, a similar behavior to Figure 11 is expected if the beam thickness is decreased which is shown in Figure 12.

Here, the results obtained by perturbation theory is compared to the results of numerical method for both cases with and without considering the middle layer stretching effect. According to equation (69), the effect of middle layer stretching is observed as a term with coefficient of $K_{cubic}$ , so by setting this term equal to zero the behavior of the system in the case without considering middle layer stretching effect is obtained. As shown in Figure 13, the solid line represents the solution of the system by the numerical method and the dashed line represents the solution by multiple scale perturbation method. It shows that the dynamic behavior of the system in both numerical and analytical methods are completely similar to each other. Therefore, it may be said that the results of numerical methods shown in previous figures are confirmed.

If the effect of middle layer stretching effect is taken into account, and this nonlinear term is weak a good agreement is expected between the results obtained by numerical and analytical method. It has been shown in Figure 14 that the numerical method predicts increasing the maximum amplitude by an increase of the cylinder diameter. It means that the difference between the response obtained by perturbation method and numerical method must be increased by an increase of the value of cylinder diameter because more displacements more nonlinear stretching. This prediction is verified in Figure 14.

9. Conclusions

In this paper the effect of middle layer stretching effect on the dynamic behavior of an energy harvester clamped-clamped piezoelectric beam has been investigated. Two foam cylinders attached together as a dumbbell has been assumed mounted on the middle part of the beam for vortex-induced vibrations. The coupled differential equations governing on the structure oscillation, harvested voltage and fluid lift force has been established applying the Hamilton’s principle, Gauss law and wake oscillator model. Then the obtained equations has been discretized using Galerkin method. The mode shape of free vibration of system has been used as comparison function in discretizing process. A stability analysis on the linear system has been done, and a domain of instability which is approximately corresponds to the lock-in phenomena has been obtained. The nonlinear discretized equations governing the dynamic of system have been solved by both numerical and analytical perturbation method. The effect of geometrical parameters of system such as cylinder length, beam length and beam thickness has been investigated on the amount of lock-in domain and harvested voltage with and without considering the stretching effect. It has been shown that for large values of cylinder diameter the difference between numerical and theoretical results increases due to increasing the middle layer stretching effect. It was shown that the middle layer stretching effect has a significant effect on the lock-in phenomena, and the output electric voltage must be evaluated considering this effect. It has been shown that when the diameter of the cylinder is large and middle layer stretching is considered, then the slope of increasing amplitude at beginning of the lock-in domain is lower than the state without considering the stretching effect, but the lengths of lock in domain is larger. It has been shown that the starting of the lock in domain are the same for both systems with and without stretching effect.

Footnotes

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Mobin Alimanesh

Mehdi Zamanian

References

Abdelkefi

Hajj

Nayfeh

(2012) Power harvesting from transverse galloping of square cylinder. Nonlinear Dynamics 70: 1355–1363.

Akaydin

Elvin

Andreopoulos

(2010) Energy harvesting from highly unsteady fluid flows using piezoelectric materials. Journal of Intelligent Material Systems and Structures 21(13): 1263–1278.

Blevins

(1977) Flow-Induced Vibration. New York, NY: Van Nostrand Reinhold.

Covaci

Gontean

(2020) Piezoelectric energy harvesting solutions: A review. Sensors 20(12): 3512.

Dahiya

Valle

(2013) Fundamentals of piezoelectricity. In: Robotic Tactile Sensing. Dordrecht, The Netherlands: Springer. pp.195–245.

Dai

Abdelkefi

Yang

, et al. (2016) Orientation of bluff body for designing efficient energy harvesters from vortex-induced vibrations. Applied Physics Letters 108(5): 053902.

Erturk

Inman

(2011) Piezoelectric Energy Harvesting. New Delhi: John Wiley & Sons. pp.1–15.

Ewere

Wang

(2014) Performance of galloping piezoelectric energy harvesters. Journal of Intelligent Material Systems and Structures 25(14): 1693–1704.

Facchinetti

de Langre

Biolley

(2004) Coupling of structure and wake oscillators in vortex-induced vibrations. Journal of Fluids and Structures 19(2): 123–140.

10.

Firouzi

Abbasi

Sendur

, et al. (2023) Enhancing the performance of Piezoelectric Energy Harvester under electrostatic actuation using a robust metaheuristic algorithm. Engineering Applications of Artificial Intelligence 118: 105619.

11.

Gao

Shih

(2013) Flow energy harvesting using piezoelectric cantilevers with cylindrical extension. IRE Transactions on Industrial Electronics 60(3): 1116–1118.

12.

Hou

Shan

, et al. (2022) A broadband piezo-electromagnetic hybrid energy harvester under combined vortex-induced and base excitations. Mechanical Systems and Signal Processing 171: 108963.

13.

Hou

Shan

Zhang

, et al. (2020) Design and modeling of a magnetic-coupling monostable piezoelectric energy harvester under vortex-induced vibration. IEEE Access 8: 108913–108927.

14.

Jia

Shan

Upadrashta

, et al. (2018) Modeling and analysis of upright piezoelectric energy harvester under aerodynamic vortex-induced vibration. Micromachines 9(12): 667.

15.

Ren

Shang

, et al. (2023) Dynamics and performance evaluation of a vortex-induced vibration energy harvester with hybrid bluff body. Smart Materials and Structures 32(4): 045016.

16.

, et al. (2022a) Geometrically nonlinear model of piezoelectric wind energy harvesters based on vortex-induced vibration and galloping. Smart Materials and Structures 31: 105019.

17.

Liu

, et al. (2022b) Numerical research on a vortex shedding induced piezoelectric-electromagnetic energy harvester. Journal of Intelligent Material Systems and Structures 33(1): 105–120.

18.

Zhang

Litak

, et al. (2024) Periodic solutions and frequency lock-in of vortex-induced vibration energy harvesters with nonlinear stiffness. Journal of Sound and Vibration 568: 117952.

19.

Mackowski

Williamson

CHK

(2013) An experimental investigation of vortex-induced vibration with nonlinear restoring forces. Physics of Fluids 25(8): 087101.

20.

Mahmoodi

Jalili

(2007) Non-linear vibrations and frequency response analysis of piezoelectrically driven microcantilevers. International Journal of Non-Linear Mechanics 42(4): 577–587.

21.

Zhang

, et al. (2023) Dynamic modeling and analysis of a tri-stable vortex-induced vibration energy harvester. Mechanical Systems and Signal Processing 187: 109924.

22.

Mehmood

Abdelkefi

Hajj

, et al. (2013) Piezoelectric energy harvesting from vortex-induced vibrations of circular cylinder. Journal of Sound and Vibration 332(19): 4656–4667.

23.

Nayfeh

Pai

(2004) Linear and Nonlinear Structural Mechanics. New York: John Wiley & Sons Inc.

24.

Metrikine

(2020) A single van der pol wake oscillator model for coupled cross-flow and in-line vortex-induced vibrations. Ocean Engineering 196(196): 106732.

25.

Rezaei

Talebitooti

(2019) Wideband PZT energy harvesting from the wake of a bluff body in varying flow speeds. International Journal of Mechanical Sciences 163: 105135.

26.

Rezaei

Zamanian

(2017) A two-dimensional vibration analysis of piezoelectrically actuated microbeam with nonideal boundary conditions. Physica E Low-dimensional Systems and Nanostructures 85: 285–293.

27.

Sarpkaya

(1979) Vortex-induced oscillations: A selective review. Journal of Applied Mechanics 46: 241–258.

28.

Seyed-Aghazadeh

Samandari

Dulac

(2020) Flow-induced vibration of inherently nonlinear structures with applications in energy harvesting. Physics of Fluids 32(7): 071701.

29.

Skop

Balasubramanian

(1997) A new twist on an old model for vortex-excited vibrations. Journal of Fluids and Structures 11(4): 395–412.

30.

Skop

Griffin

(1973) A model for the vortex-excited resonant response of bluff cylinders. Journal of Sound and Vibration 27(2): 225–233.

31.

Skop

Luo

(2001) An inverse-direct method for predicting the vortex-induced vibrations of cylinders in uniform and nonuniform flows. Journal of Fluids and Structures 15(6): 867–884.

32.

Sui

Zhang

Yang

, et al. (2022) Modeling and experimental investigation of magnetically coupling bending-torsion piezoelectric energy harvester based on vortex-induced vibration. Journal of Intelligent Material Systems and Structures 33(9): 1147–1160.

33.

Sun

Seok

(2019) Development of the optimal bluff body for wind energy harvesting using the synergetic effect of coupled vortex induced vibration and galloping phenomena. International Journal of Mechanical Sciences 156: 435–445.

34.

Lin

(2020) Design and analysis of a vortex-induced bi-directional piezoelectric energy harvester. International Journal of Mechanical Sciences 173: 105457.

35.

Wang

(2021) Development of a non-linear bi-directional vortex-induced piezoelectric energy harvester with magnetic interaction. Sensors 21: 2299.

36.

Wang

Gao

, et al. (2020) Nonlinear analysis of piezoelectric wind energy harvesters with different geometrical shapes. Archive of Applied Mechanics 90: 721–736.

37.

Williamson

CHK

Govardhan

(2004) Vortex-induced vibrations. Annual Review of Fluid Mechanics 36(1): 413–455.

38.

Williamson

CHK

Govardhan

(2008) A brief review of recent results in vortex-induced vibrations. Journal of Wind Engineering and Industrial Aerodynamics 96(6-7): 713–735.

39.

Yang

(2019) A Piezoelectric Wind Energy Harvester with Interaction Between Vortex-Induced Vibration and Galloping. IEEE Sensors, pp.1–4

40.

Younis

Nayfeh

(2003) A study of the nonlinear response of a resonant microbeam to an electric actuation. Nonlinear Dynamics 31: 91–117.

41.

Zamanian

Khadem

(2010a) Nonlinear vibration of an electrically actuated microresonator tuned by combined DC piezoelectric and electric actuations. Smart Materials and Structures 19(1): 015012.

42.

Zamanian

Khadem

(2010b) Analysis of thermoelastic damping in microresonators by considering the stretching effect. International Journal of Mechanical Sciences 52: 1366–1375.

43.

Zhang

Sui

Yang

, et al. (2022) An asymmetric magnetic-coupled bending-torsion piezoelectric energy harvester: modeling and experimental investigation. Smart Materials and Structures 31: 015037.

44.

Zhang

Wang

(2016) Experimental study on piezoelectric energy harvesting from vortex-induced vibrations and wake-induced vibrations. Journal of Sensors 2016: 1–7.

45.

Zhao

Zhang

Wang

(2019) Optimization of galloping piezoelectric energy harvester with V-shaped groove in low wind speed. Energies 12(24): 4619.

46.

Zhao

Tang

Yang

(2012) Small wind energy harvesting from galloping using piezoelectric materials. In: Smart materials, adaptive structures and intelligent systems, Vol. 45103, pp.919–927. American Society of Mechanical Engineers.

47.

Zhao

Tang

Yang

(2016) Synchronized charge extraction in galloping piezoelectric energy harvesting. Journal of Intelligent Material Systems and Structures 27(4): 453–468.