Green׳s function method for piezoelectric energy harvesting beams

Journal of Sound and Vibration 333 (2014) 3092–3108

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Green's function method for piezoelectric energy harvesting beams Amir H. Danesh-Yazdi, Niell Elvin, Yiannis Andreopoulos n Department of Mechanical Engineering, The City College of the City University of New York, New York, NY 10031, United States

a r t i c l e in f o

abstract

Article history: Received 1 August 2013 Received in revised form 7 February 2014 Accepted 19 February 2014 Handling Editor: M.P. Cartmell Available online 27 March 2014

The development of validated mathematical models for piezoelectric harvesters is important as it provides predictive capabilities of their performance and insight to their coupled electromechanical behavior. Advanced solutions to these models allows for more realistic parameters to be considered. In this paper, we present a Fourier Transform– Green's Function (FTGF) solution approach to the distributed parameter coupled electromechanical equations for a piezoelectric beam excited by an arbitrary external transverse force. This method, as opposed to modal analysis, allows for frequency-dependent material properties and damping coefficients to be considered. The special case of a harmonic base excitation is considered and closed-form expressions for the frequency response functions of the voltage generated by piezoelectric layer, relative tip displacement and local bending strain are obtained. Finally, the FTGF solution to these frequency response functions is compared with the modal analysis solution along with experimental data for validation. & 2014 Elsevier Ltd. All rights reserved.

1. Introduction Technological advances over the past two decades have led to a smaller power requirement for microelectronic components. In many cases, batteries may be an adequate electric power supply, but especially for small electronic components, the inconvenience of maintenance (recharge or replacement) along with the additional weight and size of a battery may not be desirable [1]. A potential solution to supplying power to microelectronic components is to harvest the energy available from their ambient environment continuously. The available kinetic (vibration) energy present in an environment requires a transduction mechanism to be converted into electricity. The three most important mechanisms for converting mechanical vibration to electricity are electromagnetic, electrostatic and piezoelectric transductions. The focus of this paper is on piezoelectric transduction, in which the application of mechanical strain to a piezoelectric material creates electric charge on its surface (the so called direct piezoelectric effect). Compared to the electromagnetic and electrostatic transduction mechanisms, piezoelectric energy harvesters have larger power densities [2], relatively high voltages and low currents [3], are not reliant on an external voltage input and are relatively easy to fabricate at the micro- and macro-scales. By far, the most common type of piezoelectric energy harvester is the cantilever beam design, which typically consists of a substructure bonded to one

n

Corresponding author. Tel.: þ1 212 650 5206. E-mail address: [email protected] (Y. Andreopoulos).

http://dx.doi.org/10.1016/j.jsv.2014.02.023 0022-460X & 2014 Elsevier Ltd. All rights reserved.

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3093

(unimorph) or two (bimorph) layers of piezoelectric material, such as PZT (lead zirconate titanate) or PVDF (polyvinylidene fluoride). The development of dependable mathematical models for piezoelectric harvesters is necessary as it allows for the prediction of the behavior of these devices and provides a path to the optimization of the maximum electrical output. The three primary modeling approaches available in the literature are the coupled single degree-of-freedom (SDOF) or lumped parameter models, the approximate distributed parameter models and the (exact) distributed parameter models [4]. The earliest attempt at a mathematical model for the energy harvesters involved modeling both the mechanical equation of motion and the electrical circuit equation as single degree-of-freedom systems [5]. Other researchers have built upon this work by coupling the two ordinary differential equations through an electromechanical coupling constant [6–9]. The advantages of the lumped parameter approach are that the SDOF equations are relatively simple to solve analytically and provide an initial insight into the behavior of the energy harvester problem. One of the greatest disadvantages of the lumped parameter approach though, is that by modeling the mechanical equation of motion as a SDOF system, only a single vibration mode (fundamental mode) is considered and the contribution from the other modes is ignored. Therefore, the SDOF approach may yield inaccurate results, especially when the tip mass attached to the beam is comparatively heavy [10]. An improvement on the lumped parameter approach is the approximate distributed parameter approach in which the mechanical equation of motion and electric circuit equation is derived from Hamilton's principle and the transverse displacement is discretized using the Rayleigh–Ritz approach [1,11]. This method, although approximate, shows good agreement with experimental data. In the distributed parameter approach, Erturk and Inman [4] derived the coupled mechanical equation of motion and electric circuit equation by assuming an Euler–Bernoulli beam and employing the piezoelectric constitutive relations. They solved the coupled electromechanical equations using modal analysis for the base excitation case. It should be pointed out that Erturk and Inman considered the viscous damping contribution due to the fluid surrounding the beam (typically air) and the damping contribution due to the viscoelasticity of the structure separately. The literature on damping of beams, especially viscous air damping, is fairly scarce. Perhaps the most comprehensive study on the physics of the air and internal damping mechanisms for a thin beam was done by Baker et al. [12]. In their paper, they found that material damping, characterized by a logarithmic decrement of free vibration decay in a vacuum, depends on beam thickness and frequency in thin metal beams. They also found that, at small amplitudes, air damping is primarily due to a drag force proportional to the instantaneous velocity and is dependent on air density and the length-tothickness ratio of the beam, but is independent of beam frequency and amplitude. Gibson and Plunkett [13] showed that the material damping, characterized by the loss factor for fiberglass specimens increased with increasing frequency in the lowamplitude range. Rogers et al. [14] found that aluminum and steel specimen beams showed linear damping characteristics, while gray-cast-iron and aluminum-filled-epoxy specimen beams exhibited nonlinear damping characteristics. They also found that air damping was a linear function of amplitude for all specimens. Nashif et al. [15] found that for acrylic, Young's modulus increases slightly and the material damping decreases with increasing frequency. Interestingly, Chiu [16] also showed that the measured Young's modulus of a free-free beam is not substantially influenced by material damping, but rather, the modulus decreases with increasing air damping. In this paper, we propose an alternative approach to solving the coupled distributed parameter electromechanical equations for a piezoelectric beam. The Fourier Transform–Green's Function (FTGF) solution approach is used to solve for the closed-form frequency spectrum expression of the voltage output and local transverse deflection of the beam under an arbitrary, time varying, external load. The frequency response functions (FRF) of the voltage, relative tip deflection of the beam and bending strain distribution are obtained for the case when the beam is under base excitation. Finally, the theoretical FRFs are compared to and validated by an experimental case study. Based on the damping literature [1–15], and to illustrate the power of the FTGF solution approach, in this study, we assume that viscous structural and air damping coefficients are frequency-dependent. This is in contrast to the constant damping coefficients limitation used in modal analysis. As a simple case study, the frequency-dependent structural and air damping coefficients are modeled so that they can be obtained from a curve-fit of the damping ratio data from experiments. 2. Background information 2.1. Modal analysis solution To solve the fully coupled electromechanical equations for a piezoelectric beam (discussed in Section 3), Erturk and Inman [4] employed modal analysis, in which, based on the expansion theorem, the solution to the transverse vibratory response of the system, w(x,t), can be represented by an absolutely and uniformly convergent series of eigenfunctions as: 1

wðx; tÞ ¼ ∑ ϕr ðxÞηr ðtÞ

(1)

r¼1

where ϕr(x) and ηr(t) are the mass normalized eigenfunction and the modal coordinate of the beam for the r-th mode, respectively. Modal analysis is an effective, elegant and easy to use tool that can specifically handle the dynamic response of structures under vibrational excitation in most cases. Modal analysis does have its limitations, though. For instance, modal analysis

3094

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

cannot be used to find the dynamic response in cases where the material properties and/or damping terms are dependent on the frequency of vibration of the beam.

