0
Research Papers

# A Distributed Parameter Electromechanical Model for Cantilevered Piezoelectric Energy HarvestersOPEN ACCESS

[+] Author and Article Information
A. Erturk1

Center for Intelligent Material Systems and Structures, Department of Mechanical Engineering and Department of Engineering Science and Mechanics, Virginia Tech, Blacksburg, VA 24061erturk@vt.edu

D. J. Inman

Center for Intelligent Material Systems and Structures, Department of Mechanical Engineering, Virginia Tech, Blacksburg, VA 24061

Here, as a common practice, the abbreviation PZT is used for a generic piezoelectric ceramic, rather than a specific material.

If the excitation due to the external damping of the medium (which is generally air) is negligible when compared to the inertial excitation (i.e., $Nrc(t)≪Nrm(t)$), one can simply set it equal to zero $(Nrc(t)=0)$.

Yet, one can obtain the value of $c$ for nonzero initial conditions both in the mechanical domain (for the beam) and in the electrical domain (for the electrical circuit). Note that, for nonzero initial conditions in the mechanical domain, the Duhamel integral given by Eq. 33 for the modal mechanical response must be modified accordingly (i.e., initial displacement and velocity terms must be introduced).

Excitation due to air damping is ignored.

1

Corresponding author.

J. Vib. Acoust 130(4), 041002 (Jun 11, 2008) (15 pages) doi:10.1115/1.2890402 History: Received August 16, 2007; Revised December 07, 2007; Published June 11, 2008

## Abstract

Cantilevered beams with piezoceramic layers have been frequently used as piezoelectric vibration energy harvesters in the past five years. The literature includes several single degree-of-freedom models, a few approximate distributed parameter models and even some incorrect approaches for predicting the electromechanical behavior of these harvesters. In this paper, we present the exact analytical solution of a cantilevered piezoelectric energy harvester with Euler–Bernoulli beam assumptions. The excitation of the harvester is assumed to be due to its base motion in the form of translation in the transverse direction with small rotation, and it is not restricted to be harmonic in time. The resulting expressions for the coupled mechanical response and the electrical outputs are then reduced for the particular case of harmonic behavior in time and closed-form exact expressions are obtained. Simple expressions for the coupled mechanical response, voltage, current, and power outputs are also presented for excitations around the modal frequencies. Finally, the model proposed is used in a parametric case study for a unimorph harvester, and important characteristics of the coupled distributed parameter system, such as short circuit and open circuit behaviors, are investigated in detail. Modal electromechanical coupling and dependence of the electrical outputs on the locations of the electrodes are also discussed with examples.

<>

## Introduction

For the past five years, there has been an explosion of research in the area of harvesting energy from ambient vibrations by using the direct piezoelectric effect (1). Research in this area involves understanding the mechanics of vibrating structures, the constitutive behavior of piezoelectric materials, and elementary circuit theory. This promising way of powering small electronic components and remote sensors has attracted researchers from different disciplines of engineering including electrical and mechanical as well as researchers from the field of material science.

The literature includes several experimental demonstrations both in macroscale (2-3) and in microscale (4) piezoelectric energy harvesting. In most of the experimental work, the harvester, which is a cantilevered composite beam with one or more piezoceramic (PZT) layers, is excited harmonically at its fundamental natural frequency for the maximum electrical output. Although most ambient vibration sources do not have harmonic behavior in time, most previous research has assumed harmonic excitation. Furthermore, in some of the works, the excitation frequency of the harvester is so high that the resulting work deviates from the main motivation of the problem since no such frequencies are available in the ambient energy. Nevertheless, especially in the case of microscale harvesters (4), natural frequencies in the kilohertz frequency band are almost inevitable due to the extremely low mass of the structure.

Many researchers have also focused on the mathematical modeling of these harvesters. A reliable mathematical model may allow studying different aspects of energy harvesting, predicting the electrical outputs, and, moreover, optimizing the harvester for the maximum electrical output for a given input. The modeling approaches in the literature include coupled single degree-of-freedom (SDOF) models (5-6), approximate distributed parameter models (with the Rayleigh–Ritz method) (2,6) as well as some distributed parameter modeling approaches (7-8) that consider a single vibration mode and ignore (or oversimplify) the backward coupling in the mechanical domain and even some incorrect modeling attempts (9) based on informal and weak mathematical modeling assumptions.

The SDOF modeling approach (or the so-called lumped parameter modeling) considers the cantilevered beam as a mass-spring-damper system, which is very convenient for coupling the mechanical part of the harvester with a simple electrical harvesting circuit. Although SDOF modeling gives initial insight into the problem by allowing simple closed-form expressions, it is just a simple approximation limited to a single vibration mode of the bender and it lacks several important aspects of the physical system, such as the dynamic mode shape and the accurate strain distribution along the bender. Since these cantilevered harvesters are excited mainly due to the motion of their base, the well-known SDOF harmonic base excitation relation taken from the elementary vibration texts has been frequently referred in the energy harvesting literature both for modeling (6,9) and studying the detailed characteristics (10) of harvesters. Some authors have used the SDOF harmonic base excitation relation just for representing the problem in their experimental work, although they did not provide any detailed mathematical modeling (4,11). It was recently shown (12) that the traditional form of the SDOF harmonic base excitation relation may yield highly inaccurate results both for transverse and longitudinal vibrations of cantilevered beams and bars depending on the tip (proof) mass to beam/bar mass ratio. When a structure is excited due to the motion of its base, the excitation source is nothing but its own inertia. As a result, the conventional effective mass of cantilevered beams and bars results in underestimation of the predicted response especially if there is low or no tip mass. Erturk and Inman (12) introduced correction factors for improving the existing SDOF harmonic base excitation model for both transverse and longitudinal vibrations of cantilevers.

As an alternative modeling approach, Sodano et al. (2) and duToit et al. (6) used the conventional combination of the variational principle (which is also referred to the Hamilton’s principle) and the Rayleigh–Ritz method based on the Euler–Bernoulli beam assumptions. This approximate modeling approach allows predicting the electromechanical response in higher vibration modes. However, it is a numerical approximation technique based on discretization of the continuous distributed parameter system and it is not the exact solution. Other than these approximate techniques, the literature also includes some analytical modeling approaches. In their modeling paper, Lu et al. (7) used a single vibration mode expression (rather than a modal expansion) in the piezoelectric constitutive relation that gives the electric displacement (13) in order to relate the electrical outputs to the mechanical mode shape. Such an approach lacks two important points: backward coupling in the mechanical domain (as coupling in the mechanical equation is not considered) and the contributions from the other vibration modes (as the correct approach implies using the expansion of all vibration modes). Therefore, such a model is approximately valid in the vicinity of that respective vibration mode, and the plots given for a wide frequency band (such as the power plot in Ref. 7) based on this modeling approach are meaningless since all the other vibration modes are missing in the model. Chen et al. (8) also presented a similar distributed parameter modeling approach where the mechanical response was represented by modal expansion, but the effect of backward piezoelectric coupling in the mechanical equation was assumed to be viscous damping. Representing the effect of electromechanical coupling in the mechanical equation by a viscous damping coefficient is a reasonable approach for a certain type of electromagnetic harvesters described by Williams and Yates (14) but not for piezoelectric harvesters. The effect of piezoelectric coupling in the mechanical domain is more sophisticated than just viscous damping as it results in variation of the natural frequencies rather than just attenuation of the motion amplitude, which is discussed in this paper extensively. Quite recently, Ajitsaria et al. (9) proposed a modeling approach for “bimorph piezoelectric cantilever for voltage generation.” Their main reference considers static modeling of piezoelectric actuators where it is reasonable to assume a constant radius of curvature over the beam length. However, Ajitsaria et al. (9) tried to combine these static modeling equations (with constant radius of curvature and static tip force) with the dynamic Euler–Bernoulli beam equation (where the radius of curvature varies) and base excitation (where there is no tip force).