2.2. Passive damping mechanisms and models in piezoelectric beams As a beam oscillates in a fluid medium, various dissipative forces act on the beam and reduce the amplitude of vibration, especially at or near the resonant frequencies. For a piezoelectric beam, these damping forces are primarily active or passive in nature. Active damping requires an external dissipation mechanism, typically by means of a power source. Passive damping, on the other hand, does not require an external power source and is due to the beam itself. More sophisticated damping mechanisms such as semi-passive [17,18] for example, have also been developed and studied, but are beyond the scope of this discussion. In the distributed parameter equation of motion, discussed in Section 3.1, two passive mechanical damping mechanisms are considered: (i) the internal damping due to the viscoelastic nature of the beam and (ii) the damping due to the fluid (usually air) surrounding the beam. These damping terms are considered separately because fluid damping acts on the absolute velocity while strain-rate damping acts on the relative velocity of the beam with respect to the base. Thus, including one damping coefficient is an oversimplification of the physics of the problem [10]. Due to the phenomenological and somewhat unpredictable nature of damping, its determination is considered more a skill based on experience and judgment rather than a rigorously researched and structured science. The amount of damping in a dynamic system is typically taken into account by considering either a loss factor in an SDOF system or damping ratios at each of the modes in a continuous system. The values of the loss factor or damping ratios are obtained either experimentally or are assigned based on experience. In the vast majority of the literature, a damping mechanism is usually assumed to satisfy Rayleigh (proportional) damping. Erturk and Inman [4], for example, assumed constant structural damping cs and air damping ca coefficients in solving the electromechanically coupled equations for the piezoelectric beam. As is the common practice in modal analysis, they lumped together the mechanical damping coefficients in the r-th mode damping ratio, ζr, to give the following relation: ζr ¼

cs Iωr ca þ 2YI 2mωr

(2)

where ωr is the r-th mode angular frequency, m is the mass per unit length of the beam, I is the second moment of area, and Y is Young's modulus. The advantage in expressing the entire damping in the system purely through the damping ratio at each mode is that the voltage and tip displacement FRFs can be expressed directly in terms of the damping ratios, as will be shown in Section 5. The primary disadvantage in this “lumped-damping“ approach is that the underlying physics of damping (i.e. cs and ca) which are inherently part of the original equation of motion (see Eq. (3)) are initially ignored in the modal solution are added back to the solution as essentially curve fit parameters by matching the damping ratios at two modes.

3. Modeling of piezoelectric energy harvesters Consider a cantilevered unimorph piezoelectric energy harvester, as shown in Fig. 1. The piezoelectric and substructure layer are assumed to be perfectly bonded. For the sake of generality, a mass Mt with a mass moment of inertia Jt is attached to the tip of the cantilever beam harvester of length L. The beam is excited by an external transverse force per unit length f(x,t). A resistive electric element Rl is attached to the piezoelectric layer. The composite cantilever is assumed to be an Euler–Bernoulli beam. The damping in the beam is assumed to comprise of internal damping in the structure due to the viscoelasticity of the material (Kelvin–Voigt damping) and the external damping due to fluid (here air) surrounding the beam.

Fig. 1. Unimorph piezoelectric energy harvester with tip mass Mt and mass moment of inertia Jt.

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3095

Fig. 2. Unimorph harvester cross section.

3.1. Mechanical equation of motion with electric coupling Utilizing Euler–Bernoulli beam theory and the piezoelectric constitutive relations, Erturk and Inman [4] found the distributed parameter governing mechanical equation of motion for the beam under general base motion. Replacing the base motion term in their expression with a general external transverse force per unit length f(x,t) (Fig. 1), the equation of motion is: YI


 ∂4 wðx; tÞ ∂2 wðx; tÞ ∂5 wðx; tÞ ∂wðx; tÞ dδðx  x1 Þ dδðx  x2 Þ þϑvðtÞ  ¼ f ðx; tÞ þ c þm þ c I s a ∂t dx dx ∂x4 ∂x4 ∂t ∂t 2

(3)

where w(x,t) is the transverse deflection of the beam, m is the mass per unit length, I is the second moment of area of the composite cross section, cs is the internal (viscous) damping coefficient of the composite due to structural viscoelasticity, ca is the viscous air damping coefficient, YI is the bending stiffness of the composite cross section (Fig. 2), ϑ is the electromechanical coupling term and v(t) is the voltage across the resistor (in Fig. 1). The piezoelectric layer is assumed to be covered with electrodes over the region x1 rx rx2. For the unimorph composite cross section shown in Fig. 2, YI is given by: 3

YI ¼

3

3

3

bs Y s ðhb  ha Þ þ bp Y p ðhc  hb Þ 3

(4)

where bs,bp and Ys,Yp are the widths and Young's moduli of the substructure and piezoelectric layers, respectively, ha is the distance from the neutral axis (N.A.) to the bottom of the substructure layer, hb is the distance from N.A. to the bottom of the piezoelectric layer and hc is the distance from the neutral axis to the top of the piezoelectric layer. The electromechanical coupling term ϑ is given by: ϑ¼ 

Y p d31 bp 2 2 ðhc hb Þ 2hp

(5)

where d31 is the piezoelectric constant and hs and hp are the thicknesses of the substructure and piezoelectric layers respectively, as shown in Fig. 2.

3.2. Electric circuit equation with mechanical coupling Erturk and Inman [4] also obtained the electric circuit equation with mechanical coupling in the piezoelectric layer as: Cp

dvðtÞ vðtÞ þ ¼ β dt Rl

Z

x2 x1

∂3 wðx; tÞ dx ∂x2 ∂t

(6)

where Cp ¼

εS33 bp ðx2 x1 Þ hp

β ¼ Y p d31 bp hpc ¼ 

(7)

2hp hpc 2

2

hc hb

ϑ

(8)

and Cp is the capacitance of the piezoelectric layer, εS33 is the permittivity at constant strain and hpc is the distance from N.A. to the center of the piezoelectric layer. Eqs. (3) and (6) are the distributed parameter electromechanical equations for a piezoelectric energy harvester beam in transverse vibration which will be solved in this paper using a novel approach. It is important to note that Eqs. (3) and (6) would be similar for a bimorph beam.

3096

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3.3. Boundary conditions Assuming the most general case of a cantilever beam with fixed end at x¼0 and a tip mass Mt with a mass moment of inertia Jt attached at x ¼L (Fig. 1), the boundary conditions for the equation of motion are: wð0; tÞ ¼ 0  ∂wðx; tÞ ∂x 

x¼0

¼0

(9) (10)


 ∂2 wðx; tÞ ∂3 wðx; tÞ ∂3 wðx; tÞ þ Jt þ cs I ¼0 YI 2 2 2 ∂x ∂x ∂t ∂t ∂x x ¼ L

(11)


 ∂3 wðx; tÞ ∂4 wðx; tÞ ∂2 wðx; tÞ  M þ c I ¼0 YI s t ∂x3 ∂x3 ∂t ∂t 2 x¼L

(12)

4. Fourier transform–Green's function (FTGF) solution approach As an alternative to modal analysis, the linear differential Eqs. (3) and (6) are solved using a two-step approach: 1. Fourier Transform: Applying a Fourier transform in time to Eqs. (3) and (6) reduces the two fully coupled partial differential equations to two fully coupled ordinary differential equations. 2. Green's Function: Applying a Green's function to the transformed mechanical equation of motion decouples the voltage and deflection. The Green's function [19,20], is a mathematical tool first incorporated by George Green in 1828, which provides a method for solving nonhomogeneous linear differential equations. The application of a Green's function approach to solve the Euler–Bernoulli beam equation in our paper is inspired by the work of Sader [21].

4.1. Fourier transform of the coupled electromechanical equations The Fourier transform in time for a function g(t) is defined as: Z 1 ^ gðωÞ  gðtÞe  iωt dt

(13)

1

Taking the Fourier transform of differential Eqs. (3) and (6) gives:

 ^ f^ ðx; ωÞ ∂4 wðx; mω2  iωca ðωÞ ϑv^ ðωÞ dδðx  x1 Þ dδðx  x2 Þ ωÞ ^  ¼ wðx; ωÞ þ  4 YIðωÞ þ iωcs IðωÞ YIðωÞ þ iωcs IðωÞ YIðωÞ þ iωcs IðωÞ dx dx ∂x v^ ðωÞ ¼ 