In this paper, we present the exact electromechanical solution of a cantilevered piezoelectric energy harvester for transverse vibrations with Euler–Bernoulli beam assumptions. The harvester beam is assumed to be excited by the motion of its base, which is represented by translation in the transverse direction and small rotation. First, the base motion is not restricted to be harmonic in time so that the general coupled equations are obtained for the mechanical response of the beam and the voltage output. Then, the resulting equations are reduced for the case of harmonic base translation with superimposed harmonic small rotation. Closed-form expressions for the voltage, current, and power outputs as well as the mechanical response are presented. The expressions are further reduced to single-mode relations for the case of excitation around a natural frequency of the bender.

We consider more sophisticated damping mechanisms in mechanical modeling of the harvester. The internal strain rate damping (i.e., Kelvin–Voigt damping) and the external air damping effects are treated more accurately by defining separate damping coefficients. The electrical circuit consists of a resistive load connected to the electrodes bracketing the PZT layer. Therefore, along with the internal capacitance of PZT, the electrical circuit is a first order $(RC)$ circuit. Although the harvester is taken as a unimorph, the formulation for a bimorph can easily be obtained by following a similar procedure.

The analytically obtained electromechanical expressions are used in a parametric case study where certain electromechanical frequency response functions (FRFs) are defined for graphical demonstration. Voltage, current, power, and relative tip motion FRFs are plotted against frequency for a wide range of load resistance, and the resulting trends in the plots are investigated in detail. Although the mathematical model assumes proportional damping, how to relax this assumption is explained briefly along with a short discussion of how to extract the internal strain rate and external air damping coefficients from experimental measurements. Short circuit and open circuit conditions of the FRFs are discussed and the electromechanical outputs are also plotted against load resistance for these extreme conditions of the resistive load. Finally, the importance of modal electromechanical coupling in energy harvesting is highlighted and its dependence on the locations of the electrodes is discussed. The relevance of the modal coupling concept and the locations of the electrodes to the issue of strain nodes of vibration modes, which was previously discussed (12) and experimentally verified (15), is demonstrated here with FRFs.

## Derivation of the Electromechanical Model

We consider the unimorph harvester shown in Fig. 1, which is simply a uniform composite Euler–Bernoulli beam consisting of a PZT layer perfectly bonded to the substructure layer. As well known, another typical harvesting configuration is a bimorph with two PZT layers and that configuration can also be modeled in a similar manner by following the below modeling procedure. The harvester shown in Fig. 1 is connected to the electrical circuit through the electrodes, which bracket the PZT layer. The electrodes are assumed to be perfectly conductive and they cover the entire surface of the PZT at the bottom and at the top (so that the electric field is uniform over the length of the beam). The simple electrical circuit consists of a resistive load only. We assume persistent excitation at the base of the harvester so that continuous electrical outputs can be extracted from the resistive load. In general, the leakage resistance of the PZT is much higher than the load resistance, which allows neglecting it in the electrical circuit as it is normally connected to the circuit in parallel to the resistive load. In Fig. 1, the capacitance of the PZT is considered as internal to the PZT rather than showing it as an external element parallel to the resistive load. As we will see later in this work, the piezoelectric constitutive relations (13) generate the electrical capacitance term. Therefore, although it is not shown in Fig. 1 as an element parallel to the resistive load, the capacitance of the PZT layer is not ignored and it will simply show up in the circuit equation.

The harvester beam is typically excited due to the motion of its base. If the translation and the small rotation of the base are denoted by $g(t)$ and $h(t)$, respectively, then the base motion $wb(x,t)$ on the beam can be represented as (16)Display Formula

$wb(x,t)=g(t)+xh(t)$
(1)
The governing equation of motion can then be written as (12)Display Formula
$∂2M(x,t)∂x2+csI∂5wrel(x,t)∂x4∂t+ca∂wrel(x,t)∂t+m∂2wrel(x,t)∂t2=−m∂2wb(x,t)∂t2−ca∂wb(x,t)∂t$
(2)
where $wrel(x,t)$ is the transverse deflection of the beam relative to its base, $M(x,t)$ is the internal bending moment, $csI$ is the equivalent damping term of the composite cross section due to structural viscoelasticity ($cs$ is the equivalent coefficient of strain rate damping and $I$ is the equivalent area moment of inertia of the composite cross section), $ca$ is the viscous air damping coefficient, and $m$ is the mass per unit length of the beam. Both of the damping mechanisms considered in the model satisfy the proportional damping criterion, hence, they are mathematically convenient for the modal analysis solution procedure (17). The strain rate damping indeed shows itself as an internal moment $(Ms=csI∂3wrel∕∂x2∂t)$ in the resulting equation of motion. For convenience, we consider it directly outside the internal moment term in Eq. 2; hence, it appears as a separate term.

The internal moment can be obtained by integrating the first moment of the stress distribution at a cross section over the cross-sectional area. The piezoelectric constitutive relations (13) give the stress-strain (and electric field) relations and they are expressed for the substructure and the PZT layers asDisplay Formula

$T1s=YsS1s$
(3)
Display Formula
$T1p=Yp(S1p−d31E3)$
(4)
respectively. Here, $T$ is the stress, $S$ is the strain, $Y$ is Young’s modulus; $d$ is the piezoelectric constant, and $E$ is the electric field. Equation 4 is obtained from the piezoelectric constitutive relation $S1=s11ET1+d31E3$, where $s11E$ is the elastic compliance at constant electric field and therefore $Yp$ is simply the reciprocal of $s11E$. Furthermore, subscript/superscript $p$ and $s$ stand for PZT and substructure layers, respectively; 1 and 3 directions are coincident with $x$ and $y$ directions, respectively (where 1 is the direction of axial strain and 3 is the direction of polarization). Then, the internal moment can be written asDisplay Formula
$M(x,t)=−∫hahbT1sbydy−∫hbhcT1pbydy$
(5)
where $b$ is the width of the beam, $ha$ is the position of the bottom of the substructure layer from the neutral axis, $hb$ is the position of the bottom of the PZT layer (therefore, top of the substructure layer) from the neutral axis, and $hc$ is the position of the top of the PZT layer from the neutral axis (see the Appendix). Expressing the bending strain in terms of radius of curvature (18) and employing Eqs. 3,4 in Eq. 5 giveDisplay Formula
$M(x,t)=∫hahbYsb∂2wrel(x,t)∂x2y2dy+∫hbhcYpb∂2wrel(x,t)∂x2y2dy−∫hbhcv(t)Ypbd31hpydy$
(6)
where the uniform electric field is written in terms of the voltage $v(t)$ across the PZT and the thickness $hp$ of the PZT $(E3(t)=−v(t)∕hp)$. Equation 6 can be reduced toDisplay Formula
$M(x,t)=YI∂2wrel(x,t)∂x2+ϑv(t)$
(7)
where $YI$ is the bending stiffness of the composite cross section given byDisplay Formula
$YI=b[Ys(hb3−ha3)+Yp(hc3−hb3)3]$
(8)
and the coupling term $ϑ$ can be written asDisplay Formula
$ϑ=−Ypd31b2hp(hc2−hb2)$
(9)
If the PZT layer and/or the electrodes do not cover the entire length of the beam but the region $x1⩽x⩽x2$, then the second term in Eq. 7 should be multiplied by $H(x−x1)−H(x−x2)$, where $H(x)$ is the Heaviside function. Note that, in energy harvesting from higher vibration modes, it becomes necessary to use segmented electrode pairs in order to avoid cancellation of the charge collected by continuous electrode pairs (15). In such a case, one needs to define separate voltage terms, which should appear in Eq. 7 as separate terms multiplied by Heaviside functions. However, in our case, we assume that the electrodes and the PZT layer cover the entire length of the beam displayed in Fig. 1 and it is convenient to rewrite Eq. 7 asDisplay Formula
$M(x,t)=YI∂2wrel(x,t)∂x2+ϑv(t)[H(x)−H(x−L)]$
(10)
where $L$ is the length of the beam. In Eq. 10, although the PZT layer and the electrodes cover the entire beam length, Heaviside functions are associated with the second term in Eq. 10 in order to ensure the survival of this term when the internal moment expression $M(x,t)$ is used in the differential equation of motion given by Eq. 2. Then, employing Eq. 10 in Eq. 2 yieldsDisplay Formula
$YI∂4wrel(x,t)∂x4+csI∂5wrel(x,t)∂x4∂t+ca∂wrel(x,t)∂t+m∂2wrel(x,t)∂t2+ϑv(t)[dδ(x)dx−dδ(x−L)dx]=−m∂2wb(x,t)∂t2−ca∂wb(x,t)∂t$
(11)
where $δ(x)$ is the Dirac delta function and it satisfies (19)Display Formula
$∫−∞∞d(n)δ(x−x0)dx(n)f(x)dx=(−1)ndf(n)(x0)dx(n)$
(12)