x ¼ x2 ^ iωβ ∂wðx; ωÞ 1=Rl þ iωC p ∂x x ¼ x1

(14)

(15)

The Fourier transformed boundary conditions are: ^ wð0; ωÞ ¼ 0

(16)

 ^ ∂wðx; ωÞ ¼0 ∂x x ¼ 0

(17)


 ^ ^ ωÞ ∂2 wðx; 2 ∂wðx; ωÞ J ω ¼0 ðYIðωÞ þiωcs IðωÞÞ t ∂x ∂x2 x¼L

(18)


 ^ ωÞ ∂3 wðx; ^ þ M t ω2 wðx; ¼0 ωÞ ðYIðωÞ þ iωcs IðωÞÞ 3 ∂x x¼L

(19)

It is noteworthy to mention that after applying the Fourier transforms, Young's moduli and damping coefficients can be considered more generally to be functions of frequency ω without significantly influencing the mathematical solution approach or the linearity of the governing differential equations. Using this fact, the bending stiffness YI(ω), the equivalent structural damping term csI(ω) and the viscous air damping coefficient ca(ω) in Eqs. (14)–(19) are all assumed to be frequency-dependent. In this section, Eqs. (3) and (6) have been transformed into two fully coupled ordinary differential Eqs. (14) and (15). As will be shown in the next section, a Green's function approach can be applied to Eq. (14) and, in conjunction with (15), lead to closed-form solutions for the voltage and beam deflection.

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3097

4.2. Green's function application to the equation of motion Replacing x by an auxiliary space coordinate ξ, multiplying (14) by an as-of-yet undetermined Green's function G(x,ξ;ω) and integrating from 0 to L gives:
Z L Z L ^ ωÞ ∂4 wðξ; mω2  iωca ðωÞ ϑv^ ðωÞ ^ Gðx; ξ; ωÞ dξ  Gðx; ξ; ωÞwðξ; ωÞdξ þ 4 YIðωÞ þ iωc YIðωÞ þ iωcs IðωÞ IðωÞ ∂ξ s 0 0
 Z L Z L dδðξ  x1 Þ dδðξ  x2 Þ 1  dξ  Gðx; ξ; ωÞ Gðx; ξ; ωÞf^ ðξ; ωÞdξ ¼ 0 (20)  dξ dξ YIðωÞ þ iωcs IðωÞ 0 0 Simplifying (20) and rearranging the terms results in:
4
  Z L ∂ Gðx; ξ; ωÞ mω2  iωca ðωÞ ϑv^ ðωÞ ^ Gðx; ξ; ωÞ dξ þ  wðξ; ωÞ 4 YIðωÞ þ iωc YIðωÞ þiωcs IðωÞ IðωÞ ∂ξ s 0
 Z L Z L dδðξ  x1 Þ dδðξ x2 Þ 1  dξ  Gðx; ξ; ωÞ Gðx; ξ; ωÞf^ ðξ; ωÞdξ ¼ 0  dξ dξ YIðωÞ þ iωcs IðωÞ 0 0 Note that an identity associated with the Dirac delta function is that: Z L ^ ^ wðξ; ωÞδðξ  xÞdξ ¼ wðx; ωÞ

(21)

(22)

0

Matching the identity in (22) with the first term in (21) leads to the Green's function differential equation:
 ∂4 Gðx; ξ; ωÞ mω2  iωca ðωÞ Gðx; ξ; ωÞ ¼ δðξ  xÞ  YIðωÞ þiωcs IðωÞ ∂ξ4 The solution to the Green's function differential equation is derived in detail in the Appendix. Thus, (21) reduces to: " #
 Z L 1 ∂Gðx; ξ; ωÞ ξ ¼ x2 ^ ϑv^ ðωÞ  Gðx; ξ; ωÞf^ ðξ; ωÞdξ wðx; ωÞ ¼  YIðωÞ þiωcs IðωÞ ∂ξ 0 ξ ¼ x1

(23)

(24)

Taking the derivative of (24) with respect to x and evaluating from x¼ x1 to x ¼x2 gives: "
#x ¼ x2
x ¼ x2 
Z L x ¼ x2 ^ ∂wðx; ωÞ ϑv^ ðωÞ ∂ ∂Gðx; ξ; ωÞ ξ ¼ x2 1 ∂ ¼ þ Gðx; ξ; ωÞf^ ðξ; ωÞdξ ∂x YIðωÞ þ iωcs IðωÞ ∂x ∂ξ YIðωÞ þ iωcs IðωÞ ∂x 0 x ¼ x1 ξ ¼ x1 x ¼ x1

(25)

x ¼ x1

Using expression (25), Eq. (15) can now be fully solved for v^ ðωÞ: hR ix ¼ x2 L ∂ ^ iωβ ∂x 0 Gðx; ξ; ωÞf ðξ; ωÞdξ x ¼ x 1 v^ ðωÞ ¼  h ix ¼ x2 x2 ∂ ð1=Rl þ iωC p ÞðYIðωÞ þ iωcs IðωÞÞ  iωβϑ ∂x ½ð∂Gðx; ξ; ωÞÞ=∂ξξξ ¼ ¼ x1

(26)

x ¼ x1

Using the result in (26), Eq. (14) reduces to: hR ix ¼ x2 L x2 ^ iωβϑ½ð∂Gðx; ξ; ωÞÞ=∂ξξξ ¼ ¼ x1 ∂=∂x 0 Gðx; ξ; ωÞf ðξ; ωÞdξ 1 x ¼ x1 ^ wðx; ωÞ ¼ h i YIðωÞ þ iωcs IðωÞ ð1=R þ iωC ÞðYIðωÞ þ iωc IðωÞÞ iωβϑ ∂ ½ð∂Gðx; ξ; ωÞÞ=∂ξξ ¼ x2 x ¼ x2 p s l ξ ¼ x1 ∂x 1 þ YIðωÞ þ iωcs IðωÞ

x ¼ x1

Z

L 0

Gðx; ξ; ωÞf^ ðξ; ωÞdξ

(27)

Eqs. (26) and (27) give the frequency spectrum of the voltage and local transverse beam deflection. In the present approach, we have obtained a closed-form solution for the displacement w(x,t) and the voltage output v(t) ^ in the Fourier space wðx; ωÞ and v^ ðωÞ, respectively, as an explicit function of the transformed forcing function f^ ðx; ωÞ. If the time histories of the voltage and deflection are desired, the inverse Fourier transform can be applied as follows: Z 1 1 vðtÞ ¼ (28) v^ ðωÞeiωt dω 2π  1 wðx; tÞ ¼

1 2π

Z

1 1

^ wðx; ωÞeiωt dω

(29)

3098

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

4.3. Advantages of the FTGF method Compared to the modal analysis solution method utilized by Erturk and Inman, the FTGF method has some distinct advantages:

 The FTGF method provides an exact solution to the fully coupled electromechanical equations for a piezoelectric beam, while the modal analysis approach yields an exact solution only when the damping is proportional.

 Eqs. (26) and (27) derived using the FTGF method can handle cases where the damping coefficients (cs and ca) and/or 



Young's moduli (Yp and Ys) are functions of frequency ω. The modal analysis approach cannot be applied when the material properties and damping coefficients are frequency-dependent. In modal analysis, for a fully coupled system of differential equations such as the electromechanical equations for a piezoelectric beam, the time response of at least one of the variables has to be assumed. Erturk and Inman [4], for example, assumed a harmonic voltage response to their solution for a piezoelectric beam excited by harmonic base motion. For the FTGF solution method, no assumption is made about the voltage and beam deflection response. The excitation force on the beam dictates the form of the voltage output due to the linear nature of Eq. (3). In the modal analysis approach, this is a limiting factor in the type of cases that can be considered since, as mentioned previously, an assumption has to be made about the time response of the voltage. Therefore, random excitation functions such as those due to turbulent flows become difficult to analyze. On the other hand, virtually any transverse excitation force can be directly placed into Eqs. (26) and (27) to yield the voltage and beam deflection response in the FTGF solution.

5. Case studies  base excitation frequency response function One of the most frequently used study cases in literature is the harvester cantilever beam is excited due to the motion of its base. From Timoshenko [22], if the base displacement is nonzero, the absolute transverse displacement of the beam can be expressed as: wðx; tÞ ¼ wb ðx; tÞ þwrel ðx; tÞ

(30)

where wb(x,t) is the base displacement and wrel(x,t) is the relative transverse displacement. Assuming that the base is under general harmonic translation, we have: wb ðtÞ ¼ Aeiωb t

(31)

where A is the amplitude of the displacement and ωb is the applied base frequency. From (31), the base acceleration can be expressed as: ab ðtÞ ¼  ω2b Aeiωb t

(32)

5.1. Unimorph cantilever beam with tip mass and moment of inertia For a unimorph cantilever beam with a tip mass Mt and mass moment of inertia Jt subjected to a base acceleration, the external forcing function reduces to: f ðx; tÞ ¼ ð½m þM t δðx  LÞω2b  ica ωb ÞAeiωb t

(33)

Replacing the Fourier transform of the forcing function (33) into Eqs. (26) and (27) and applying the inverse Fourier transform gives the voltage (per base acceleration) and the deflection (per base displacement) frequency response functions (FRFs): vðtÞ ¼ ab ðtÞ

hR L

ix ¼ x2

x2 þiM t ωb β½∂Gðx; L; ωb Þ=∂xxx ¼ ¼ x1 x ¼ x1 ξ ¼ x x2 ∂ ð1=Rl þ iωb C p ÞðYI þiωb cs IÞ  iωb βϑ ∂x ½δGðx; ξ; ωb Þ=δðξÞξ ¼ x21 xx ¼ ¼ x1

∂ iðmωb  ica Þβ ∂x

0

Gðx; ξ; ωb Þdξ

(34)


 Z ϑω2b mω2b  ica ωb L wrel ðx; tÞ ∂Gðx; ξ; ωb Þ ξ ¼ x2 vðtÞ ¼ þ Gðx; ξ; ωb Þdξ wb ðtÞ ∂ξ YI þ iωb cs I YI þ iωb cs I 0 ξ ¼ x1 ab ðtÞ þ

M t ω2b Gðx; L; ωb Þ YI þ iωb cs I

(35)

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3099

5.2. Unimorph cantilever beam with free end For a unimorph cantilever beam with a free end, Mt ¼0. Thus, the voltage (per base acceleration) and the relative deflection (per base displacement) frequency response functions (FRFs) become: hR ix ¼ x2 L ∂ iðmωb  ica Þβ ∂x Gðx; ξ; ωb Þdξ 0 vðtÞ x ¼ x1 ¼ (36) h ix ¼ x2 ab ðtÞ ð1=R þ iω C ÞðYI þ iω c IÞ iω βϑ ∂ ½∂Gðx; ξ; ωÞ=∂ξξ ¼ x2 l

b

p

b s

b

∂x

ξ ¼ x1

x ¼ x1


 Z mω2b  ica ωb L wrel ðx; tÞ ∂Gðx; ξ; ωb Þ ξ ¼ x2 vðtÞ ¼ þ Gðx; ξ; ωb Þdξ wb ðtÞ ∂ξ YI þ iωb cs I YI þ iωb cs I 0 ξ ¼ x1 ab ðtÞ ϑω2b

(37)

The frequency response of the bending strain in the beam, S1(x,t), per base acceleration is also often of interest as it provides an expression for comparing the theoretical solution with strain gauge data from experiment. The bending strain at y¼ha is: " #
3  Z ϑω2b mω2b ica ωb L ∂2 Gðx; ξ; ωb Þ S1 ðx; tÞ ha ∂ Gðx; ξ; ωb Þ vðtÞ ¼ 2 þ dξ (38) ab ðtÞ YI þiωb cs I 0 ∂x2 ∂ξ ∂x2 ωb YI þ iωb cs I ξ ¼ L ab ðtÞ

5.3. FTGF solution validation As a validation, we will compare the FTGF solution with the modal analysis FRF results as derived by Erturk and Inman [4]. The voltage per base acceleration and the relative tip deflection per base displacement FRFs are: 2 2 ∑1 ð  imωb φr γ w vðtÞ r Þ=ðωr  ωb þi2ζ r ωr ωb Þ ¼ 1 r¼1 ab ðtÞ ∑r ¼ 1 iωb χ r φr =ðω2r  ω2b þi2ζ r ωr ωb Þ þ ðð1 þ iωb τc Þ=τc Þ

(39)


  1 wrel ðx; tÞ χ r vðtÞ mωb 2 ϕr ðxÞ ¼ ∑ γw r  wb ðtÞ m ab ðtÞ ω2r  ωb 2 þ i2ζ r ωr ωb r¼1

(40)

where φr ¼ 

 d31 Y p hpc hp dϕr ðxÞx ¼ x2 εs33 ðx2  x1 Þ dx x ¼ x1 Z γw r ¼

ϕr ðxÞdx

(42)

x ¼ x2 dϕ ðxÞ χ r ¼ ϑ r  dx x ¼ x1

(43)

τc ¼

ϕr ðxÞ ¼

L

(41)

0

Rl εS33 bp ðx2  x1 Þ hp

(44)

rffiffiffiffiffiffiffi    1 λr x λr x sinh λr  sin λr λr x λr x  cos   sin sinh cosh mL L L L L cosh λr þ cos λr

Table 1 Geometrical, material and electromechanical parameters of the unimorph cantilever beam from Erturk and Inman [4]. Length of the beam, L [mm] Width of the beam, b [mm] PZT electrode start distance from the base, x1 [mm] PZT electrode end distance from the base, x2 [mm] Thickness of the substructure layer, hs [mm] Thickness of the piezoelectric layer, hp [mm] Young's modulus of the substructure layer, Ys [GPa] Young's modulus of the piezoelectric layer, Yp [GPa] Equivalent damping due to structural viscoelasticity, csI/YI [s rad  1] Damping due to air, ca/m [rad s  1] Density of the substructure layer, ρs [kg m  3] Density of the piezoelectric layer, ρp [kg m  3] Piezoelectric constant, d31 [pm V  1] Permittivity, εs33 [nF m  1]

100 20 0 100 0.5 0.4 100 66 1.2433  10  5 4.886 7165 7800  190 15.93

(45)

3100

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

Fig. 3. (a) Voltage and (b) relative tip displacement comparison between the FTGF and Erturk and Inman (modal analysis) solutions for Rl ¼ 1 MΩ.

ωr ¼

λ2r

sffiffiffiffiffiffiffiffiffi YI mL4

(46)

and the dimensionless frequency numbers λr are determined from: 1 þ cos λr cosh λr ¼ 0 The frequency response of the bending strain S1(x,t) per base acceleration at y¼ha is found from Eq. (40):
  2 ω2b S1 ðx; tÞ ha 1 vðtÞ d ϕr ðxÞ ¼ 2 ∑ mγ w r  χr 2 2 ab ðtÞ ab ðtÞ ωr  ωb þ i2ζ r ωr ωb dx2 ωb r ¼ 1

(47)

(48)