Equation 11 is the mechanical equation of motion with electrical coupling. In order to obtain the electrical circuit equation with mechanical coupling, one should consider the following piezoelectric constitutive relation:Display Formula

$D3=d31T1+ε33TE3$
(13)
where $D3$ is the electric displacement and $ε33T$ is the permittivity at constant stress. If we rearrange Eq. 13 to express axial stress $T1$ in terms of bending strain $S1$ and Young’s modulus of PZT $(Yp=1∕s11E)$, permittivity component must be replaced by permittivity at constant strain, which is given by $ε33S=ε33T−d312Yp$(13). After writing the electric field in the PZT in terms of the voltage across the PZT (that is, $E3(t)=−v(t)∕hp$), one obtainsDisplay Formula
$D3(x,t)=d31YpS1(x,t)−ε33Sv(t)hp$
(14)
The average bending strain at position $x$ and time $t$ in the PZT layer can be expressed as a function of the distance $hpc$ of the center of the PZT layer (in thickness direction) to the neutral axis (see the Appendix) and curvature of the beam at position $x$ and time $t$,Display Formula
$S1(x,t)=−hpc∂2wrel(x,t)∂x2$
(15)
Therefore, Eq. 14 becomesDisplay Formula
$D3(x,t)=−d31Yphpc∂2wrel(x,t)∂x2−ε33Sv(t)hp$
(16)

The electric charge $q(t)$ developed in the PZT (and collected by the electrodes) can be obtained by integrating the electric displacement over the electrode area as (13)Display Formula

$q(t)=∫AD⋅ndA=−∫x=0L(d31Yphpcb∂2wrel(x,t)∂x2+ε33Sbv(t)hp)dx$
(17)
where $D$ is the vector of electric displacements and $n̰$ is the unit outward normal. Clearly, the nonzero terms of these vectors are the ones in 3 direction (i.e., in the $y$-direction). Then, the current generated by the PZT can be given byDisplay Formula
$i(t)=dq(t)dt=−∫x=0Ld31Yphpcb∂3wrel(x,t)∂x2∂tdx−ε33SbLhpdv(t)dt$
(18)
Here, the current generated is a function with two components: The first component is due to the vibratory motion of the beam and the second component includes the voltage across the PZT. It is obvious from Eq. 18 that the second term is due to the static capacitance of the PZT. In the literature, the term $ε33SbL∕hp$ is called the capacitance of the PZT layer (in general, it is denoted by $Cp$) and this capacitance is connected to the resistive load $(Rl)$ shown in Fig. 1 in parallel although, physically, the capacitance of the PZT is internal to the PZT itself. Since the current expression given by Eq. 18 includes the capacitance information of the PZT, in this model, it is convenient to connect the PZT directly to the resistive load as a current source without any external capacitive element as we did in Fig. 1. Then, the voltage across the resistive load is simplyDisplay Formula
$v(t)=Rli(t)=−Rl[∫x=0Ld31Yphpcb∂3wrel(x,t)∂x2∂tdx+ε33SbLhpdv(t)dt]$
(19)
or alternately the electrical circuit equation can be written asDisplay Formula
$ε33SbLhpdv(t)dt+v(t)Rl=−∫x=0Ld31Yphpcb∂3wrel(x,t)∂x2∂tdx$
(20)
where $v(t)$ is the voltage across the resistive load.

Equations 11,20 are the distributed parameter electromechanical equations for a cantilevered piezoelectric energy harvester in transverse vibrations.

## General Transient Base Motions

The aim of this section is to provide the formal solution procedure of the electromechanically coupled equations of a harvester beam, which are Eqs. 11,20. The relative vibratory motion of the beam can be represented by an absolutely and uniformly convergent series of the eigenfunctions asDisplay Formula

$wrel(x,t)=∑r=1∞ϕr(x)ηr(t)$
(21)
where $ϕr(x)$ and $ηr(t)$ are the mass normalized eigenfunction and the modal coordinate of the clamped-free beam for the $rth$ mode, respectively. Since the system is proportionally damped, the eigenfunctions denoted by $ϕr(x)$ are indeed the mass normalized eigenfunctions of the corresponding undamped free vibration problem (17) given byDisplay Formula
$ϕr(x)=1mL[coshλrLx−cosλrLx−σr(sinhλrLx−sinλrLx)]$
(22)
where the $λr$’s are the dimensionless frequency numbers obtained from the characteristic equation given byDisplay Formula
$1+cosλcoshλ=0$
(23)
and $σr$ is expressed asDisplay Formula
$σr=sinhλr−sinλrcoshλr+cosλr$
(24)
Note that, for the sake of completeness, the tip mass (or the so-called proof mass) is excluded and the above expressions are for a cantilevered beam without a tip mass. The effect of a tip mass on the resulting formulation of the base excitation problem as well as the respective expressions for the eigenfunctions and the characteristic equation can be found in a recent paper by Erturk and Inman (20).

The mass normalized form of the eigenfunctions given by Eq. 22 satisfies the following orthogonality conditions:Display Formula

$∫x=0Lmϕs(x)ϕr(x)dx=δrs,∫x=0LYIϕs(x)d4ϕr(x)dx4dx=ωr2δrs$
(25)
where $δrs$ is the Kronecker delta, defined as being equal to unity for $s=r$ and equal to zero for $s≠r$, and $ωr$ is the undamped natural frequency of the $rth$ mode given byDisplay Formula
$ωr=λr2YImL4$
(26)

Using Eq. 21 in the partial differential equation of motion along with the orthogonality conditions given by Eq. 25, the electromechanically coupled ordinary differential equation for the modal response of the beam can be obtained asDisplay Formula

$d2ηr(t)dt2+2ζrωrdηr(t)dt+ωr2ηr(t)+χrv(t)=Nr(t)$
(27)
whereDisplay Formula
$χr=ϑ∣dϕr(x)dx∣x=L$
(28)
is the modal coupling term andDisplay Formula
$ζr=csIωr2YI+ca2mωr$
(29)
is the mechanical damping ratio that includes the effects of both strain rate damping and viscous air damping (12). It is clear from Eq. 29 that strain rate damping is assumed to be proportional to the bending stiffness of the beam, whereas air damping is assumed to be proportional to the mass per unit length of the beam.

The modal mechanical forcing function $Nr(t)$ can be expressed asDisplay Formula

$Nr(t)=Nrm(t)+Nrc(t)$
(30)
Here, the components of mechanical excitation (which are the inertial and the damping excitation terms) are given by the following expressions, respectively;
$Nrm(t)=−m(γrwd2g(t)dt2+γrθd2h(t)dt2),$
$Nrc(t)=−ca(γrwdg(t)dt+γrθdh(t)dt)$
whereDisplay Formula
$γrw=∫x=0Lϕr(x)dx,γrθ=∫x=0Lxϕr(x)dx$
(32)

The modal response (that is, the solution of the ordinary differential equation given by Eq. 27) can be expressed using the Duhamel integral asDisplay Formula

$ηr(t)=1ωrd∫τ=0t[Nr(τ)−χrv(τ)]e−ζrωr(t−τ)sinωrd(t−τ)dτ$
(33)
where $ωrd=ωr1−ζr2$ is the damped natural frequency of the $rth$ mode. Clearly, obtaining $ηr(t)$ from Eq. 33 implies knowing the voltage $v(t)$ across the PZT.

Consider Eq. 20, which is a first order ordinary differential equation for $v(t)$, whose forcing term is a function of the relative vibratory motion. Using Eq. 21 in Eq. 20 yieldsDisplay Formula

$dv(t)dt+hpRlε33SbLv(t)=∑r=1∞φrdηr(t)dt$
(34)
whereDisplay Formula
$φr=−d31Yphpchpε33SL∫x=0Ld2ϕr(x)dx2dx=−d31Yphpchpε33SL∣dϕr(x)dx∣x=L$
(35)
Equation 34 can be solved for $v(t)$ by using the following integrating factor:Display Formula
$Ψ(t)=et∕τc$
(36)
where $τc$ is the time constant of the circuit given byDisplay Formula
$τc=Rlε33SbLhp$
(37)
Multiplying both sides of Eq. 34 by the integration factor and solving the resulting equation yield the following expression for the voltage across the resistive load:Display Formula
$v(t)=e−t∕τc(∫et∕τc∑r=1∞φrdηr(t)dtdt+c)$
(38)
where $c$ is an arbitrary constant that depends on the initial value of the voltage $v(t)$ across the resistive load as well as the initial velocity of the beam. However, the form of Eq. 33 assumes the initial displacement and the velocity of the beam to be zero (a more complicated form of it could handle the initial displacement and velocity terms). Therefore, since we assume zero initial conditions for the mechanical domain, the term $c$ in Eq. 38 depends only on the initial value of the voltage across the resistive load. Hence, if the initial voltage is zero $(v(0)=0)$, simply, $c=0$ and this is what we assume for the sake of simplicity.

Then, if the summation sign and the integral are switched in Eq. 38, and for $c=0$, one obtainsDisplay Formula

$v(t)=e−t∕τc∑r=1∞φr∫et∕τcdηr(t)dtdt$
(39)

Now, consider the resulting coupled equations (in time domain), which are given by Eqs. 33,39. One can either eliminate the voltage term $v(t)$ or the modal mechanical response term $ηr(t)$ in order to obtain a single equation in $v(t)$ or $ηr(t)$. Eliminating $ηr(t)$ and obtaining a single expression for $v(t)$ is preferable for the energy harvesting problem. When the modal mechanical response is eliminated in Eq. 39, one obtains the following implicit expression for the voltage across the resistive load:Display Formula

$v(t)=e−t∕τc∑r=1∞φrωrd∫et∕τcddt(∫τ=0t[Nr(τ)−χrv(τ)]e−ζrωr(t−τ)sinωrd(t−τ)dτ)dt$
(40)
When the voltage term is eliminated in Eq. 33, one obtains the following implicit equation for the modal response:Display Formula
$ηr(t)=1ωrd∫τ=0t[Nr(τ)−χre−τ∕τc∑r=1∞φr∫eτ∕τcdηr(τ)dτdτ]e−ζrωr(t−τ)sinωrd(t−τ)dτ$
(41)
The solution of Eq. 41 for $ηr(t)$ can be used in Eq. 21 to obtain the coupled physical response of the beam.

## Harmonic Base Motion (Translation With Small Rotation) at an Arbitrary Frequency

In this section, we assume the translation and small rotation of the beam to be harmonic in time (i.e., $g(t)=Y0ejωt$ and $h(t)=θ0ejωt$, where $Y0$ and $θ0$ are the amplitudes of the base translation and rotation, respectively, $ω$ is the driving frequency, and $j=−1$ is the unit imaginary number) and obtain closed-form expressions for steady state voltage output and beam response. Since the electromechanical system is assumed to be linear, we expect the voltage output to be harmonic also in the form of $v(t)=V0ejωt$ (where $V0$ is the amplitude of the harmonic voltage across the resistive load). Then, the modal equation of motion given by Eq. 27 can be reduced toDisplay Formula

$d2ηr(t)dt2+2ζrωrdηr(t)dt+ωr2ηr(t)+χrV0ejωt=mω2(γrwY0+γrθθ0)ejωt$
(42)
The steady state solution of Eq. 42 can be expressed asDisplay Formula
$ηr(t)=[mω2(γrwY0+γrθθ0)−χrV0]ejωtωr2−ω2+j2ζrωrω$
(43)
It is important to note that the input in the energy harvesting problem is the base motion $wb(x,t)=(Y0+xθ0)ejωt$ and the output is the voltage across the resistive load $v(t)=V0ejωt$. Although they look like two separate forcing terms in Eq. 43, $v(t)$ is indeed induced due to $wb(x,t)$, but it acts back as a forcing term in the mechanical equation due to the electromechanical coupling. Furthermore, $v(t)$ and $wb(x,t)$ do not have to be in phase, which makes $V0$ complex valued although $Y0$ and $θ0$ are real valued.

For the electrical circuit, from Eq. 34, one obtainsDisplay Formula

$(1+jωτcτc)V0ejωt=∑r=1∞φrdηr(t)dt$
(44)
Then, using Eq. 43 in Eq. 44, one can eliminate $ηr(t)$ and obtain the following implicit expression for the voltage amplitude $V0$ across the load resistance:Display Formula
$(1+jωτcτc)V0=∑r=1∞φrjω[mω2(γrwY0+γrθθ0)−χrV0]ωr2−ω2+j2ζrωrω$
(45)
However, it is now straightforward to express $V0$ explicitly asDisplay Formula
$V0=∑r=1∞jmω3φr(γrwY0+γrθθ0)ωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc$
(46)
which is the amplitude of the voltage across the resistive load. Therefore, the voltage response $v(t)$ (across the resistive load) due to the harmonic base motion $wb(x,t)=(Y0+xθ0)ejωt$ can be expressed asDisplay Formula
$v(t)=∑r=1∞jmω3φr(γrwY0+γrθθ0)ejωtωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc$
(47)
Note that, after obtaining the time history of the voltage across the resistive load, one can easily obtain the current generated by the PZT by using $i(t)=v(t)∕Rl$. Furthermore, the instantaneous power output can be expressed using the well-known relation $P(t)=v2(t)∕Rl$. It is worthwhile to mention that, due to the similarity between the energy harvesting problem (for persistent vibrations) and the piezoelectric shunt damping problem, some of the equations derived here may also be applicable to the multimode shunt damping of a cantilevered beam (see, for instance, Moheimani et al. (21)). Note that the resistive load connected to the leads of the electrodes in the circuit can easily be replaced by the impedance of a more general (but linear) $RLC$ circuit and we can obtain closed-form expressions for harmonic excitations as long as the linearity of the electromechanical system is preserved.

It is also possible to eliminate the voltage term in Eq. 43 in order to obtain the modal mechanical response of the beam asDisplay Formula

$ηr(t)=[(γrwY0+γrθθ0)−χr(∑r=1∞jωφr(γrwY0+γrθθ0)ωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc)]mω2ejωtωr2−ω2+j2ζrωrω$
(48)
which can be employed in Eq. 21 to obtain the vibratory re-sponse of the beam relative to its base asDisplay Formula
$wrel(x,t)=∑r=1∞[(γrwY0+γrθθ0)−χr(∑r=1∞jωφr(γrwY0+γrθθ0)ωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc)]mω2ϕr(x)ejωtωr2−ω2+j2ζrωrω$
(49)

## Harmonic Base Translation at an Arbitrary Frequency

In the energy harvesting literature, the base excitation is considered as harmonic translation in the transverse direction and the beam is assumed to be not rotating. In such a case (i.e., when $h(t)=0$), the resulting expressions for the steady state voltage response can be reduced toDisplay Formula

$v(t)=∑r=1∞jmω3φrγrwωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτcY0ejωt$
(50)
and the vibratory response of the beam relative to its base can be obtained fromDisplay Formula
$wrel(x,t)=∑r=1∞[γrw−χr(∑r=1∞jωφrγrwωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc)]mω2ϕr(x)Y0ejωtωr2−ω2+j2ζrωrω$
(51)

## Harmonic Base Translation Around a Natural Frequency

If the beam is excited around the natural frequency of the $rth$ mode, the main contributions in the summation signs appearing in Eqs. 50,51 are from the $rth$ mode. In most cases, the mode of interest is the fundamental vibration mode of the harvester (i.e., $r=1$). Therefore, it is a useful practice to consider the beam to be excited around $ω1$. The reduced expression for the voltage across the resistive load can be written asDisplay Formula

$v*(t)=jτcmω3φ1γ1wjωτcχ1φ1+(1+jωτc)(ω12−ω2+j2ζ1ω1ω)Y0ejωt$
(52)
and the reduced expression for vibratory response of the beam relative to its base isDisplay Formula
$wrel*(x,t)=(1+jωτc)mω2γ1wϕ1(x)jωτcχ1φ1+(1+jωτc)(ω12−ω2+j2ζ1ω1ω)Y0ejωt$
(53)
where the superscript “$*$” denotes that the respective expression is reduced for excitation around a specific vibration mode (which is the fundamental mode in this case).