Using the material, geometrical and electromechanical parameters of Erturk and Inman's piezoelectric cantilever beam listed in Table 1, we validate the FTGF approach by comparing the voltage and relative tip displacement FRFs to their modal analysis solution with three and ten modes, as shown in Fig. 3. It is important to note that since the FTGF solution is exact, the modal analysis solution approaches this result under the proportional damping assumption by summing over more modes in Eq. (1). 6. Case studies – passive damping modeling One of the main advantages of the FTGF solution technique compared to the modal analysis solution is its ability to handle cases where the damping coefficients cs and ca and/or Young's moduli Yp and Ys are dependent on the vibration frequency of the beam. Our focus here is on the frequency-dependence of the damping coefficients. It is noteworthy to mention that several researchers have attempted to modify the modal analysis solution to include the effects of frequency on the damping ratio [23] and dynamic stiffness [24]. However, most of the research focuses on discrete systems, while we consider a continuous beam. The FTGF approach is ideally suited for cases where a model for the frequency-dependent damping coefficient exists since the method can be readily extended to provide a solution for the voltage and beam deflection responses. However, in most cases, a preliminary experiment is typically necessary to estimate the effects of damping forces on the physics of the problem. As a sample case study, we assume Young's moduli Yp and Ys to be constant and, in order to maintain the linear nature of Eq. (3), we assume the damping coefficients cs and ca are functions of the angular frequency, ω, only. Since the purpose of this paper is to introduce a new approach to solve Eqs. (3) and (6) and not necessarily to study the effects of each individual damping term, we generate functions cs(ω) and ca(ω) by curve-fitting the experimental data. In this approach, Eq. (2) is modified so that the damping ratio ζ is a continuous function of ω and the following definitions for the damping terms are used: cs ðωÞI ¼ K 1

YI ω

ca ðωÞ ¼ K 2 ðωÞmω

(49) (50)

where K1 is a constant and K2 (ω) is an as-of-yet undetermined function of ω. Replacing (49) and (50) into (2) results in: K 1 þ K 2 ðωÞ ¼ 2ζðωÞ

(51)

The unknown terms K1 and K2 (ω) can be obtained from a curve-fit of the damping ratio data obtained from experiments. The term K1 corresponds to the constant in the curve-fit equation and K2 (ω) corresponds to the parts of the curve-fit equation that are functions of ω. It is once again important to stress that the damping relations (49) and (50) may not be representative of the behavior of the actual damping coefficients and are only chosen in this manner for mathematical convenience.

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3101

7. Experimental setup and techniques The theoretical FTGF solution approach and the frequency-dependent passive damping coefficients were also validated experimentally and compared to the modal analysis solution. The experimental setup is shown in Fig. 4. The unimorph cantilever beam consists of a thin piezoelectric PVDF layer from MSI-USA, Inc. of dimensions 30 mm  12.7 mm  28 μm bonded on top of a much thicker acrylic beam of dimensions 20.1 cm  15.9 mm  1.65 mm, as shown in Fig. 5. The beam is held in place by an aluminum grip, which is bolted onto the mechanical shaker. A quarterbridge strain gauge is bonded to the bottom of the acrylic beam and is located 4 mm from the base. A PCB393B04 accelerometer (PCB Piezotronics) is placed on top of the grip. Also, a Micro-Epsilon optoNCDT 1402-10 proximity sensor is positioned so as to measure the tip displacement of the beam. The experiment involves performing a frequency sweep over the first three transverse vibration modes of the beam under harmonic base excitation, with more data taken at the resonant frequencies than the frequencies in between modes. At each frequency, the voltage response of the accelerometer (representing base acceleration), the strain gauge (representing the strain at xs ¼4 mm from the base), the piezoelectric material (PVDF) and proximity sensor (representing the absolute tip displacement) were recorded. As an example, the calibrated accelerometer, strain, tip displacement and PVDF voltage signals are shown in Fig. 6 at the 2nd resonance mode (82.00 Hz). Moreover, we also captured the first two resonance frequencies with a v710 Phantom high-speed camera (Vision Research), as shown in Fig. 7. For the unimorph cantilever beam discussed in the experiment, the geometrical, material and electromechanical parameters are listed in Table 2. The manner in which each value is obtained is also indicated. In Table 2, the values of the geometrical parameters, load resistance and capacitance of the PVDF are measured, the material and electromechanical properties of the PVDF layer are obtained from the manufacturer's datasheet and Young's modulus of the acrylic beam is fitted to the first resonance mode of the bending strain per base acceleration amplitude, |S1(xs,t)/ab(t)|, FRF. It should be noted that acrylic's density and modulus can fall within a relatively wide range and thus one of these parameters need to be experimentally curve-fitted. 8. Results and discussions In this section, we present the comparison between the voltage, tip deflection and bending strain FRFs obtained from the theoretical FTGF solution given by Eqs. (36, 37 and 38), theoretical modal analysis solution given by Eqs. (39, 40, 48) and experimental results for a cantilever beam vibrating under harmonic base excitation.

Fig. 4. Experimental setup for cantilever beam vibration harvesting.

Fig. 5. Schematic layout for the cantilever energy harvester experiment with relevant dimensions (in mm).

3102

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

Fig. 6. (a) Base acceleration, (b) strain at xs, (c) PVDF voltage and (d) absolute tip displacement signals for a vibration frequency of 82.00 Hz.

Fig. 7. Beam vibration at (a) 1st resonance frequency at 12.97 Hz and (b) 2nd resonance frequency at 82.00 Hz.

Table 2 Geometrical, material and electromechanical parameters of the unimorph cantilever beam: (c) calculated, (d) manufacturer's datasheet, (f) fitted to experiment, (m) measured. Length of acrylic beam, L [cm] Width of acrylic beam, bs [mm] Width of PVDF layer, bp [mm] Thickness of the acrylic beam, hs [mm] Thickness of PVDF layer, hp [mm] Distance of the start of PVDF from the base, x1 [mm] Distance of the end of PVDF from the base, x2 [mm] Distance of the center of the strain gauge from the base, xs [mm] Young's modulus of acrylic, Ys [GPa] Young's modulus of PVDF, Yp [GPa] Density of acrylic, ρs [kg m  3] Density of PVDF, ρp [kg m  3] Piezoelectric constant of PVDF, d31 [pm V  1] Permittivity of PVDF at constant strain, εs33 [pF m  1] Internal capacitance of PVDF layer, Cp [nF] Load resistance, Rl [MΩ]

20.1(m) 15.9(m) 12.7(m) 1.65(m) 28(m) 10(m) 40(m) 4(m) 4.8(f) 3.0(d) 1220(d) 1780(d) 25(d) 106.2(d) 1.45(c) 1.3(m) 10(m)

8.1. Damping ratio and damping coefficients determination Before proceeding with comparing the theoretical and experimental FRFs, we will provide an estimate of the amount of damping at each mode. In our study, the damping ratio at each mode, ζr, is determined experimentally from the bending strain per base acceleration amplitude, |S1(xs,t)/ab(t)| FRF using the half-power bandwidth method (Nashif et al. [15], for example). The damping ratios for the first three modes and the resulting curve-fit equation are shown in Fig. 8. The modal frequencies and damping ratios for the three modes are shown in Table 3.

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3103

Fig. 8. Damping ratio with respect to the angular frequency for the first three modes.

Table 3 Mode frequency and damping ratio. Mode

Frequency [Hz]

Damping ratio [%]

1 2 3

12.85 82.0 234.5

4.94 3.65 2.69

In our frequency-dependent passive damping model, the damping ratio Eq. (51) is assumed to contain both a constant K1 and a frequency-dependent K2 (ω) component. Due to the behavior of the modal damping ratios in Fig. 8, we fit a curve to the data of the form: ζðωÞ ¼ a0 þ a1 e  a2 ω where a0, a1, and a2 are constants. In this experiment, a0 ¼0.02456, a1 ¼0.02846 and a2 ¼0.00169 s rad (52) and (51) and substituting the resulting expressions K1 and K2 (ω) into Eqs. (49), (50) results in:

(52) 1

. Matching Eqs.

cs ðωÞI ¼ ð0:0491YIÞω  1

(53)

ca ðωÞ ¼ ð0:0569 mÞωe  0:00169ω

(54)