Equation 52 can be rewritten asDisplay Formula

$v*(t)=∣V0*∣ej(ωt+Φv)$
(54)
where the amplitude of the voltage output is (assuming $Y0>0$)Display Formula
$∣V0*∣=τcmω3Y0∣φ1γ1w∣[ω12−ω2(1+2τcζ1ω1)]2+[2ζ1ω1ω+τcω(χ1φ1+ω12−ω2)]2$
(55)
and the phase angle between the base displacement and the reduced voltage output is simplyDisplay Formula
$Φv=π2sgn(φ1γ1w)−tan−1(2ζ1ω1ω+τcω(χ1φ1+ω12−ω2)ω12−ω2(1+2τcζ1ω1))$
(56)
where “sgn” is the signum function.

Similarly, Eq. 53 can be expressed asDisplay Formula

$wrel*(x,t)=∣Wrel*(x)∣ej(ωt+Φw)$
(57)
where the amplitude of the vibratory motion at point $x$ (relative to the base of the beam) can be written asDisplay Formula
$∣Wrel*(x)∣=mω2Y0∣γ1wϕ1(x)∣1+(ωτc)2[ω12−ω2(1+2τcζ1ω1)]2+[2ζ1ω1ω+τcω(χ1φ1+ω12−ω2)]2$
(58)
and the phase angle between the base displacement and the relative displacement at point $x$ is simplyDisplay Formula
$Φw=tan−1(ωτcγ1wϕ1(x)γ1wϕ1(x))−tan−1(2ζ1ω1ω+τcω(χ1φ1+ω12−ω2)ω12−ω2(1+2τcζ1ω1))$
(59)

The reduced expression for the current flow through the resistive load can then be obtained asDisplay Formula

$i*(t)=∣I0*∣ej(ωt+Φi)$
(60)
whereDisplay Formula
$∣I0*∣=Cpmω3Y0∣φ1γ1w∣[ω12−ω2(1+2τcζ1ω1)]2+[2ζ1ω1ω+τcω(χ1φ1+ω12−ω2)]2$
(61)
and $Φi=Φv$, i.e., the voltage across the resistive load and the current flow through it are in phase, reasonably. In Eq. 61, $Cp$ is the capacitance of the PZT due to $Cp=ε33SbL∕hp$ as mentioned previously and $Cp=τc∕Rl$ from Eq. 37.

The amplitude of the power output (reduced for the fundamental vibration mode) can then be expressed asDisplay Formula

$∣P0*∣=Cpτc(mφ1γ1wω3Y0)2[ω12−ω2(1+2τcζ1ω1)]2+[2ζ1ω1ω+τcω(χ1φ1+ω12−ω2)]2$
(62)

It is important to note that Eqs. 52,53,54,55,56,57,58,59,60,61,62 are valid only at the excitation frequencies around the first natural frequency of the harvester beam. Replacing the subscript “1” by $r$ in Eqs. 52,53,54,55,56,57,58,59,60,61,62, one can obtain the respective electrical and mechanical expressions for excitation frequencies around the $rth$ natural frequency.

## Parametric Case Study for a Unimorph Harvester

In this section, we analyze the unimorph harvester shown in Fig. 2 by using the proposed distributed parameter model. The geometric, material, and the electromechanical parameters of the harvester are given in Table 1. The conductive electrodes bracketing the PZT layer (not shown in Fig. 2) are assumed to be covering the entire length of the harvester beam. The excitation of the harvester is due to the harmonic translation $wb(t)=Y0ejωt$ of its base in the transverse direction and we are interested in the steady state dynamic behavior of the system. Rather than specifying certain $Y0$ and $ω$ values for the input, it is preferred to obtain the results in terms of these parameters so that it becomes possible to represent the electromechanical outputs as FRFs.

Let the frequency range of interest be $0–1000Hz$. Then, for the parameters given in Table 1, it is straightforward to show that the uncoupled harvester has three natural frequencies lying in this frequency range. The first part of our analysis focuses on the following four important FRFs: voltage across the resistive load, current flow through the resistive load, electrical power output, and the relative motion transmissibility from the base to the tip of the beam. The first three FRFs are for the electrical domain, and clearly, the third one (power output) depends on the first two (voltage and current). The last FRF gives the ratio of the vibration amplitude at the tip of the beam (relative to its base) to the amplitude of the input base translation. Therefore, only the mechanical FRF gives some idea about the electromechanical effects on the beam caused by the energy harvesting circuit at different vibration modes. It is worthwhile to mention that distributed parameter modeling is not limited to the motion at the tip of the beam as the analysis presented here allows predicting the coupled vibratory response at any point along the beam.

Due to the electromechanical coupling, each vibration mode has a short circuit resonance frequency $ωrsc$ (for $Rl→0$) and an open circuit resonance frequency $ωroc$ (for $Rl→∞$), where subscript $r$ stands for the mode number. These frequencies are defined for the extreme cases of the resistive load, and for its moderate values, reasonably, the resonance frequency of the respective vibration mode takes values between $ωrsc$ and $ωroc$. In the below analysis, after plotting each FRF against the excitation frequency, we investigate the short circuit and open circuit behaviors of the outputs by plotting them against the load resistance.

###### Identification of Mechanical Damping Coefficients

Before presenting the resulting FRFs and discussing their respective trends, it may be appropriate to add a few words on mechanical damping and evaluation of the mechanical damping coefficients from experimental measurements. With the form of the differential equation of motion given by Eq. 2, we assumed two separate damping terms for the internal (strain rate) and the external viscous (air) damping mechanisms, which are $csI$ and $ca$, respectively. This is the formal way of treating the problem and a detailed discussion is available in Ref. 12. As mentioned previously, $csI$ is assumed to be stiffness proportional and $ca$ is assumed to be mass proportional. Evaluating these damping coefficients from experimental measurements requires knowing the frequencies and damping ratios of two separate vibration modes (say, modes $m$ and $n$) (22). Then, for $r=m$ and $r=n$, one has two linear algebraic equations to solve for the unknowns $csI$ and $ca$ (see Eq. 29). Suppose that, for the harvester of our interest (Fig. 2), we obtain the damping ratios of the first two modes as $ζ1=0.010$ and $ζ2=0.013$ from experimental modal analysis under short circuit conditions $(Rl→0)$ and this is what we assume for our parametric study. We also know that $ω1=300.4rad∕s$ and $ω2=1882.5rad∕s$. Then, the two equations coming from Eq. 29 give the proportionality constants as $csI∕YI=1.2433×10−5s∕rad$ and $ca∕m=4.886rad∕s$. Once these proportionality constants are used in the mathematical model, the rest of the modal damping ratios are automatically set to the following numbers: $ζ3=0.033$, $ζ4=0.064$, $ζ5=0.106$, and so on. However, one should always keep in mind that proportional damping is a convenient mathematical modeling assumption and the physical system may not agree with this assumption. In other words, the damping ratios of higher vibration modes may not converge to the above values. At this point, it is important to note that instead of using $csI$ and $ca$ in the distributed parameter model we propose, one can always use the modal damping ratios ($ζr$’s) obtained experimentally directly in the modal expansions appearing in the resulting electromechanical expressions. For example, if the first three modes are of interest, one obtains the damping ratios of these modes experimentally and employs them in the summation terms of the resulting equations by just considering these three modes $(n=3)$. This way allows a relaxation in the proportional damping assumption and one does not need to obtain the values of $csI$ and $ca$.

###### Frequency Response of Voltage Output

Having discussed some important points regarding the mechanical damping and its accurate evaluation, we start investigating the resulting electromechanical FRFs. Since the value of load resistance $Rl$ is an important parameter that shapes the dynamic behavior of the system, we plot the FRFs for five different orders of magnitude of $Rl$ ranging from $102Ωto106Ω$. As discussed previously, short circuit behavior is expected for low values of load resistance $(Rl→0)$, whereas the system is expected to move toward the open circuit conditions for large values of load resistance $(Rl→∞)$.