As mentioned in Section 6, the frequency-dependence of the air and structural damping coefficients constructed here is purely for mathematically convenience. Admittedly, Eqs. (53) and (54) may not be physically valid over the entire frequency range, as demonstrated by the fact that, for example, the structural damping coefficient is unbounded at very low frequencies. 8.2. Voltage per base acceleration FRF One of the simplest FRFs that can be obtained experimentally is the magnitude of the voltage across the resistive load Rl per base acceleration, |v(t)/ab(t)|. The comparison between the magnitude of the experimental |v(t)/ab(t)| FRF and the theoretical FTGF solution FRF (36) with frequency-dependent damping coefficients (53), (54) and the theoretical modal analysis solution (39) with the modal damping ratios from Table 3 is shown in Fig. 9. The agreement between the theoretical approaches at all three modes is very good, because our frequency-dependent damping coefficients are modeled based on the curve-fitting of the damping ratios, meaning that at the modes, the damping coefficients correspond exactly to the damping ratios in Table 3. The comparison between the magnitude of the voltage per base acceleration FRF at each mode is presented in Table 4. From Table 4, the FTGF theoretical solution overestimates the voltage by 16, 21 and 80 percent, correspondingly, at the three resonant frequencies. This is expected since the stress transmitted between the bonding layer that attaches the PVDF layer to the acrylic beam is ignored in the distributed parameter model, i.e. the PVDF and acrylic are assumed to be perfectly bonded [25]. 8.3. Bending strain per base acceleration FRF Another FRF that can be obtained experimentally is the magnitude of the bending strain at x¼ xs per base acceleration, |S1(xs,t)/ab(t)|, using the strain gauge and accelerometer data. The comparison between the magnitude of the experimental FRF and the theoretical FTGF solution FRF (38) and the theoretical modal analysis solution (48) is shown in Fig. 10. It is important to note that the experimental strain FRF data at 40 and 160 Hz are not included due to the very low signal-tonoise ratio of the strain signal at these frequencies.

3104

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

Fig. 9. Voltage per base acceleration FRF for the first three modes. Table 4 Magnitude of the voltage per base acceleration comparison between the experimental data and FTGF solution. Mode

|v(t)/ab(t)| (experiment) [V s2 m  1]

|v(t)/ab(t)| (FTGF) [V s2 m  1]

1 2 3

1.98 0.155 0.0065

2.305 0.188 0.0117

The agreement between the theoretical approaches at all three modes is excellent, similar to the voltage FRF case. The theoretical FTGF and modal analysis solutions also match the experimental data extremely well at the first mode. The magnitudes of the theoretical FRFs are also in good agreement with the experimental data at the second and third modes, but, as in the voltage to base acceleration FRF, there is a clear shift in frequency. The theoretical modes occur at a lower frequency compared to the experimental modes. This frequency shift, which becomes larger at higher modes, can primarily be attributed to the fact that in both the FTGF and modal analysis solution, Young's modulus of the acrylic beam and the PVDF layer, Ys and Yp, are assumed to be constant. As mentioned in Section 1, Young's modulus of acrylic, for example, increases slightly with increasing frequency. The mode occurs at a higher frequency with a larger Young's modulus, which would allow for both the second and third modes to be matched exactly. In our case, to prevent overcomplicating the problem, Young's moduli are assumed to be constant, but as mentioned previously, the FTGF solution approach can also handle frequency-dependent Young's moduli; this is not the case for modal analysis. Overall, though, the FTGF solution is in very good agreement with the experimental data considering the fact that a constant Young's modulus for the acrylic at the first mode is the only parameter that is fitted to the data. It should be noted that as shown in Fig. 1, the piezoelectric layer is assumed to cover the entire beam length (even if the electrodes do not). However, for the case of a thin piezoelectric layer (as is the case in our experiment), including or neglecting the piezoelectric layer provides a negligible difference in bending stiffness, and thus the derived equations can be used directly. For the parameters mentioned in Table 2, the bending stiffness from Eq. (4) is YI ¼0.0287 N m2; if the effect of the piezoelectric layer is ignored, the bending stiffness is YI¼ 0.0279 N m2. The difference between the two values is less than 3 percent and the voltage and bending strain FRFs are not discernibly influenced by using either value. Therefore, for this experiment, Eq. (4) provides a reasonable estimate of the bending stiffness for the entire beam.

8.4. Theoretical relative tip displacement per base displacement FRF comparison for different damping models Finally, we consider the magnitude of the relative tip displacement per base displacement FRF, |wrel(L,t)/wb(t)|. Due to the difficulty in experimentally measuring the relative displacement of the beam's tip, we will focus only on the theoretical FRFs. The behavior of the relative tip displacement FRF depends on the damping model that is employed in the theory. In Fig. 11, we compare the following theoretical FRFs: (i) FTGF solution (37) with frequency-dependent damping coefficients (53), (54), (ii) modal analysis solution (40) using the damping ratios for the first three modes determined experimentally, (iii) the modal analysis solution (40) for the first ten modes determined from the curve-fit damping ratio Eq. (52), and (iv) FTGF solution (37) with Rayleigh (proportional) damping using the first two experimentally determined modes (i.e. csI ¼3.35  10  6 kg m3 s  1 and ca ¼0.2357 kg m  1 s  1 from Eq. (2)). The agreement between all four solutions is excellent for the first two modes. At the third mode, though, the FTGF solution with Rayleigh damping differs from the other three solutions, which is expected, since only the first two modes are matched exactly in proportional damping. It is clear from Fig. 11 that, for this experimental setup, the proportional damping model overestimates the damping ratio at the third mode. The major difference between the four theoretical solutions occurs at the off-resonance minimum amplitude frequency between the second and third modes. This is referred to as the antiresonance frequency, at which the relative tip displacement is less than 1 percent of the base displacement [4]. Comparing the FTGF solution with frequency-dependent

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3105

Fig. 10. Bending strain per base acceleration FRF for the first three modes.

Fig. 11. Relative tip displacement per base displacement FRF comparison for the first three modes using different damping models.

damping coefficients and the modal analysis solution with three modes, we see that the FTGF solution predicts a much smaller antiresonance magnitude. When ten modes (by extrapolating the damping ratio values from Eq. (52)) are considered in the modal analysis solution, the magnitude of the FRF decreases, but there is a pronounced shift in the antiresonance frequency. The behavior of these three cases can be compared to their proportional damping counterpart, which we depicted in Fig. 3(b) as a validation for the FTGF method. In Fig. 3(b), the FTGF solution predicts a smaller antiresonance magnitude compared to the modal analysis solution with three modes, just like in Fig. 11. However, when ten modes are considered in Fig. 3(b), the modal analysis solution virtually matches the FTGF solution, even at the antiresonance frequency. This difference between the theoretical FRFs with respect to proportional and non-proportional damping can be traced to eq. (2), which relates the damping ratio at each mode to the damping coefficients. Under the proportional damping assumption, cs and ca are constants that are determined from the frequencies and damping ratios from the first two modes. The damping ratios of the higher modes are then automatically set by Eq. (2). The modal analysis solution, as defined in Eq. (1), provides an exact solution when infinitely-many modes are known, but considering a finite number of modes suffices. The FTGF solution is exact for a known cs and ca. Since cs and ca set the higher mode damping ratios in proportional damping, the modal analysis and FTGF solutions would match if enough modes are considered. When the damping is non-proportional, though, the modal analysis and FTGF approach lead to different solutions, especially at off-resonance frequencies. In modal analysis, as the name implies, the only pertinent information to the solution is the frequency and damping ratio at the modes. If the damping ratios at higher modes are determined experimentally, they can be applied into Eqs. (39), (40) and (48) directly. The resulting FRFs will yield the correct magnitude at the modes, but the values at off-resonance frequencies may not correspond to experimental results, as demonstrated by the modal analysis solution with ten modes in Fig. 11. On the other hand, the FTGF approach is a direct solution to the electromechanical equations and as long as the cs(ω) and ca(ω) functions correspond to the experimentally determined coefficients, the solution is expected to be exact, even at the off-resonance frequencies. This is the reason that we expect the experimental data to show better agreement with the FTGF solution with frequency-dependent damping coefficients than with the other three solutions in Fig. 11. Accurate experimental results at the antiresonance frequencies (which is beyond our current experimental capability) need to be performed to validate the damping models especially at the anti-resonances.

9. Conclusions In the present work, we have developed and demonstrated the applicability of the Fourier Transform–Green's Function (FTGF) solution approach to the distributed parameter coupled electromechanical equations for a piezoelectric cantilever beam under transverse vibration. In the most general case, the Euler–Bernoulli beam is assumed to be excited by an arbitrary external transverse force with a tip mass and a mass moment of inertia. The viscous air and internal damping

3106

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

coefficients are considered separately in the mechanical equation of motion and, along with Young's moduli of the piezoelectric and substructure layer, are assumed to be frequency-dependent in this solution method. The FTGF approach involves first taking the Fourier transform of the mechanical equation of motion, the electric circuit equation and the boundary conditions. A Green's function is then applied to the transformed equation of motion. This allows us to obtain the frequency spectrum of the voltage generated by the piezoelectric layer and the local deflection of the beam. The time histories of both the voltage and beam deflection can be then found by applying an inverse Fourier transform. The special case of a harmonic base excitation is considered and closed-form expressions for the FRFs of the voltage per base acceleration, relative tip displacement per base displacement and the bending strain per base acceleration are obtained. The validity of the FTGF approach is verified experimentally. The experimental setup consisted of a thin PVDF layer bonded to an acrylic cantilever beam and excited by a mechanical shaker. An accelerometer was used to measure the base acceleration of the beam and a strain gauge was attached to the bottom of the beam close to the base. The experiment involved performing a frequency sweep over the first three modes of the beam under harmonic transverse base excitation. At each frequency, the voltage response of the accelerometer at the base of the beam, the strain gauge, the PVDF layer and proximity sensor were recorded. In our experimental case study, we assume Young's moduli to be constant and the viscous air and internal damping coefficients to be modeled as functions of the angular frequency and related to a curve-fit equation of damping. Young's moduli are assumed constant to prevent further curve-fitting, but, as mentioned, they can be easily expanded to be frequency-dependent. The damping ratios for the first three modes are determined from the magnitude of the experimental bending strain per base acceleration using the half-power bandwidth method. The damping ratios are then curve-fitted to an equation as a function of the angular frequency and this equation is then used to obtain explicit relations for the viscous air and internal damping coefficients based on the aforementioned mathematical model. In this study, we also compare the voltage and bending strain FRFs obtained from the theoretical FTGF and modal analysis solutions with experimental results. In both FRFs, we observe that the two theoretical solutions are in very good agreement at the three modes, but are discernibly different at off-resonance frequencies. Another observation is that a shift in the frequency between the theoretical results and the experimental data is in the second and third modes, which also indicates that the modulus of the acrylic beam is frequency dependent. We also discuss the reason behind why the theoretical models overestimate the amplitude of the voltage FRF compared to experimental data. The theoretical relative tip displacement FRFs are also compared with each other for different damping models and the differences between the FTGF and modal analysis solution at the antiresonance frequency is discussed and explained. Finally, the FTGF approach needs to be considered in the context of the assumptions involved and the limitations are dictated by the restrictions one faces in taking Fourier transforms and applying a Green's function. The Fourier transform that we consider in (13) by nature leads only to the steady-state solution, in which enough time is assumed to have passed so that the effect of the initial conditions are dampened out of the solution. If the transient solution is desired, the Laplace transform has to be employed instead. It is important to note that the inclusion of the initial conditions would make the solutions (26), (27) much more involved. Also, the application of the Fourier transform in time should lead to an explicit transform. For example, Baker et al. [12] mention that in large amplitude vibrations, air damping is due to a drag force proportional to the square of the local beam velocity, |∂w/∂t|(∂w/∂t). The Fourier transform of this ^ nonlinear term in time cannot be explicitly expressed in terms of wðx; ωÞ. Similarly, the application of a Green's function to ^ the equation of motion is limited mathematically by the ability to isolate wðx; ωÞ by applying integration by parts in Eq. (20). This excludes virtually every form of nonlinearity from consideration by this method. Therefore, one of the requirements for applying the FTGF method is for the governing equation to be linear. In cases where the governing equation is nonlinear, such as when damping in large amplitude vibrations is considered, approximation techniques such as perturbation methods and various numerical approaches can be employed.

Acknowledgments The present work is supported by the National Science Foundation under Grant no. CBET #1033117. We would also like to thank Mr. Oleg Goushcha for his assistance in setting up and conducting the experiments.

Appendix A. Solution to Green's function differential equation In Section 4.2, the voltage and beam deflection were found in terms of the still- undetermined Green's function. However, the Green's function differential equation was found to be in (23): ∂4 Gðx; ξ; ωÞ  B4 Gðx; ξ; ωÞ ¼ δðξ  xÞ ∂ξ4 where B  ððmω2  iωca ðωÞÞ=ðYIðωÞ þ iωcs IðωÞÞÞ1=4 .

(A.1)

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

3107

The boundary conditions for the problem are found during the simplifying of (20): Gðx; 0; ωÞ ¼ 0

(A.2)

∂Gðx; 0; ωÞ ¼0 ∂ξ

(A.3)

∂2 Gðx; L; ωÞ ∂Gðx; L; ωÞ ¼ J t ω2 ∂ξ ∂ξ2

(A.4)

∂3 Gðx; L; ωÞ ¼  M t ω2 Gðx; L; ωÞ ∂ξ3

(A.5)

ðYIðωÞ þ iωcs IðωÞÞ

ðYIðωÞ þ iωcs IðωÞÞ

Since the differential Eq. (A.1) contains the Dirac delta function on the right hand side, the solution to the Green's function is not continuous. In fact, the domain can be divided into two segments: 0 rξ rx and x rξrL. Thus, four continuity conditions at ξ ¼x should be introduced along with boundary conditions (A.2)–(A.5) to fully solve (A.1). Three of the continuity conditions involve the continuity of G, ∂G/∂ξ and ∂2G/∂ξ2 at ξ ¼x  and ξ ¼x þ . The fourth continuity condition can be obtained by integrating the original differential equation from ξ ¼x  to ξ¼x þ and simplifying. Thus, the four continuity conditions can be listed as: xþε lim½Gðx; ξ; ωÞξξ ¼ ¼ xε ¼ 0

(A.6)


 ∂Gðx; ξ; ωÞ ξ ¼ x þ ε lim ¼0 ∂ξ ε-0 ξ ¼ xε

(A.7)


2 ξ ¼ x þ ε ∂ Gðx; ξ; ωÞ ¼0 lim ε-0 ∂ξ2 ξ ¼ xε

(A.8)


3 ξ ¼ x þ ε ∂ Gðx; ξ; ωÞ ¼1 lim 3 ε-0 ∂ξ ξ ¼ xε

(A.9)

ε-0

To construct the Green's function, the boundary (A.2)–(A.5) and continuity conditions (A.6)–(A.9) are applied to the general solution of the corresponding homogeneous differential equation of (A.1), in which the right hand side is set to zero. The Green's function is then found to be: 8 ð cos ðBξÞ C 2 sinhðBξÞ þ C 3 sin ðBξÞÞ > > > > > ½ð cos ðBxÞ  coshðBxÞÞ þ ðC 1 þC 2 Þð sin ðBxÞ  sinhðBxÞÞ > > > : 0 r x r ξ rL > >  ðcoshðBξÞ þ C 1 sinhðBξÞ þC 2 sin ðBξÞÞ > > > > > > < ½ð cos ðBxÞ  coshðBxÞÞ þ ðC 3 C 2 Þð sin ðBxÞ  sinhðBxÞÞ (A.10) Gðx; ξ; ωÞ ¼ C 4 ð cos ðBxÞ C 2 sinhðBxÞ þ C 3 sin ðBxÞÞ > > > > > ½ð cos ðBξÞ  coshðBξÞÞ þ ðC 1 þC 2 Þð sin ðBξÞ  sinhðBξÞÞ > > > : 0 r ξ rx r L > > >  ðcoshðBxÞ þ C 1 sinhðBxÞ þC 2 sin ðBxÞÞ > > > > > : ½ð cos ðBξÞ  coshðBξÞÞ þ ðC 3 C 2 Þð sin ðBξÞ  sinhðBξÞÞ where ½B4 ðYIðωÞ þ iωcs IðωÞÞ2  J t M t ω4 ½ cos ðBLÞcoshðBLÞ  sin ðBLÞsinhðBLÞ C1 ¼

 2BðYIðωÞ þ iωcs IðωÞÞω2 ½B2 J t cos ðBLÞsinhðBLÞ þM t sin ðBLÞcoshðBLÞ ½B4 ðYIðωÞ þ iωcs IðωÞÞ2  J t M t ω4 ½ sin ðBLÞcoshðBLÞ  cos ðBLÞsinhðBLÞ