The voltage FRF is described as the ratio of the voltage output to the base acceleration. Therefore, this FRF can be extracted from Eq. 50 asDisplay Formula

$v(t)−ω2Y0ejωt=∑r=1∞−jmωφrγrwωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc$
(63)

The modulus of the “voltage FRF” given by Eq. 63 is plotted in Fig. 3. As can be seen from Fig. 3, the amplitude of the voltage output increases monotonically with increasing load resistance for all excitation frequencies. Furthermore, with increasing load resistance, the resonance frequency of each vibration mode moves from the short circuit resonance frequency $(ωrsc)$ to the open circuit resonance frequency $(ωroc)$. As an example, the enlarged view of the first resonance is displayed in Fig. 3. For excitation at the first vibration mode, the maximum voltage is obtained for $ω1sc=47.8Hz$ when $Rl=102Ω$. However, when $Rl=106Ω$ is used, the resonance frequency of the fundamental vibration mode becomes $ω1oc=48.8Hz$. For every excitation frequency, the maximum voltage output is obtained when the system is close to open circuit conditions. A similar trend (the existence of open circuit and short circuit resonance frequencies) is valid for all vibration modes. Table 2 shows the frequencies of the short circuit and the open circuit resonances for the first three vibration modes.

As mentioned, the voltage output increases monotonically with increasing load resistance at every excitation frequency. The two important excitation frequencies of the fundamental vibration mode are the short circuit and the open circuit resonance frequencies, which are $ω1sc=47.8Hz$ and $ω1oc=48.8Hz$, respectively. Variation of the voltage output with load resistance for excitations at these frequencies is shown in Fig. 4. As can be seen from Fig. 4, for low values of load resistance, the voltage outputs at these two particular excitation frequencies increase with the same slope (in log-log scale) and the voltage output at the short circuit resonance frequency is higher since the system is close to short circuit conditions. However, the curves intersect at a certain value of load resistance (around $39.8kΩ$) and for the values of load resistance higher than the value at the intersection point, the voltage output at the open circuit resonance frequency is higher expectedly. The voltage output becomes less sensitive to the variations in the load resistance at open circuit conditions (i.e., for very large values of load resistance).

###### Frequency Response of Current Output

The current FRF can easily be obtained by dividing the voltage FRF by the load resistance as follows:Display Formula

$i(t)−ω2Y0ejωt=v(t)−Rlω2Y0ejωt=∑r=1∞−jmωφrγrwωr2−ω2+j2ζrωrωRl(∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc)$
(64)

The modulus of the current FRF is plotted against frequency in Fig. 5. Unlike the voltage FRF shown in Fig. 3, the amplitude of the current decreases with increasing load resistance. Indeed, this is the opposite of the voltage behavior shown in Fig. 3 but the behavior is still monotonic. For every excitation frequency, the maximum value of the current is obtained when the system is close to short circuit conditions.

Figure 6 shows the current output as a function of load resistance for excitations at the short circuit and the open circuit resonance frequencies of the first mode. It is clear from Fig. 6 that the current is highly insensitive to the variations of the load resistance at the range of its low values (i.e., the slope is almost zero). In this relatively low load resistance region, the current output is higher at the short circuit resonance frequency, as in the case of voltage (Fig. 4), since the system is close to short circuit conditions. Then, current starts decreasing with further increase in load resistance, and at a certain value of load resistance (again, around $39.8kΩ$), the curves intersect. For the values of load resistance higher than the value at this intersection point, the current output at the open circuit resonance frequency becomes higher since the system approaches the open circuit conditions.

###### Frequency Response of Power Output

The FRF of the power output is simply the product of the voltage and the current FRFs given by Eqs. 63,64, respectively. Therefore, unlike the voltage and the current FRFs, the power output FRF is defined as the power divided by the square of the base acceleration. The modulus of the power output FRF is displayed in Fig. 7. Since it is the product of two FRFs that have the opposite behaviors with increasing load resistance, the behavior of the power output FRF with load resistance is more interesting than the previous two electrical outputs and it deserves more discussion. It is clear from Fig. 7 that the power output FRF does not exhibit a monotonic behavior with increasing (or decreasing) load resistance. Among the sample values of the load resistance considered in this work, the value of maximum power output for the first vibration mode corresponds to $Rl=105Ω$ (see the first enlarged view in Fig. 7) at frequency $ω=48.66Hz$, which is expectedly around the open circuit resonance frequency of the first mode (Table 2) for this relatively large value of load resistance. Considering the second vibration mode (see the second enlarged view in Fig. 7), one observes that the maximum power output is obtained for $Rl=104Ω$ at frequency $ω=300.88Hz$ and this frequency has a moderate value between the short circuit and the open circuit resonance frequencies of the second vibration mode (Table 2) since the respective resistive load also has a moderate value. According to Fig. 7, among the sample values of load resistance employed in the analysis, the maximum power output for excitation at the third vibration mode corresponds to $Rl=103Ω$ at frequency $ω=838.34Hz$, which is close to the short circuit resonance frequency of the third vibration mode (Table 2) as this resistive load has a relatively low value. One should note that the values of the load resistance we use in this analysis are taken arbitrarily to observe the general trends. Therefore, the maximum power outputs obtained from each vibration mode are for these sample values and they are not necessarily the maximum possible (or the optimized) power outputs. It is a straightforward practice to obtain the optimum resistive load and its respective resonance frequency for each vibration mode and it is beyond the discussion of this section, which aims to address more general points. Another interesting point to mention is the switching between the curves of different values of load resistance, which results in intersections between the FRFs. These intersections are observed not only around the resonance frequencies (see the first enlarged view, e.g., the curves for $Rl=104Ω$ and $Rl=105Ω$ intersect at $48.19Hz$) but also they are observed at the off-resonance frequencies (e.g., the curves for $Rl=103Ω$ and $Rl=104Ω$ intersect at $193.68Hz$). At these intersection frequencies, the two respective load resistance values yield the same power output.

We further investigate the variation of power output with load resistance for excitations at the short circuit and open circuit resonance frequencies of the first vibration mode through Fig. 8. It can be remembered from Figs.  46 that the voltage and the current outputs obtained at the short circuit resonance frequency are higher than the ones obtained at the open circuit resonance frequency up to a certain load resistance (approximately $39.8kΩ$ in this case) after which the opposite is valid. Since the power output is simply the product of the voltage and current, this observation is also valid for the power versus load resistance curves. As can be seen form Fig. 8, we have the same intersection point $(Rl≅39.8kΩ)$ and the power output at the short circuit resonance frequency is higher before this point, whereas the power output at the open circuit resonance frequency is higher after this point. The trend in the low load resistance region is similar to that of the voltage output where the graphs increase with the same slope (in log-log scale) with increasing load resistance.

More importantly, since the behavior of power with changing load resistance is not monotonic, both of the power graphs shown in Fig. 8 display peak values, which correspond to the optimum values of load resistance. When the optimum values of load resistance are used for each of the cases (short circuit and open circuit excitations), both of them yield the same power output. Considering Figs.  46, it can be observed that neither voltages nor currents are identical at these optimum values of load resistance for excitations at short circuit and open circuit resonance frequencies. However, the products of voltage and current for both cases are identical so that the power outputs for these resistive loads are identical for the short circuit and open circuit resonance frequency excitations separately.

###### Frequency Response of Beam Vibration

Now, we define the relative tip motion FRF (or the relative motion transmissibility function), which is the ratio of the vibration (displacement) amplitude at the tip of the beam (relative to the base) to the amplitude of the base displacement. Therefore, this mechanical FRF can be expressed using Eq. 51 asDisplay Formula

$wrel(L,t)Y0ejωt=∑r=1∞[γrw−χr(∑r=1∞jωφrγrwωr2−ω2+j2ζrωrω∑r=1∞jωχrφrωr2−ω2+j2ζrωrω+1+jωτcτc)]mω2ϕr(L)ωr2−ω2+j2ζrωrω$
(65)

Note that one could as well define the relative motion transmissibility FRF for any other point (say, for point $x1$) throughout the beam (by simply setting $ϕr(x1)$ at the right hand side of Eq. 65). However, the motion at the tip of the beam is of particular interest, because it is the position of the maximum transverse displacement for the practical vibration modes. As a consequence, the vibratory motion at the tip of the beam plays an important role while deciding on the volume of the harvester.