! !

(A.11)

!

(A.12)

þ 2BðYIðωÞ þ iωcs IðωÞÞω2 ½B2 J t cos ðBLÞcoshðBLÞ þ M t sin ðBLÞsinhðBLÞ C2 ¼

B4 ðYIðωÞ þ iωcs IðωÞÞ2 þ J t Mt ω4 4

½B ðYIðωÞ þ iωcs IðωÞÞ2  J t M t ω4 ½ sin ðBLÞcoshðBLÞ  cos ðBLÞsinhðBLÞ þ 2BðYIðωÞ þ iωcs IðωÞÞω2 ½B2 J t cos ðBLÞcoshðBLÞ þ M t sin ðBLÞsinhðBLÞ

 ½B4 ðYIðωÞ þ iωcs IðωÞÞ2  J t Mt ω4 ½ cos ðBLÞcoshðBLÞ þ sin ðBLÞsinhðBLÞ C3 ¼

þ 2BðYIðωÞ þ iωcs IðωÞÞω2 ½B2 J t sin ðBLÞcoshðBLÞ M t cos ðBLÞsinhðBLÞ ½B4 ðYIðωÞ þ iωcs IðωÞÞ2  J t M t ω4 ½ sin ðBLÞcoshðBLÞ  cos ðBLÞsinhðBLÞ þ 2BðYIðωÞ þ iωcs IðωÞÞω2 ½B2 J t cos ðBLÞcoshðBLÞ þ M t sin ðBLÞsinhðBLÞ

!

!

(A.13)

3108

A.H. Danesh-Yazdi et al. / Journal of Sound and Vibration 333 (2014) 3092–3108

C4 ¼

1 2B3 ½C 1 þ2C 2  C 3 

(A.14)

If the end is free, i.e. the tip mass and tip mass moment of inertia are Mt ¼Jt ¼0, the constants (A.11)–(A.13) become: C1 ¼

cos ðBLÞcoshðBLÞ  sin ðBLÞsinhðBLÞ sin ðBLÞcoshðBLÞ  cos ðBLÞsinhðBLÞ

(A.15)

C2 ¼

1 sin ðBLÞcoshðBLÞ  cos ðBLÞsinhðBLÞ

(A.16)

C3 ¼ 

cos ðBLÞcoshðBLÞ þ sin ðBLÞsinhðBLÞ sin ðBLÞcoshðBLÞ  cos ðBLÞsinhðBLÞ

(A.17)

References [1] H.A. Sodano, D.J. Inman, G. Park, Estimation of electric charge output for piezoelectric energy harvesting, Strain Journal 40 (2004) 49–58. [2] A. Erturk, Electromechanical modeling of piezoelectric energy harvester, PhD Thesis, Virginia Polytechnic Institute and State University, 2009. [3] H.D. Akaydin, N. Elvin, Y. Andreopoulos, Energy harvesting from highly unsteady fluid flows using piezoelectric materials, Journal of Intelligent Material Systems and Structures 21 (2010) 1263–1278, http://dx.doi.org/10.1177/1045389X10366317. [4] A. Erturk, D.J. Inman, A distributed parameter electromechanical model for cantilevered piezoelectric energy harvesters, Journal of Vibration and Acoustics ;130http://dx.doi.org/10.1115/1.2890402. [5] C.B. Williams, R.B. Yates, Computer Communications 26 (2003) 1131–1144, http://dx.doi.org/10.1016/S0140-3664(02)00248-7. [6] S. Roundy, P.K. Wright, J.M. Rabaey, A study of low level vibrations as a power source for wireless sensor nodes, Computer Communications 26 (2003) 1131–1144, http://dx.doi.org/10.1016/S0140-3664(02)00248-7. [7] N.E. duToit, B.L. Wardle, Experimental verification of models of models for microfabricated piezoelectric vibration energy harvesters, AIAA Journal 45 (2007) 1126–1137, http://dx.doi.org/10.2514/1.25047. [8] S. Priya, Advances in energy harvesting using low profile piezoelectric transducers, Journal of Electroceramics 19 (2007) 167–184, http://dx.doi.org/ 10.1007/s10832-007-9043-4. [9] S. Adhikari, M.I. Friswell, D.J. Inman, Piezoelectric energy harvesting from broadband random vibrations, Smart Materials and Structures 18 (2009) 115005, http://dx.doi.org/10.1088/0964-1726/18/11/115005. [10] A. Erturk, D.J. Inman, On mechanical modeling of cantilevered piezoelectric vibration energy harvesters, Journal of Intelligent Material Systems and Structures 19 (2008) 1311–1325, http://dx.doi.org/10.1177/1045389X07085639. [11] N.G. Elvin, A.A. Elvin, General equivalent circuit model for piezoelectric generators, Journal of Intelligent Material Systems and Structures 20 (2009) 3–9, http://dx.doi.org/10.1177/1045389X08089957. [12] W.E. Baker, W.E. Woolam, D. Young, Air and internal damping of thin cantilever beams, International Journal of Mechanical Sciences 9 (1967) 743–766. [13] R.F. Gibson, R. Plunkett, A forced-vibration technique for measurement of material damping, Experimental Mechanics 17 (1977) 297–302. [14] J.D. Rogers, L.W. Zachary, K.G. McConnell, Damping characterization of a filled epoxy used for dynamic structural modeling, Experimental Mechanics 26 (1986) 283–291. [15] A.D. Nashif, D.I.G. Jones, J.P. Henderson, Vibration Damping, John Wiley & Sons, New York, 1985. [16] C.C. Chiu, Influence of air and material damping on dynamic elastic modulus measurement, Journal of Materials Science 26 (1991) 4910–4916. [17] H. Shen, H. Ji, J. Qiu, K. Zhu, A semi-passive vibration damping system powered by harvested energy, International Journal of Applied Electromagnetics and Mechanics 31 (2009) 219–233. [18] E. Lefeuvre, A. Badel, L. Petit, C. Richard, D. Guyomar, Semi-passive piezoelectric structural damping by synchronized switching on voltage sources, Journal of Intelligent Material Systems and Structures 17 (2006) 653–660, http://dx.doi.org/10.1177/1045389X06055810. [19] M. Abu-Hilal, Forced vibration of Euler–Bernoulli beams by means of dynamic Green functions, Journal of Sound and Vibration 267 (2003) 191–207. [20] F.B. Hildebrand, Advanced Calculus for Applications, 2nd ed. Prentice-Hall, New Jersey, 1976. [21] J.E. Sader, Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope, Journal of Applied Physics 84 (1998). [22] S. Timoshenko, D. Young, W. Weaver, Vibration Problems in Engineering, Wiley, New York, 1974. [23] S. Adhikari, Damping modelling using generalized proportional damping, Journal of Sound and Vibration 293 (2006) 156–170, http://dx.doi.org/ 10.1016/j.jsv.2005.09.034. [24] A.S. Plouin, E. Balmes, Pseudo-modal representations of large models with viscoelastic behavior, Proceedings of the16th International Modal Analysis Conference on Society for Experimental Mechanics, Inc, 2, 1998, pp. 1440–1446. [25] A. Nechibvute, A.R. Akande, P.V.C. Luhanga, Modelling of a PZT beam for voltage generation, Pertanika Journal of Science and Technology 19 (2011) 259–271.