Figure 9 shows the modulus of the relative tip motion FRF against excitation frequency. Considering the main graph, it is not easy to distinguish between the FRFs for different values of load resistance. However, as can be seen from the enlarged views in Fig. 9, there are considerable variations around the resonance frequencies. We observe the same short circuit and open circuit resonance frequency behaviors. Note that the uncoupled (purely mechanical) FRF is also provided, and expectedly, as $Rl→0$, the coupled FRF converges to the uncoupled FRF.

As the value of load resistance is increased from $Rl=102Ω$ to $Rl=105Ω$, the vibration amplitude at the short circuit resonance frequency $(47.8Hz)$ decreases considerably (by a factor of about 2.5). However, when $Rl$ is further increased to $106Ω$, the amplitude of vibrations at this frequency starts increasing. Therefore, it can be concluded that the vibration amplitude at a frequency does not necessarily show a monotonic behavior with increasing/decreasing load resistance, as in the case of the power output FRF. If we investigate the vibration amplitude at the open circuit resonance frequency $(48.8Hz)$ as $Rl$ is increased from $102Ω$ to $106Ω$, we see that the vibration amplitude first starts decreasing smoothly and then it starts increasing sharply.

Figure 1 allows investigating the variation of relative tip displacement amplitude with load resistance more clearly. As can be seen from Fig. 1, the relative tip vibration is insensitive to variations of load resistance in the low load resistance region and reasonably the vibration amplitude at the short circuit resonance frequency is higher in this region. As the load resistance is further increased, due to the electromechanical effects, the vibration amplitude at the short circuit resonance frequency is attenuated. One should be aware of the fact that this attenuation in the vibration amplitude at the short circuit resonance frequency is indeed more complicated than just damping. Considering the first enlarged view in Fig. 9, one can see that the peak moves from $47.8Hzto48.8Hz$. This is the reason of the attenuation of vibration amplitude at $47.8Hz$. At this point, one should expect some increase in the vibration amplitude at $48.8Hz$ and this is what we observe both in Figs.  910. As the load resistance increases, the peak, which used to be at the short circuit resonance frequency, moves toward the open circuit resonance frequency, causing not only an attenuation at the former frequency but also an increase at the latter frequency.

It is also worthwhile to investigate the power versus load resistance (Fig. 8) and the relative tip displacement versus load resistance (Fig. 1) trends simultaneously. At the short circuit resonance frequency, as the load resistance is increased gradually, the power output increases until the load resistance takes its optimum value. At the same time, the vibration amplitude is attenuated considerably. Further increase in the load resistance reduces the power output, which is associated with a slight increase in the vibration amplitude. At the open circuit resonance frequency, an increase in the load resistance first reduces the vibration amplitude slightly and the power increases at the same time. Then, in the vicinity of the optimum load resistance for excitation at open circuit resonance frequency, the vibration amplitude starts increasing and it increases with increasing load resistance with a high rate. The important observation is the opposite trend in the relative tip displacement for excitations at the short circuit and the open circuit resonance frequencies around their respective optimum resistive loads.

It should be mentioned that the relative tip motion FRF exhibits antiresonance frequencies. In the frequency range $0–1000Hz$, an antiresonance frequency is captured at $536.7Hz$ for short circuit conditions, which moves to $538.6Hz$ for open circuit conditions. It can be seen from Fig. 9 that, at this antiresonance frequency, the relative displacement at the tip of the beam is less than 1% of the base displacement. Note that this is the formal definition of an antiresonance frequency of an FRF in the vibration literature and this frequency should not be confused with the open circuit frequency of a vibration mode, which duToit et al. (6) call the antiresonance frequency in their SDOF model. Although it is not possible to see from Fig. 9, the FRFs are not identical at the off-resonance frequencies and they are slightly different. The intersections between the curves observed for the power FRF are also observed here in the relative tip motion FRF.

## Electromechanical Coupling and the Effect of Strain Nodes

Consider the general form of the circuit equation given by Eq. 34. The forcing term at the right hand side of this equation determines the amplitude of the electrical output. Other than the modal mechanical response (in the form of modal velocity response), the forcing term depends on the modal coupling term$φr$. As far as the mechanical equation of motion is concerned, the coupling information is included in the term $χr$, as can be seen from Eq. 27. Although $φr$ and $χr$ are defined separately in the electrical and the mechanical equations for convenience, they are not too different terms. Both $φr$ and $χr$ depend on the geometric parameters of the beam, Young’s modulus, and the electromechanical parameters of the PZT, and more importantly, the bending slope eigenfunction evaluated at the boundaries of the electrodes. If the type of the PZT and the geometry of the composite harvester beam are already decided, the last term related to the bending slope plays a very important role in shaping the coupled dynamics of the electromechanical system.

In the main derivation, since we assumed the electrodes to be covering the entire length of the beam, the parameter related to the bending slope eigenfunction was reduced to the slope at the free end of the beam (because the slope at the root of the beam is already zero). If the electrodes cover only the arbitrary region $x1⩽x⩽x2$, the term related to the bending slope eigenfunction appearing in the coupling term expressions, Eqs. 28,35, must be modified as follows:Display Formula

$∣dϕr(x)dx∣x=L→∣dϕr(x)dx∣x=x1x=x2$
(66)
which describes the dependence of electromechanical coupling on the locations of the electrodes for harvesting energy from the $rth$ vibration mode. If the continuous electrodes are located such that the slopes at their boundaries are identical for a vibration mode, one should expect zero or very low electromechanical coupling for vibrations at that mode. It should be kept in mind that good electromechanical coupling yields good electrical outputs in energy harvesting. Therefore, for vibrations at a particular mode shape, the aim must be to keep the bending slope difference at the boundaries of the electrodes maximum so that the maximum electrical output is extracted from the system.

The physical reason is the existence of the strain nodes where the distribution of the bending strain throughout the length of the beam changes sign, i.e., its phase changes by $180deg$(12). Strain nodes exist in all vibration modes other than the fundamental vibration mode of a cantilevered beam. Only in the first vibration mode, the strain distribution along the length of the beam is in phase. In all higher vibration modes, the curvature eigenfunction $(d2ϕr(x)∕dx2)$ of the beam (which is a measure of bending strain) changes sign as we move from the clamped end to the free end of the beam. We observe from Eq. 35 that the electromechanical coupling is proportional to the integral of the curvature over the length of the electrodes. This is why we end up with the bending slope eigenfunction $(dϕr(x)∕dx)$ evaluated at the boundaries of the electrodes. Due to the sign change in the integrand, cancellation occurs and the electromechanical coupling and consequently the electrical outputs are reduced drastically. In order to avoid this cancellation, segmented electrode pairs must be used to collect the electric charge developed at the opposite sides of strain nodes. The resulting electrical outputs of these electrode pairs will be out of phase by $180deg.$ and their leads must be combined accordingly for the maximum electrical output in the circuit (15).

We provide two simple demonstrations from our parametric case study given in the previous section for a specific resistive load $(Rl=105Ω)$. The first vibration mode has no strain nodes, whereas the second vibration mode has one strain node at $x=0.216L$, and the third vibration mode has two strain nodes at $x=0.132L$ and $x=0.497L$ (remember that $x$ is measured from the clamped end of the beam and $L$ is the length of the beam). Furthermore, for the second mode shape, it can be obtained that the slope at $x=0.471L$ is zero. For the third mode shape, one of the two positions where the slope is zero is $x=0.291L$(12). Note that all these numbers are for an Euler–Bernoulli beam without a tip mass, and existence of a tip mass alters these numerical values (15).

Figure 1 shows three voltage FRFs for different electrode configurations and the attention is given to the second vibration mode. In the first configuration, the continuous electrode pair covers the entire surface of the beam $(0, which was already the case in our main discussion. We know that it gives the best electrical output for the first vibration mode since the strain distribution along the length of the beam is in phase for this mode. However, we also know that the second vibration mode has a strain node at $x=0.216L$ and covering this point by continuous electrodes should cause some cancellation. In order to see this, we introduce the second configuration where the electrodes cover the region between the strain node of the second mode and the free end of the beam (that is, $0.216L). As can be seen from Fig. 1, this configuration improves the voltage output from the second vibration mode (by a factor of more than 1.4). If desired, another electrode pair can be used to cover the region $0 and these two outputs (which are $180deg$ out of phase) can be combined for the best output from the second vibration mode. The third configuration covers the region $0. The major strain distributions between $0 and $0.216L cancel each other, and as can be seen from Fig. 1, the resulting voltage output from the second vibration mode is reduced drastically for this configuration (by a factor of more than 70). The slight contribution comes from residues of the neighboring modes, particularly from the first mode. It should also be mentioned that, among these configurations, the first and the third vibration modes give the best voltage output for the first configuration where the entire surface of the beam is covered with continuous electrodes.

Next, consider Fig. 1 where the attention is given to the third vibration mode. Again, we consider three electrode configurations. The continuous electrode pair covers the entire length of the beam in the first configuration $(0. Although this configuration gives the highest voltage output for the first vibration mode, this is not the case for the third mode, which has two strain nodes at $x=0.132L$ and $x=0.497L$. In the second configuration, the electrodes cover the region between these strain nodes, which is $0.132L. As can be seen from Fig. 1, this configuration improves the voltage output from the third mode (by a factor of more than 1.2). One might as well obtain the voltage outputs from the remaining regions $0 and $0.497L, and combine all three outputs to obtain the maximum voltage output. Note that the voltage output from the electrode pair covering $0.132L will be $180deg$ out of phase and the outputs of the other two regions will be in phase. In the third configuration, the electrodes cover the region $0 for demonstrating the cancellation in the third mode clearly. Since the strain distributions in $0 and $0.132L cancel each other, the voltage output is attenuated by a factor of more than 13. The slight electrical output is due to the contributions from the neighboring modes, especially from the second vibration mode. One should also notice the antiresonance frequencies that show up in the FRF for the second and the third electrode configurations.

## Conclusions

In this paper, a distributed parameter electromechanical model for cantilevered piezoelectric harvesters is derived. The analytical formulation of the coupled system is based on Euler–Bernoulli beam assumptions. The harvester beam is assumed to be excited due to the translational motion of its base in the transverse direction with superimposed small rotation. In mechanical modeling, the internal strain rate damping (i.e., Kelvin–Voigt damping) and the external air damping are treated more accurately by defining separate damping coefficients.

The electromechanical equations are first derived for general transient base motions. In other words, the base motion is not restricted to be harmonic in time so that general coupled expressions for the mechanical response and voltage output are obtained. Then, the electromechanically coupled equations are reduced for the case of harmonic base translation with small rotation, and closed-form expressions are presented for the voltage, current, and power outputs as well as the coupled mechanical response of the harvester. The resulting equations are further reduced for the case of excitation around a natural frequency.

The analytically obtained expressions are then used in a parametric case study. In order to observe the frequency response behavior of the electrical outputs and the relative tip motion of the harvester, the FRFs which relate the voltage, current, power, and relative tip motion to the base motion are identified. These FRFs are plotted against frequency for a wide range of load resistance. Short circuit and open circuit conditions of the system are discussed. For a better understanding of short circuit and open circuit conditions, the electromechanical outputs are also plotted against load resistance for these two extreme cases of the resistive load. The mathematical modeling is based on the assumption of proportional damping (strain rate damping is assumed to be stiffness proportional, whereas air damping is assumed to be mass proportional). However, after describing how to identify the strain rate and air damping coefficients from experimental measurements, a simple relaxation is described for handling experimentally obtained nonproportional damping in the modal expansion.

Finally, electromechanical coupling and its relevance to the locations of the electrodes and the strain nodes are discussed. Once the geometry and the material of the bender are decided, the locations of the electrodes become important in determining the magnitude of the modal electromechanical coupling, which determines the magnitude of the electrical outputs. The issue of strain nodes of vibration mode shapes and cancellation of the electrical outputs due to covering the strain nodes of higher vibration modes with continuous electrodes are described with examples. Suggestions are made for increasing the electromechanical coupling and therefore the electrical outputs of the harvester.

## Acknowledgements

The authors gratefully acknowledge the support of the Air Force Office of Scientific Research MURI under Grant No. F9550-06-1-0326 “Energy Harvesting and Storage Systems for Future Air Force Vehicles” monitored by Dr. B. L. Lee.

## Appendix

Consider Fig. 1, which displays the cross section of the composite unimorph beam of Fig. 1. The width of the beam is denoted by $b$, the thickness of the PZT layer is $hp$, and the thickness of the substructure layer is $hs$. The procedure of finding the position of the neutral axis of a composite cross section is described in elementary strength of material texts (e.g., Timoshenko and Young (18)) and it requires transforming the cross section to a homogeneous cross section of single Young’s modulus (see Fig. 1). We take the PZT as the material of the transformed cross section and define the ratio of Young’s moduli as $n=Ys∕Yp$. In the transformed cross section, the width of the substructure layer is increased if $Ys>Yp$ or it is reduced if $Yp>Ys$. For demonstration, the typical case $Ys>Yp$ is assumed in Fig. 1 (which is also the case in our parametric case study) so that widening occurs in the substructure layer. Then, it is a simple practice to obtain the parameters defined in Fig. 1 (and therefore the position of the neutral axis) in terms of the parameters of Fig. 1 and Young’s moduli ratio $n$ as follows:

$hpa=hp2+2nhphs+nhs22(hp+nhs),hsa=hp2+2hphs+nhs22(hp+nhs),$
$hpc=nhs(hp+hs)2(hp+nhs)$
where $hpa$ is the distance from the top of the PZT layer to the neutral axis, $hsa$ is the distance from the bottom of the substructure layer to the neutral axis, and $hpc$ is the distance from the center of the PZT layer to the neutral axis. Note that the geometric parameters used in the main formulation ($ha$, $hb$, and $hc$) describe positions from the neutral axis rather than distances. Therefore, they can be expressed asDisplay Formula
$ha=−hsa,hb=hpa−hp,hc=hpa$
(A2)
It can be remembered from the main text that $ha$ is the position of the bottom of the substructure layer from the neutral axis (always negative), $hb$ is the position of the bottom of the PZT layer from the neutral axis (positive or negative), and $hc$ is the position of the top of the PZT layer from the neutral axis (always positive).

## References

View article in PDF format.

## Figures

Figure 3

Voltage FRF for five different values of load resistance (with the enlarged view of Mode 1 resonance showing the short circuit and open circuit behaviors)

Figure 4

Variation of voltage output with load resistance for base excitations at the short circuit and open circuit resonance frequencies of the first vibration mode

Figure 5

Current FRF for five different values of load resistance (with the enlarged view of Mode 1 resonance showing the short circuit and open circuit behaviors)

Figure 6

Variation of current output with load resistance for base excitations at the short circuit and open circuit resonance frequencies of the first vibration mode

Figure 7

Power output FRF for five different values of load resistance (with the enlarged views of Mode 1 and Mode 2 resonances showing the short circuit and open circuit behaviors)

Figure 8

Variation of power output with load resistance for base excitations at the short circuit and open circuit resonance frequencies of the first vibration mode

Figure 1

Unimorph piezoelectric energy harvester under translational and small rotational base motions

Figure 2

(a) The unimorph piezoelectric energy harvester used for the parametric case study and (b) a detail from its cross section (dimensions are in millimeters)

Figure 9

Relative tip motion FRF for the uncoupled system and for the coupled system with five different values of load resistance (with the enlarged views of Mode 1 and Mode 2 resonances showing the short circuit and open circuit behaviors)

Figure 10

Variation of relative tip displacement to base displacement ratio with load resistance for base excitations at the short circuit and open circuit resonance frequencies of the first vibration mode

Figure 11

Effect of the location of continuous electrode pair on the voltage FRF (focusing on the vibrations around the second mode)

Figure 12

Effect of the location of continuous electrode pair on the voltage FRF (focusing on the vibrations around the third mode)

Figure 13

(a) Cross section of a unimorph harvester and (b) the transformed cross section

## Tables

Table 1
Geometric, material, and electromechanical parameters of the sample harvester
Table 2
Short circuit and open circuit resonance frequencies of the harvester for the first three vibration modes

## Discussions

Some tools below are only available to our subscribers or users with an online account.

### Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections