Article

Journal of Intelligent Material Systemsand Structures23(2) 183–199� The Author(s) 2011Reprints and permissions:sagepub.co.uk/journalsPermissions.navDOI: 10.1177/1045389X11432656jim.sagepub.com

Nonlinear nonconservative behaviorand modeling of piezoelectric energyharvesters including proof mass effects

Samuel C Stanton1, Alper Erturk3, Brian P Mann2, Earl H Dowell2 andDaniel J Inman4

AbstractNonlinear piezoelectric effects in flexural energy harvesters have recently been demonstrated for drive amplitudes wellwithin the scope of anticipated vibration environments for power generation. In addition to strong softening effects,steady-state oscillations are highly damped as well. Nonlinear fluid damping was previously employed to successfullymodel drive dependent decreases in frequency response due to the high-velocity oscillations, but this article instead har-monizes with a body of literature concerning weakly excited piezoelectric actuators by modeling nonlinear damping withnonconservative piezoelectric constitutive relations. Thus, material damping is presumed dominant over losses due tofluid-structure interactions. Cantilevers consisted of lead zirconate titanate (PZT)-5A and PZT-5H are studied, and theaddition of successively larger proof masses is shown to precipitate nonlinear resonances at much lower base excitationthresholds while increasing the influence of higher-order nonlinearities. Parameter identification results using a multiplescales perturbation solution suggest that empirical trends are primarily due to higher-order elastic and dissipation nonli-nearities and that modeling linear electromechanical coupling is sufficient. This article concludes with the guidelines forwhich utilization of a nonlinear model is preferred.

KeywordsEnergy harvesting, piezoelectric, structural health monitoring, sensor, ferroelectric, actuator

Introduction

Piezoelectric energy harvesting has matured as a topicof consistent theoretical and experimental investigation.In particular, the bimorph harvester has been prolifi-cally studied. This type of energy harvester is typicallya cantilevered beam consisted of an electrically neutralsubstrate material with symmetric electroelastic lami-nates. A proof mass may be attached to tune the reso-nant frequency of the device upon coupling to a powerconditioning circuitry to extract the optimal amount ofpower from an ambient source of mechanical motion.

It is widely accepted that the linear piezoelectric con-stitutive relations are valid for weakly excited harvesters.However, several recent studies have shown otherwise(Stanton et al., 2010a, 2010b; Stanton and Mann). Thisarticle provides a comprehensive modeling approachand experiment to investigate the validity of modelinglinear behavior within piezoceramics. To first set thecontext, however, we review the relevant history ofpiezoelectric harvester modeling.

Sodano et al. (2004) and DuToit et al. (2005) wereamong the first to publish theoretical models for the

linear dynamics of piezoelectric energy harvesters withexperimental validation. Sodano et al. (2004) general-ized several approaches for piezoelectric sensors andactuators (Crawley and Anderson, 1990; Hagood et al.,1990; Smits et al., 1991) to an energy-harvesting schemeand showed that damping effects of power harvestingfollow that of a resistive shunt. DuToit et al. (2005)closely followed this work with a lengthy and detailedexposition encompassing not only a variationalapproach to modeling but also various microelectrome-chanical system (MEMS) fabrication and design issues,

1US Army Research Office, Durham NC, USA2Department of Mechanical Engineering, Duke University, Durham, NC,

USA3G. W. Woodruff School of Mechanical Engineering, Georgia Institute of

Technology, Atlanta, GA 30332-04054Department of Aerospace Engineering, University of Michigan, Ann

Arbor, MI

Corresponding author:

Samuel C Stanton, US Army Research Office, Durham NC USA

Email: samuel.c.stanton2.civ@mail.mil

resonance versus antiresonance, and nontrivial proofmass kinematics. Later, instead of the Rayleigh–Ritzapproach in the study by Sodano et al. (2004) andDuToit et al. (2005), Jiang et al. (2005) solved theeigenvalue problem for a flexural energy harvester witha point mass approximation from the governing partialdifferential equations where the electromechanical cou-pling appears in the natural boundary conditions. Theexact electromechanical resonance point is derivedfrom a determinant equation, and the ensuing analysisfocused on optimal impedance loads, bimorph geome-try, and proof mass size for increased power densities.Erturk and Inman (2008a, 2009a) later published linearmodels for both unimorph and bimorph harvesterstransmitting a proof mass and both a translational anda small rotary base excitation. The same authorsaddressed errors stemming from single-degree-of-free-dom (SDOF) models (Erturk and Inman, 2009b) andother questionable results (Erturk and Inman, 2008b,2008c).

While the preceding articles have formed a basis foraccurate energy harvester models, several recent articleshave explicitly focused on proof mass effects. Kim et al.(2010) published experimental results confirming thevalidity of the theory presented in the study by DuToitet al. (2005) with a macroscopic device and two differentnonsymmetric proof masses. Yu et al. (2010) also stud-ied nonsymmetric proof masses but added nonvanishinglongitudinal displacement and rotary inertia within thecantilever itself in the theoretical model. Finite elementmethods were used to solve the equations of motion,and a point mass assumption, as expected, inaccuratelycaptured the resonance. Interestingly, the experimentaldata from both articles show poor correlation in someinstances with the linear modeling. In fact, Kim et al.(2010) suggest that nonlinear piezoelectric effects werethe culprit, while Yu et al. (2010) do not hypothesizewhy their data exhibit a clear softening and reducedamplitude effects reminiscent of nonlinear damping (seeFigure 7 in the study by Yu et al. (2010)).

Piezoelectric nonlinearity in a harvester was first stud-ied by Hu et al. (2006). A numerical analysis of a piezo-electric plate showed a hardening response due tothickness/shear mode vibration. The physical basis fortheir model was the cubic theory of nonlinear electroelas-ticity (Maugin, 1985; Tiersten, 1993; Yang, 2005) butwith linear dissipation and linear coupling. Nonlinearcoupling in the context of piezoelectric beam actuatorswas modeled by Von Wagner and Hagendorn (2002),and the nonlinear response to weak electric fields wasexplained by a combination of nonlinear elastic and cou-pling phenomena. Triplett and Quinn (2009) theoreticallystudied the influence of nonlinear strain-dependent cou-pling in harvesters. While the investigation utilized para-meters that yielded hardening responses, harvesterpiezoceramics typically exhibit softening curves (Priya etal., 2001; Stanton et al., 2010a, 2010b). Nevertheless, the

nonlinear coupling could indeed influence already non-linear resonance curves due to elastic effects. In all afore-mentioned analyses, however, linear mechanicaldissipation is presumed. Stanton et al. (2010a, 2010b)studied the nonlinear response of bimorph harvestersboth theoretically and experimentally and modeled non-linear damping through quadratic air drag. However,without a detailed analysis of the fluid loading in theexperiments, it is difficult to maintain confidence thatthis is the primary source of dissipation.

Motivated by empirical trends and the existing liter-ature concerning weakly excited piezoelectric materials,this article seeks to comprehensively model both non-linear elasticity and dissipation in piezoelectric energyharvesters. The section titled ‘‘Nonlinear and noncon-servative modeling’’ provides a derivation of a SDOFmodel from Rayleigh–Ritz methods for the piezoelec-tric harvester to include higher-order piezoelectricrelations and nonconservative work expressions.Experimental validation is presented next and followedby metrics for which utilization of a nonlinear model ispreferred in the section titled ‘‘Linear versus nonlinearmodeling.’’ The section titled ‘‘Summary and conclu-sions’’ provides a summary and a list of our main con-clusions along with the directions for future work.

Nonlinear and nonconservative modeling

This section employs the Rayleigh–Ritz procedure toderive a SDOF model for the harvester. In the Rayleigh–Ritz procedure, a modal expansion of the spatiotemporaldeflection of a cantilever beam is substituted into theappropriate conservative and nonconservative energyexpressions, and subsequent application of Hamilton’sextended principle yields governing ordinary differentialequations. Provided that our experimental devices areexcited within close proximity of the first resonance ofthe coupled electromechanical system, a single-modeapproximation is presumed to be sufficiently accurate formodeling and identification purposes.

Figure 1 illustrates a prototypical harvester con-sisted of a bimorph cantilever beam with an attachedmass. The harvester is composed of a brass shim ofthickness hp with symmetric piezoelectric laminatesof thickness hp each. The transverse deflection isdescribed by the coordinate w x, tð Þ and is subject to abase displacement z tð Þ. A proof mass is affixed suchthat the center of mass is at the tip of the beam. Thekinetic energy is distributed through the volume ofthe cantilever as

T =1

2rb

ðVb

_w + _zð Þ2dVb +1

2rp

ðVp

_w + _zð Þ2dVp ð1Þ

where w x, tð Þ is the transverse deflection, rb, rp, Vb, andVp are the respective material densities and total volume

184 Journal of Intelligent Material Systems and Structures 23(2)

for the brass substrate and piezoelectric laminates. Thecross section of the cantilever is rectangular with widthb and total thickness of hs + 2hp, where hs is the thick-ness of the substrate, and hp is the thickness of the indi-vidual piezoelectric laminates. To formulate the kineticenergy of the rigid body at the boundary, consider aposition vector from the base of the cantilever to centerof mass to the proof mass. Considering Figure 2, thisvector may be written as

rMt= L + oxw9 L, tð Þ� �

e1 + z + w L, tð Þ � ox + hp + 12hs

� �� �e3

ð2Þ

where the kinematic assumptions reflect negligible long-itudinal motion and first-order rotations (i.e. cos u’1

and sin u’u). The proof mass kinetic energy follows as

T Mt=

1

2_rTMt

Mt _rMt+

1

2Jyy _w92

� �����s = L

=1

2Mt _w + _zð Þ2 +

1

2Jo _w92

� �����s = L

ð3Þ

where Jyy is the inertia of the proof mass about its cen-ter of gravity and Jo = Jyy + 2Mto

2x is the rotary inertia

of the proof mass that incorporates parallel axis theo-rem as a result of expanding _rT

MtMt _rMt

. The total kineticenergy of the harvester-mass system is found by sum-ming equations (1) and (3).

The stored energy within the cantilever is composedof the bending potential of the brass shim and the elec-troelastic enthalpy of the piezoelectric laminates. Storedenergy within the brass shim is presumed linear and

varies in proportion to the square of the longitudinalstrain

U =1

2cxx, b

ðVb

S2x dVb ð4Þ

where cxx, b is the brass elastic modulus and Sx denotesthe strain along the x-axis.

The conservative potential of the piezoceramics isexpressed in terms of an electric enthalpy function. Forlinear piezoelectric behavior, this function consists ofpurely elastic and dielectric terms to describe themechanical and electrical behaviors along with a cou-pling parameter that yields electromechanical effects.Mathematically, this is written in a reduced form for abeam-like device as

H=1

2

ðVp

cExxS2

x � 2ezxSxEz � eSzzE

2z

� �dVp ð5Þ

where cExx is the elastic modulus, ezx is the piezoelectric

coupling coefficient, eSzz is the dielectric constant, and

Ez is the generated electric field perpendicular to thestrain. However, this form is inadequate to explainempirical trends. Experimental results within the litera-ture relating to piezoelectric plates, rods, and beamsdriven by weak electric fields demonstrate nonlineareffects typically characterized by frequency-dependentamplitude responses, nonlinear damping, and evenbifurcation phenomena. Figure 3 shows the generaltrend in the voltage response normalized by the driving

R

e1

z(t)

Piezoceramic

Brass

Proof Mass

e3

w(L,t)

Figure 1. Illustration of the energy harvester studied. The cantilever is excited by a base motion z tð Þ and transversely deflects a dis-tance of w x, tð Þ.

Stanton et al. 185

amplitude for one of the experimentally identified can-tilevers. All experimental tests resulted in softening fre-quency response curves with complicated dampingcharacteristics. To ascertain a first principles–basedmodel for the energy harvester, a fully nonlinear exten-sion of equation (5) following Samal et al. (2006a) wasapplied that included third- and fourth-order nonlinearcoupling and dielectric effects in addition to higher-order structural nonlinearities. In a similar fashion,Stanton et al. (2010a, 2010b) considered nonlinear con-stitutive relations for modeling piezoelectric harvest-ers with dielectric nonlinearities disregarded butcoupling nonlinearities retained. In these studies,experimental identification using a nonlinear leastsquares algorithm converged upon several localminima indicating nonvanishing nonlinear coupling,although a global minimum was unattainable. In thepresent analysis, exhaustive simulation studies includ-ing proof mass effects revealed that numerical trendsresulting from conservative nonlinear coupling andnonlinear dielectric effects were counter to the experi-mental measurements when a cubic structural nonli-nearity was kept to explain the softening frequencyresponse. This result is in agreement with the cubictheory of electroelasticity and reduces model com-plexity to that of an electromechanical Duffing oscil-lator with linear coupling. For high accelerationvalues, the inclusion of a proof mass suggested atrend indicative of higher-order structural nonlineari-ties (fifth- and seventh-order harmonics) along withweakened nonlinear damping effects attributable tothe inertial influence of the proof mass. To model theempirical trends in this study, we employ a nonlinearelectroelastic enthalpy expression that includeshigher-order elastic nonlinearities up to fifth orderbut disregards higher-order coupling and electricaleffects. Accordingly, we now have

H=1

2

ðVp

cExxSx � 2ezxSxEz � eS

zzE2z +

2

3c3S3

x +1

2c4S4

x

+2

5c5S5

x +1

3c6S6

x )dVp

ð6Þ

where c3, c4, c5 and c6 are nonlinear stiffness coeffi-cients. From equation (6), the nonlinear constitutiverelations are derived from the compatibility equations

Tx =∂H∂Sx

= cExxSx � ezxEz + c3S2

x + c4S3x + c5S4

x + c6S6x ð7Þ

and

Dz = � ∂H∂Ez

= ezxSx + ezzEZ ð8Þ

where Tx is the stress and Dz is the electric displacement.With this formulation the following relations hold

∂2H∂Sx∂Ez

=∂Tx

∂Ez

= � ∂Dz

∂Sx

=∂2H

∂Ez∂Sx

ð9Þ

and the necessary and sufficient thermodynamic condi-tions for the existence of an enthalpy function are satis-fied. Due to the symmetry of the bimorph, oddnonlinearities in the enthalpy function will cancel.Thus, presuming perfect symmetry in the bimorphrestricts nonlinear effects to cubic and fifth-order struc-tural effects in the equation of motion.

The gravitational potential of the beam is also mod-eled for completeness and to account for induced asym-metry, although this was found to be negligible. Futureinvestigations with heavier proof masses may inducestatic deformation of the beam, and the potential fieldfor modeling is given by

Ug = g rb

ðVb

w + z + oxð Þdx + rp

ðVp

w + z + oxð Þdx

264

375

+ Mtg w L, tð Þ + z½ � ð10Þ

where the datum is taken as the center of gravity of theproof mass, and we have again presumed small-anglerotations in determining the potential of the proofmass.

Consistent experimental observation of decreasedfrequency response amplitudes to higher driving ampli-tudes in a variety of piezoelectric cantilevers by theauthors necessitates modeling nonlinear damping. Inaddition to conservative energy considerations, thisresult suggests nonvanishing virtual work contributionsto the oscillatory dynamics of the beam as well. In thepast, quadratic fluid damping successfully modeled thisphenomena (Stanton et al., 2010a, 2010b), but theunderlying physical mechanism of this influence at thescale is a relatively open question. Although air damp-ing is a dominant source of dissipation in

ox

L

z(t)+w(L,t)

x

θ

(ox+hs+hp/2)cos( )θ

ox

oxsin( )θ

Figure 2. Closer view of the proof mass kinematics for deter-mining the position vector from the origin to the center of mass.

186 Journal of Intelligent Material Systems and Structures 23(2)

microcantilever systems (Zhang and Meng, 2005), thelarger amplitudes in our macroscopic systems may ren-der air damping less of a primary contributor to the dis-sipative influences despite the high-response velocityoscillations. In this study, we model damping thatincreases with the response amplitude in an alternativemanner by presuming material dissipation dominatesair damping. Note, however, that quadratic air dampingis critical for limiting the growth of oscillations in canti-lever harvesters engaging geometric and inertial nonli-nearity (Daqaq et al., 2009). Considering theexperimental observations within the study by Samal etal. (2006b) of nonlaminate piezoceramics, the source ofnonlinear damping in this article is presumed to be dueto viscoelastic effects within the piezoceramics. Thebrass substrate, however, experiences a virtual work inproportion to the time rate of change in strain of theform

dWb =

ðVp

hb_SxdSxdVp ð11Þ

and is attributable to structural damping, losses at theclamp location, and thermoelastic dissipation, amongother sources (Kamel et al., 2010).

Ikeda (1996) introduced dissipation within electroe-lastic media by extending the constitutive equationswith time derivatives of the conservative basis variables.Later, Von Wagner (2003) employed this method to anonlinear set of constitutive relations, and his techniqueis applied here. Accordingly, we extend the constitutiverelations as

Tx =∂H∂Sx

+ Tnc ð12aÞ

and

Dz = � ∂H∂Ez

+ Dnc ð12bÞ

where Tnc and Dnc denote nonconservative stress andelectric displacement that contain both linear and non-linear dissipative terms defined by

Tnc = hp_Sx � g _Ez + za

_S3

x

� �+ zb

_S5

x

� �ð13aÞ

and

Dnc = g _Sx + v _Ez ð13bÞ

Thus, nonvanishing virtual work due to piezoelectriclosses can be incorporated into Hamilton’s extendedprinciple through

dWp =

ðVp

TncdSx + DncdEzð ÞdVp ð14Þ

With all conservative and nonconservative energyformulations defined, we are in a position to apply var-iational calculus (Crandall et al., 1968) to derive theequations of motion. First, the kinetic energy (equa-tions (1) and (3)), potential energy (equations (4) to (6)and (10)), and the virtual work (equations (11) and(14)) can be further reduced by considering a single-mode approximation of the Euler–Bernoulli strain

Sx = � zr tð Þf00 xð Þ ð15Þ

where ()9 is shorthand for d=dx, r tð Þ is the displace-ment, and f xð Þ is the first spatial mode of a fixed-freecantilever with proof mass boundary conditions (seeequations (12) and (16) within Erturk and Inman(2009a) where Jo = It). The eigenvalues and modes aswithin the study by DuToit et al. (2005) were also cal-culated with a proof-mass/cantilever overlap of ox butwere found to offer no improvement to the center ofmass approximation.

The electric field can be expressed in terms of an elec-tric potential function u x, tð Þ for the upper and lowerlaminates as

Euz = � u9 x, tð Þ

hp

andElz =

u9 x, tð Þhp

ð16Þ

where the change in sign is due to the opposite polingdirections as a consequence of a series connectionbetween bimorph layers. Next, by defining theLagrangian as

L = T + T Mt�H� Ug ð17Þ

Frequency (Hz)

Nor

mal

ized

Vol

tage

Increasing Base Excitation

Figure 3. Illustration of the empirical trend in nonlinear reso-nance curves for the piezoelectric harvester. The output voltageis normalized by the drive amplitude to better portray the non-linear response.

Stanton et al. 187

incorporating equations (15) and (16), and integratingover the volume of the beam, we obtain a single-modeapproximation for the harvester dynamics. The varia-tion of the Lagrangian becomes

dL’ M _r + Fz _zð Þd_r + �Kr � Gr3 � Hr6 +Yu9� Fg

� �dr

+ � ~Yr � ~Cu9� �

du9 ð18Þ

where the constant coefficients are defined in Appendix1. Similarly, the virtual work is

dW’ Da + Dbr2 + Dcr4� �

_r � ~Lu9h i

dr + ~L_r + xu9h i

du9

ð19Þ

where linear and nonlinear structural dissipation, con-version losses, and dielectric losses are present. Theelectrical network introduces electronic losses across aresistive load to eliminate complications arising fromvoltage conditioning for a storage circuit. Thisimproves identification of nonlinear beam parametersand serves to represent the equivalent impedance of amore complex storage network. Since the voltage is dis-tributed across two laminate layers, we define the elec-tric field in terms of the circuit flux linkage u = 1

2x _l tð Þ

to facilitate unified variational modeling. Hamilton’sextended principle for electromechanical systems istherefore

ðt1t0

dL+ dW + Idlð Þdt = 0 ð20Þ

where dL and dW are as before and Idl is a generalizedcurrent. For a linear resistor, this is given by

Idl =1

R

_ldl ð21Þ

where R is the load resistance. Applying the calculus ofvariations and collecting terms common in dr and dl

yield a nonlinear equation of motion for the displace-ment of the cantilever

M€r + Da + Dbr2 + Dcr4� �

_r + Kr + Gr3 + Hr5

�YV � Dd_V = � Fg + Fz€z tð Þ ð22Þ

and a linear equation for the harvesting circuit

De€V + C _V +

1

RV +Y_r + Dd€r = 0 ð23Þ

where we have also replaced the time derivative of theflux-linkage coordinate _l with voltage V for conveni-ence. The form for the base acceleration €z tð Þ is general,but for this particular study, it is harmonic with an exci-tation frequency near fundamental resonance such that

€z tð Þ= Z cos Oet, where Z is the amplitude and Oe is thefrequency of excitation.

A dimensionless form for equations (22) and (23)may be derived by substituting a characteristic time,length, and voltage

Tc =

ffiffiffiffiffiM

K

r, Lc =

ffiffiffiffiffiffiffiK

jGj

s, Vc =

Lc

Tc

ffiffiffiffiffiM

C

rð24Þ

into the equations of motion. Accordingly, the remain-der of this article examines the following dimensionlessequations for the mechanical domain

€x + ma + mbx2 + mcx4� �

_x + x� x3 + bx5

�uv� k _v = � fg + fz cos Ot ð25aÞ

and electrical domain

me€v + _v + mdv + u _x + k€x = 0 ð25bÞ

where x tð Þ= r tð Þ=Lc, v tð Þ= V tð Þ=Vc, and the sign pre-ceding x3 are negative for convenience due to the soften-ing frequency response of the beam. The dimensionlessparameters are defined as

ma =Da

M

ffiffiffiffiffiM

K

r, mb =

Db

jGj

ffiffiffiffiffiK

M

r, mc =

DcK

G2

ffiffiffiffiffiK

M

r,

md =1

RC

ffiffiffiffiffiM

K

r, me =

Dd

C

ffiffiffiffiffiK

M

r, b =

HK

G2

u =YffiffiffiffiffiffiffiCKp , k =

Dc

MC, fg =

Fg

KjGj

and

fz = F Gj jZ 1

K2

ffiffiffiffiffiffiffiK

jGj

s

while the dimensionless excitation is O =Oe

ffiffiffiffiffiffiffiffiffiffiffiM=K

p.

Multiple time-scaling perturbation solution

This section develops an analytic solution for harvesterdynamics based on perturbation methods to enable para-meter identification. In particular, the method of multiplescales is employed such that a weakly nonlinear behaviorcan be predicted by separating the dimensionless displa-cement and voltage into a third-order summation ofterms proportional to a small expansion parameter e andcorresponding separate time scales as

x tð Þ=X1

k = 0

ekxk T0, T1ð Þ and v tð Þ=X1

k = 0

ekvk T0, T1ð Þ

ð26Þ

and time derivative operators become

188 Journal of Intelligent Material Systems and Structures 23(2)

d

dt=

∂

∂T0

+ e∂

∂T1

+O e2� �

ð27Þ

and

d2

dt2=

∂2

∂T20

+ 2e∂2

∂T0∂T1

+O e2� �

ð28Þ

Following Masana and Daqaq (2011), we order thecoupling parameter u in the mechanical dynamics toappear at higher-order perturbation corrections butto an order of zero in the electrical network. Thisallows for the electrical circuit to be studied as a dri-ven linear system in terms of the zeroth-order pertur-bation correction for displacement. We also balanceall damping, nonlinear stiffness terms, and harmonicforcing to appear at the first-order correction. Theconstant force due to gravity is ordered to O 1ð Þ, sothat these effects are not lost in the perturbation solu-tion. Furthermore, extensive simulation trials indi-cated that for the voltages and strains induced in ourexperimental devices, the terms encompassing dissipa-tive coupling and dielectric losses were negligible.Hence, we choose to order these terms such that thederived solvability condition does not reflect theseinfluences. As a result, our focus is now on the fol-lowing system of equations:

€x + x + fg = � e ma + mbx2 + mcx4� �

_x�

�x3 + bx5 + uv� fb cos Ot�+ e2mc _v ð29aÞ

and

_v + mev = � u _x� e me€v + mc€xð Þ ð29bÞ

where e serves as a bookkeeping parameter that will beset to unity at the end of our analysis. Multiple time-scaling perturbation solution requires the excitation fre-quency O to instead be expressed in terms of the near-ness to the unit resonance by

O= 1 + es ð30Þ

where es is a small detuning away from resonance.Substituting equation (26) into equations (29a) to (29b)and collecting terms independent of e gives

D21x0 + x0 = � fg ð31aÞ

and

D1v0 + mdv0 = � uD0x0 ð31bÞ

where Dnk denotes the nth partial derivative with respect

to the kth time scale (i.e. Dnk = ∂Tn

k =∂Tnk ). Solving for x0

yields simple harmonic motion offset by fg for the dis-placement, while the voltage equation yields a forcedresponse due to x0 with an exponential decay. Hence,

in steady state, we have solutions for the lowest ordercontributions

x0 = A T1ð ÞejT0 � fg + c:c: ð32aÞ

and

v0 = � j

j + me

uA T1ð ÞejT0 + c:c: ð32bÞ

where j =ffiffiffiffiffiffiffi�1p

, A(T1) is a complex amplitude in theslower time scale and c.c. indicates a complex conju-gate. Substituting the order zero solution into the equa-tions at first order in e gives the first perturbationcorrection

D20x1 + x1 = � 2D0D1x0 + ma + mbx2

0 + mcx4� �

D0x0

��x3

0 + bx50 � uv0 + fz cos T0 + sT1ð Þ� ð33Þ

Upon consideration of the zeroth-order perturbationcorrections to the displacement and voltage (equations(32a) and (32b)) in equation (33), the solvability condi-tion is derived by disregarding higher harmonics andsetting secular terms to zero. Hence, we have

0 = Pr � jPið ÞA + Qr � jQið ÞA2�A

� Sr + jSið ÞA3�A2 � 2jA9 � 1

2fze

jsT1 ð34Þ

where the over bar denotes a complex conjugate, ()9now represents a derivative with respect to the slowertime scale T1, and the model parameters are collectedas

Pr = f 2g 3� 5bf 2

g

� �� u2

1 + m2d

ð35Þ

Pi = ma + mbf 2g + mcf 4

g +mdu2

1 + m2d

ð36Þ

Qr = 3� 30bf 2g ð37Þ

Qi = mb + 6mcf 2g ð38Þ

Sr = 10b ð39Þ

and

Si = 2mc ð40Þ

Substituting a polar expression for the complexamplitude

A =1

2aejc ð41Þ

where a and c indicate a real-valued amplitude andphase, respectively, and separating the right-hand side ofequation (34) into real and imaginary components gives

Stanton et al. 189

ac9 = � 1

2Pra�

1

8Qra

3 +1

32Sra

5 +1

2fz cos sT1 � cð Þ

ð42aÞ

and

a9 = � 1

2Pia�

1

8Qia

3 � 1

32Sia

5 � 1

2fz sin sT1 � cð Þ

ð42bÞ

Transformation to an autonomous system isachieved by letting f = sT1 � c, and the steady-statefrequency response curve is derived by setting all timederivatives in the system above to zero. Accordingly,steady-state solutions obey the algebraic system

fz cosf = 2s + Prð Þa +1

4Qra

3 � 1

16Sra

5 ð43aÞ

� fz sinf = Pia +1

4Qia

3 +1

16Sia

5 ð43bÞ

where the phase can be eliminated by squaring and add-ing both equations. The steady-state magnitude andphase are therefore

a2 =f 2z

2s + Pr + 14Qra2 � 1

8Sra4

� �2+ Pi + 1

4Qia2 + 1

8Sia4

� �2

ð44Þ

and

tan f =Pi + 1

4Qia

2 + 18Sia

4

2s + Pr + 14Qra2 � 1

8Sra4

ð45Þ

Alternatively, one can solve for the frequency detun-ing as

s = � 1

2Pr �

1

8Qra

2 +1

16Sra

4

6f 2z

a2� Pi +

1

4Qia

2 +1

8Sia

4

2" #1=2

ð46Þ

The harvester’s nonlinear resonance has severalinteresting features in comparison to a classic Duffingoscillator. First, as a coupled electromechanical system,the system resonant frequency is shifted a value Pr

away from the purely mechanical cantilever resonance.In view of equation (35), a linear harvester is shifted(u2=1 + m2

d) as the energy transfers through u betweenthe electrical and mechanical domains. By consideringa static gravitational load, it is also not surprising thatconsideration of this preload also shifts resonance awayfrom that predicted by Euler–Bernoulli beam theory.Second, the linear dissipation term Pi is augmented bythe electrical dissipation through the dimensionlessload md . Including fifth order, structural effects arecaptured through the term Sr, and third- and fifth-order dissipation appears within Qi and Si, respectively.

The manifestation of these results in the mathematics isin harmony with current experimental studies andthose within the literature (Stanton et al., 2010a,2010b).

Although equation (44) is tenth order in a, only fewsolutions will have no imaginary component and arehence physically meaningful. With only a cubic nonli-nearity, the existence of three solutions is a harbingerof a saddle node bifurcation and jump phenomena inthe harvester response. Because our experimental trialsdid not indicate jump phenomena would occur for therange of drive amplitudes studied, we complete our per-turbation analysis and forgo a discussion on linear sta-bility of the steady-state motion. Despite the fifth-ordernonlinearity, determining the stability of steady-statemotion may be accomplished by straightforward exten-sion of methods for typical Duffing oscillators. Whilethe study by Nayfeh and Mook (1979), among manyothers, details procedures for doing so, a good stabilityanalysis and discussion in the context of a nonlinearharvester can be found in a recent article by Daqaqet al. (2009).

Experimental investigation andidentification

This section describes the series of experiments thatwere conducted to assess the validity of the proposedmodel and to enable parameter identification. Twobimorph cantilevers, one symmetrically laminated withPZT-5H and another with PZT-5A (model numbersT226-H4-203X and T226-A4-203X, respectively, andmanufactured by Piezo Systems Inc.), were separatelyinvestigated. However, comparative analysis was facili-tated by affixing each bimorph to a clamp with similaroverhang lengths and proof masses and across the sameelectrical impedance load. A picture of a typical experi-ment is shown in Figure 4. Base motion was recordedby an accelerometer (PCB Piezotronics U352C67), andthe cantilever’s transverse tip velocity relative to thefixed reference frame is measured using a laser vibrom-eter (Polytec OFV353 laser head with OFV3001 vib-rometer). The cantilevers were never unclamped duringthe course of the experimental investigation andremained secured to the clamping apparatus as shownin Figure 4, while successive proof masses were affixed.Table 1 lists the plane-stress elastic, piezoelectric, anddielectric properties of each beam, as well as the proofmass properties. The dimensionless mass ratios

a =Mt

mLð47Þ

where m = b(rshs + 2rphp) is the mass per unit length ofthe cantilever. The values of m in the experiments were0, 0.298, and 0.976 for PZT-5H and 0, 0.293, and 0.972for PZT-5A. In all, 6 different cantilever configurations

190 Journal of Intelligent Material Systems and Structures 23(2)

were tested over 5 distinct frequencies and for 11 differ-ent base excitation values for PZT-5H and 9 base exci-tation values for PZT-5A.

Isolating the source of nonlinearity to piezoelectriceffects was paramount in designing the experiment.Hence, electrical impedance due to a resistive load in theelectronics was utilized so as not to invigorate higher-order harmonics stemming from nonlinear AC to DCelectronic conversion or switching techniques. In allexperiments, R = 100 Kohm. Also, the cantilever dimen-sions were selected such that ensuing small amplitudeoscillations despite large drive amplitudes ensures thevalidity of linear Euler–Bernoulli beam theory. In otherwords, no geometric nonlinearities such as nonlinearinertia or third-order structural nonlinearity (which ishardening in the first vibration mode of a cantilever)interfere with the already nonlinear cantilever response.

Linear domain frequency response functions werefirst measured for each case to identify resonance andlinear damping parameters. Linear damping is pre-sumed proportional to the cantilever resonance suchthat Da = 2zvn, where z is a dimensionless damping fac-tor and vn is the uncoupled fundamental frequency ofthe beam. Afterward, highly sampled (50 kHz) steady-state oscillations near resonance were recorded forincreasing levels of base excitation. The frequency con-tent of the data was analyzed via a power spectral den-sity (PSD) function (in MATLAB), and the magnitudeof the harmonics was found by Parseval’s theorem as

X = 2

ðvb

va

PX vð Þdv

24

35

1=2

ð48Þ

where X is the signal amplitude, PX is the PSD, and v

is the frequency. The magnitude of the acceleration, tipvelocity, and voltage signals at the first excitation fre-quency was observed to be of the orders of magnitudegreater than odd integer multiples of the same

frequency as expected, and reconstruction of the timeseries signal was coincident with the recorded timeseries for all driving amplitudes. Hence, a single-modeapproximation is indeed sufficient to model thedynamics of the harvester. Average base accelerationlevels for the PZT-5H cantilever (with and without aproof mass) were 60 mg, 145 mg, 230 mg, 310 mg,430 mg, 560 mg, 840 mg, 1.12g, 1.4g, 1.7g, and 2g,while the PZT-5A cantilever was excited at 50 mg,100 mg, 200 mg, 300 mg, 400 mg, 800 mg, 1.2g, 1.6g,and 2g.

The nonlinear resonance curve given by equation(46) is used in conjunction with experimental data to fitparameters in accordance with a nonlinear least squaresoptimization algorithm. This has been previously donein the study by Stanton et al. (2010a, 2010b), and weagain stress the importance of simultaneously studyingthe response in both the mechanical and electricaldomains when exploring the variety of local minimasuch an algorithm may converge upon. Good experi-mental data fitting in the voltage domain does notimply success in matching the displacement amplitude.For example, in the course of the present investigation,a fully nonlinear enthalpy function was studied with upto fourth-order electromechanical and dielectric effects.Initially, it was found that higher-order electromechani-cal coupling proportional to S2

x E2x generated excellent

agreement in the voltage domain only to be negated byextremely suppressed displacement curves. Hence, alloptimization routines were examined for agreement inboth voltage and mechanical response. After a widerange of analytical and numerical investigation span-ning all the possible nonlinearities, it was found thatwith the new nonlinear structural damping model,parameter identification with conservative

Table 1. Geometric and material properties of the bimorphcantilevers.

Parameter PZT-5Hbeam

PZT-5Abeam

Length (L) (mm) 24.06 23.82Width (b) (mm) 6.4 6.4Brass thickness (hb) (mm) 0.140 0.140Brass mass density (rb) (kg/m3) 9000 9000Brass elastic modulus (cxx, b) (GPa) 105 105PZT layer thickness (hp) (mm) 0.265 (each) 0.265 (each)PZT mass density (rp) (kg/m3) 7500 7750PZTelastic modulus (cp

xx) (GPa) 60.6 61Piezoelectric constant (ezx) (C/m2) 216.1 210.1Permittivity constant (es

zz) (nF/m) 25.55 13.27Mass 1 total mass (Mt) (kg) 0.24 0.24Mass 2 total mass (Mt) (kg) 0.787 0.795Mass 1 side length (2ox) (mm) 3.2 3.2Mass 2 side length (2ox) (mm) 4.7 4.7Mass 1 rotary inertia (Jo) (kg m2) 1:2903310�9 1:2903310�9

Mass 2 rotary inertia (Jo) (kg m2) 8:7386310�9 8:6549310�9

(1) Accelerometer

(2) Piezoelectric Beam

(3) Proof Mass

(4) Dissipative Loads

1

2

3

4

Figure 4. Experimental setup for all six test configurations.

Stanton et al. 191

nonlinearities due to elasticity alone was sufficient toachieve good experimental fits. This is plausible in viewof the weak electric field generation (at most 25–30 V/mm). Weak fields of similar magnitudes examined inseveral studies concerning piezoceramics for ultrasonicmotor applications (Parashar and Von Wagner, 2004;Samal et al., 2006b) have also not been strong enoughto evoke nonlinear coupling and dielectric effects butmaintain strong elastic and damping nonlinearity.

Figure 5 shows the strong model agreement for can-tilever beams with no attached proof mass. Properparameter identification is assured given the theoreticalagreement in both the voltage response and tip displace-ment. Each graph presents several nonlinear resonancecurves to demonstrate the softening and nonlineardamping trends. Figure 5(a) and (c) displays the resultsfor PZT-5H, and Figure 5(b) and (d) displays theresults for PZT-5A. To highlight the importance ofnonlinear modeling, the frequency response functionsfor the linear model are also shown. Both cantileversbegin to exhibit nonlinear resonance curves beyondZ’0:6g, a point beyond which the nonlinear dampingand elasticity conspire to shift and limit the frequencyresponse well away from that predicted by a linearmodel. With the nonconservative modeling, parameteridentification yields a larger fourth-order elasticity coef-ficient for both PZT-5H and PZT-5A than the valuesreported in prior studies employing a nonlinear air dragmodel (Stanton et al., 2010a, 2010b). However, in keep-ing with previous results, PZT-5A is verified to main-tain a larger softening effect in comparison to PZT-5H,as indicated by the larger value for c4 in Table 2. Forthe full range of experimental tests, modeling cubicnonlinearities due to x3 and mbx2 _x were found to be suf-ficient for both cantilevers with subtle deviation nearZ = 2g. Accordingly, with no proof mass, a nonlinearlydamped Duffing oscillator is accurate for a large rangeof base acceleration levels. The section titled ‘‘Linearversus nonlinear modeling’’ provides a quantitativeapproach for determining the accuracy of the nonlinearmodel.

With the addition of a proof mass, higher-ordernonlinearities begin to influence the harvester motion.In fact, in their analysis of a piezoelectric shell, Samalet al. (2006b) noted a likely need for damping termslarger than third order. While our derivation includedfifth-order damping and elasticity, it was determinedthat disregarding seventh-order nonlinearities began toreduce the model accuracy for the highest accelerationlevels and the largest proof mass. Considering realisticapplications are not likely to experience such largeacceleration drives, we instead focus on the results forwhich the fifth-order model retains validity. Figure 6gives the results for the PZT-5H beam for two proofmasses and the same five base acceleration levels. Theincreased inertia in the second proof mass case gener-ates a larger amplitude motion and hence higher

voltages in comparison to the first proof mass. Thenonlinear resonance curves coincident with the experi-mental data include terms identified proportional to b

and mc. Overlayed on the same plot are nonlinear reso-nance curves for the model with b = mc = 0. For Z .1g,the resonance curves begin to shift to the right and thecubic model breaks down. This is especially true whena is closer to unity. The increased inertia also assuagesthe severity of the third-order damping that was sostrong in the no-proof-mass case. This is modeled asamplitude-dependent damping through a higher-orderterm proportional to mc. The fifth-order influences areexpounded for the PZT-5A cantilever, which alreadyexperiences a stronger third-order nonlinearity. Figure7 mirrors the third-order to fifth-order model compari-son but for a different set of excitation values. This isbecause for a = 0:972 and Z .1:2g, we found that thedata could only be well characterized by a more com-plex seventh-order dynamical system. Again, due tothe continued low electric field generation, simulationresults encompassing nonlinear coupling, dielectriclosses, and electromechanical losses could not ade-quately match the empirical trends. Hence, seventh-order and higher harmonics in the frequency analysisof the time series data are suspected to also be due tomaterial elasticity. The threshold for which the cubicmodel breakdown occurs with the inclusion of aproof mass for PZT-5A is lower in comparison toPZT-5H. Figure 7(b) indicates that the cubic modelhas a much more limited range of validity as a

increases. Symmetry breaking and longitudinalstretching due to the tip mass orientation is not sus-pected as a cause for the breakdown of the fifth-ordermodel. Quadratic nonlinearities would only furthersoften the response (Nayfeh and Mook, 1979), whichis counter to the empirical results, and the heaviestproof mass did not induce initial curvature or astrong zeroth harmonic in the experimental timeseries. However, nonideal clamp conditions couldalso induce the increasing shift back toward linearresonance by generating a potential function with aweak dead-zone nonlinearity. To rule out hardeningdue to geometric effects, the equations of motionwere rederived with nonlinear Euler–Bernoulli beamtheory only to find that the nonlinear strain assump-tion introduced nonlinear terms of an order of magni-tude too small to be of any consequence.

One issue that arose in the course of this investiga-tion, even for the highly linear oscillations, is that non-linear damping parameters were not global and aremoderately influenced by the addition of a proof mass.The best-fit dimensionless mechanical damping con-stants for all cases are listed in Table 3. While proofmasses are known to require modification of lineardamping parameters, a more systematic approach fordetermining nonlinear damping in piezoelectric energyharvesters is a topic for further research in the future.

192 Journal of Intelligent Material Systems and Structures 23(2)

This would be in the same vein as the damping analysisin the study by Kim et al. (2010).

Linear versus nonlinear modeling

Guidelines for which simpler linear modeling is suffi-cient before nonlinear effects can no longer be ignoredare desirable. Hence, this section develops a quantita-tive metric for which application of a nonlinear modelis more appropriate for predicting harvester perfor-mance. The approach suggested herein studies the con-tributions of the lowest order nonlinear terms in theperturbation solution, which in this case are cubic, todetermine a point for which seeking a simple linear

model becomes a poor choice. We choose to disregardthe influence of gravity, provided the lack of inducedcurvature wrought by the addition of a proof mass,and to eliminate coupling terms in the perturbationsolution. Accordingly, we begin by estimating the con-dition for which the third-order nonlinearities will be ofthe same order as the linear terms. In view of equations(42a) and (42b), this occurs when

1

2Prj ja =

1

8Qrj ja3 ð49Þ

which translates to a critical dimensionless displace-ment a0 and voltage v0 defined by

a20 = 4

Pr

Qr

�������� ð50aÞ

and

v20 = 4

u2

1 + md

Pr

Qr

�������� ð50bÞ

500 520 540 560 5800

1

2

3

4

5

6

7

8

9

VR

(V)

Excitation Frequency (Hz)

500 520 540 560 5800

5

10

15

20

25

30

35

40

45

50

55

Dis

plac

emen

t(µ

m)

Excitation Frequency (Hz)

500 520 540 560 580 6000

1

2

3

4

5

6

7

8

9

VR

(V)

Excitation Frequency (Hz)

500 520 540 560 580 6000

5

10

15

20

25

30

35

40

45

50

55

Dis

plac

emen

t(µ

m)

Excitation Frequency (Hz)

(a)

(c) (d)

(b)

Figure 5. Experimental data points (circles), identified nonlinear resonance curves (black lines), and the linear model predictions(gray lines) for (a and c) PZT-5H and (b and d) PZT-5A bimorph with no attached proof mass. The excitation levels shown for PZT-5Hare 0.3 g, 0.56 g, 1.12 g, and 1.6 g, and for PZT-5A, the excitation levels are 0.2 g, 0.4 g, 0.8 g, 1.2 g, and 1.6 g.

Table 2. Identified higher-order elasticity coefficients for PZT-5H and PZT-5A.

PZT-5H PZT-5A

c4 �9:608631017 �9:772731017

c6 9:695031025 1:470031026

Stanton et al. 193

For PZT-5H, with no proof mass excited near reso-nance (Oe = 540Hz), solving equation (50a) givesa0 = 0:3199 and equation (50b) yields v0 = 0:0887.However, as shown in Figure 8, this value is well abovethe point from which the more accurate nonlinear solu-tion begins to deviate from the linear solution. Thismeans that while the nonlinear and linear terms are notof the same order, nonlinear effects can still be invigo-rated before a0, and hence, v0 reaches the critical value.

The next step toward guiding application of the non-linear nonconservative model is to define a suitablemeasure for which the linear model’s overpredictionbecomes unacceptable. Since the generated power is ofprimary importance in evaluating a harvester’s perfor-mance, we calculated the error in generated powerbetween the linear and nonlinear models according to

EP =jPlin � Pnonlinj

Plinð51Þ

where Plin and Pnonlin denote the dimensional powerdissipated by the resistive load. Figure 9 gives the trendin percent error between the linear and fully nonlinearmodels for the PZT-5H beam with no proof mass as afunction of the dimensionless drive and at the samedriving frequency as before (540 Hz). The 10% and

20% error thresholds are crossed when fz = 0:0054 andfz = 0:0074, respectively. The corresponding nonlinearamplitude response as calculated from equation (44)indicates that a minimum 10% error in generatedpower occurs for acrit0 � 0:123 and vcrit0 � 0:034, while aminimum 20% error occurs for acrit0 � 0:159 andvcrit0 � 0:044. Considering the conditions for which bothlinear and nonlinear terms in the multiple scales modu-lation equations were of the same order (equations(50a) and (50b)), the 10% error threshold is approxi-mately crossed when

acrit5H, 10’

3

4

Pr

Qr

��������

1=2

and vcrit5H, 10’3

4

uffiffiffiffiffiffiffiffiffiffiffiffiffi1 + md

p !

Pr

Qr

��������

1=2

ð52aÞ

while the 20% error threshold is given by

acrit5H, 20’

Pr

Qr

��������

1=2

and vcrit5H, 20’uffiffiffiffiffiffiffiffiffiffiffiffiffi

1 + md

p !

Pr

Qr

��������

1=2

ð52bÞ

Figure 10 verifies an expectation that lower excita-tion thresholds will accumulate larger error due to thestronger softening effect of PZT-5A. Using the same

200 210 220 230 240 250 260 2700

2

4

6

8

10

12

14

VR

(V)

Excitation Frequency (Hz)320 340 360 380 400

0

2

4

6

8

10

12

14V

R(V

)

Excitation Frequency (Hz)

(a) (b)

Figure 6. Experimental data points (circles), identified nonlinear resonance curves with fifth-order nonlinearities (black lines), andthe cubic model predictions (gray lines) for the PZT-5H bimorph with (a) Mt = 0:24 kg and (b) Mt = 0:787 kg. The excitation levelsshown are 0.145 g, 0.31 g, 0.56 g, 1.12 g, and 1.7 g.

194 Journal of Intelligent Material Systems and Structures 23(2)

reasoning for a mathematical development as before,we find that for PZT-5A, the 10% error threshold isapproximately crossed when

acrit5A, 10’

2

3

Pr

Qr

��������

1=2

and vcrit5A, 10’2

3

uffiffiffiffiffiffiffiffiffiffiffiffiffi1 + md

p !

Pr

Qr

��������

1=2

ð53aÞ

while the 20% error threshold, similar to that for PZT-5H, is found to also be scaled by one

acrit5A, 20’ j Pr

Qr

j 1=2

and vcrit5A, 20’uffiffiffiffiffiffiffiffiffiffiffiffiffi

1 + md

p !

j Pr

Qr

j 1=2

ð53bÞ

The addition of a proof mass changes the aboveguidelines. For example, the scaling factor for the 10%error threshold for the PZT-5A beam with the heaviestproof mass was found to be about 6=5 as opposed to2=3 for the no-proof-mass case. Similarly, the second

320 340 360 380 4000

1

2

3

4

5

6

7

8

9

10V

R(V

)

Excitation Frequency (Hz)200 210 220 230 240 250 260 2700

1

2

3

4

5

6

7

8

9

10

VR

(V)

Excitation Frequency (Hz)

(a) (b)

Figure 7. Experimental data points (circles), identified nonlinear resonance curves with fifth-order nonlinearities (black lines), andthe cubic model predictions (gray lines) for the PZT-5A bimorph with (a) Mt = 0:24 kg and (b) Mt = 0:795 kg. The excitation levels in(a) are 0.4 g, 0.8 g, 1.2 g, 1.6 g, and 2 g, while in (b), we show 0.3 g, 0.4 g, 0.8 g, and 1.2 g.

0 0.002 0.004 0.006 0.008 0.01 0.0120

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

fz

vss

Figure 8. Experimental data points (red circles), nonlinear pre-diction (black line), and linear theory (gray line) for the voltageresponse at O= 540Hz for the PZT-5H bimorph with no proofmass.

Table 3. Identified dimensionless damping parameters for all sixcantilever configurations.

Cantilever configuration ma mb mc

PZT-5H, no proof mass 0.0086 0.895 22.3PZT-5H, proof mass 1 0.0052 1.149 23.5PZT-5H, proof mass 2 0.008 0.977 22.7PZT-5A, no proof mass 0.0088 1.021 0PZT-5A, proof mass 1 0.0114 0.893 0PZT-5A, proof mass 2 0.0125 0.453 0.1

Stanton et al. 195

proof mass attached to the PZT-5H results in 10%error for a scaling factor of 9=10 instead of 3=4 for theno-proof-mass case. This further supports the experi-mental and theoretical findings that a proof mass canhave a strong influence on the ensuing nonlinearresponse of a piezoelectric beam. As may also beexpected, the change in scaling factor for PZT-5A is fargreater than for PZT-5H due to the larger materialsoftness.

To better visualize the influence of a proof mass onthe nonlinear dynamics, Figures 11 and 12 show howboth the base acceleration and proof mass influence theerror in generated power from the linear model. In bothgraphs, there is an approximately linear relationshipbetween the harvester-mass ratio and base acceleration

and the resulting linear model error in generated power.Figures 11 and 12 also show the significantly strongernonlinear response characteristics of PZT-5A, wherethe error thresholds between the more accurate non-linear model and the standard linear model are muchlower in comparison to PZT-5H.

These figures are of particular value for definingparameter regions for which a linear model is sufficientor when nonlinear modeling is more appropriate. ForEP < 10%, a linear model will provide a simple determi-nation of the harvester response and is thus adequate.Depending on one’s willingness to accept increasingerror, the appropriateness of a linear approach willvary. For a PZT-5H cantilever with no proof mass, alinear model will be sufficient up to about Z = 0:8 g,whereas a linear model for the same cantilever consistedof PZT-5A will become unsatisfactory when Z = 0:5 g.In all cases, successively heavier proof masses introducepiezoelectric nonlinearity and greatly reduce the para-meter space for which the validity of a linear model isan acceptable modeling assumption. The linear rela-tionships were determined by a regression analysis andare of the form

Z = p1a + p0 ð54Þ

where a is the mass ratio from equation (47) and Z isthe base acceleration. Numerical parameters for p1 andp0 are listed in Tables 4 and 5. For a given base accel-eration Z, mass ratio a, and error tolerance EP, theregression lines provide a guide for which linear piezo-electric constitutive relations will accurately model theharvester response.

Summary and conclusions

This article thoroughly investigated the influence oftwo types of common piezoelectric materials on theresponse of prototypical energy-harvesting beams. Thedevices were excited from low to high g base accelera-tion and for two different proof masses. Nonlinearresonance curves were observed and modeled from thefirst principles with a dissipation formulation due toviscoelectroelasticity. The Rayleigh–Ritz method wasemployed to directly derive a system of ordinary differ-ential equations through Hamilton’s extended princi-ple. The dynamical system was solved analytically bythe method of multiple scales, and the ensuing non-linear frequency response function was used to providea basis for parameter identification in conjunction withsteady-state laboratory measurements. Following theanalysis of description of the experiment and findings,a framework for determining when nonlinear modelingbecomes necessary was provided.

Several important observations and conclusions aresummarized as follows:

0 0.002 0.004 0.006 0.008 0.01 0.012 0.0140

5

10

15

20

25

30

35

40

fz

E P

PZT−5A, proof mass 2PZT−5A, no proof massPZT−5A, proof mass 1

Figure 10. Error in predicted power due to linear piezoelectri-city assumptions for the PZT-5A bimorph.

0 0.002 0.004 0.006 0.008 0.01 0.012 0.0140

5

10

15

20

25

30

35

40

fz

E P

PZT−5H, no proof massPZT−5H, proof mass 1PZT−5H, proof mass 2

Figure 9. Error in predicted power due to linear piezoelectri-city assumptions for the PZT-5H bimorph.

196 Journal of Intelligent Material Systems and Structures 23(2)

1. From a kinematic perspective, linear Euler–Bernoulli beam theory is appropriate for ourexperiments and for a wide range of bimorphharvesters, especially if uniformly laminated. Athird-order strain assumption along with aninextensibility condition (Da Silva and Glynn,1978; Stanton and Mann) was modeled only tofind that geometric nonlinearities were too smallto have any influence.

2. Nonlinear damping is critical and can be wellexplained and modeled by nonconservativepiezoelectric constitutive relations. However,more research is required to determine the inter-play between material losses and fluid-structureinteraction. This is an important undertaking inview of the fact that nonlinear damping wasobserved to be reducing the response amplitudeat resonance before cubic nonlinearities begin toshift the resonance curves to the left.

3. Elasticity is the primary source of nonlinearity.That is, higher-order coupling and dielectric effectscan be neglected for such weak electric fields.

4. With no proof mass, nonlinear behavior is suffi-ciently modeled by a cubic nonlinearity, third-order damping, and linear coupling. Even witha proof mass, this model is sufficient for lowbase accelerations. However, higher drivingamplitudes render fifth-order nonlinearities

important. This is especially true for PZT-5Aand for heavier proof masses.

Since most harvesters incorporate a proof mass forresonant frequency tunability, the results of this investiga-tion and the scaling analysis of the section titled ‘‘Linearversus nonlinear modeling’’ indicate that nonlinear elasti-city and damping become increasingly important. Thedissipation model utilized in this study was selected overair drag in view of the nonlinear resonance curves of

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Acc

eler

atio

n(g

)

EP ≤ 40%EP ≤ 30%EP ≤ 20%EP ≤ 10%

Mt/mL

Figure 12. Influence of the mass ratio a and base accelerationon the linear model error for the predicted power of the PZT-5A bimorph.

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Mt/mL

Acc

eler

atio

n(g

)

EP ≤ 40%

EP ≤ 30%

EP ≤ 20%

EP ≤ 10%

Figure 11. Influence of the mass ratio a and base accelerationon the linear model error for the predicted power of the PZT-5H bimorph.

Table 4. Regression coefficients for PZT-5H.

Regressioncoefficients

EP < 10% EP < 20% EP < 30% EP < 40%

p1 24.1707 25.2880 26.4307 27.7441p0 7.9644 10.9745 13.9298 17.4131

Table 5. Regression coefficients for PZT-5A.

Regressioncoefficients

EP < 10% EP < 20% EP < 30% EP < 40%

p1 22.3145 22.9230 23.5569 24.2415p0 4.9598 6.0000 7.0836 8.2538

Stanton et al. 197

piezoelectric oscillators not subject to large deformations(Samal et al., 2006b; Von Wagner, 2003). However, thehigh oscillation velocities may be rendering both quadra-tic air damping and nonlinear structural dissipation rele-vant. Correctly determining the physical mechanism ofthe harvester dissipation aside from the harvesting pro-cess is important for accurate modeling as illustrated inour results. Methods within several articles concerningatomic force microscopes immersed in fluids may be ofbenefit in investigating and modeling damping in piezo-electric energy harvesters (Chon et al., 2000; Cole andClark; Sader, 1998) along with the recent work of Kamelet al. (2010). Furthermore, the presumption of linear pro-portional damping in a distributed parameter systemitself, though the widespread norm, may require reexami-nation (Banks et al., 1998).

Future work is also in order pertaining to the non-linear coupling. While this article only examined oneimpedance load such that low voltages were generated,larger impedance loads translate to greater electricpotential fields in the electroelastic laminates and maythus render higher-order electrical effects and nonlinearcoupling important. Hence, another direction for futureresearch would be to investigate the influence of resis-tive loads on the nonlinear response of piezoelectric har-vesters to determine at what point nonlinear couplingbegins to influence the nonlinear resonance curves.

Funding

Financial support from the Duke University Center forTheoretical and Mathematical Sciences and Dr Ronald Joslinthrough an Office of Naval Research (ONR) YoungInvestigator Award is gratefully acknowledged as is supportfrom AFOSR MURI Grant No. F-9550-06-1-0326, ‘‘EnergyHarvesting and Storage Systems for Future Air ForceVehicles,’’ monitored by Dr B.L. Lee.

References

Banks HT, Luo Z, Bergman LA, et al. (1998) On the existence

of normal modes of damped discrete-continuous systems.

Journal of Applied Mechanics 65(4): 980–990.Chon J, Mulvaney P and Sader J (2000) Experimental valida-

tion of theoretical models for the frequency response of

atomic force microscope cantilever beams immersed in

fluids. Journal of Applied Physics 87(8): 3978.Cole DG and Clark RL. (2007) Fluid-structure interaction in

atomic force microscope cantilever dynamics and thermal

response, J. Appl. Phys. 101, 034303.Crandall S, Karnopp D, Kurtz E, et al. (1968) Dynamics of

Mechanical and Electromechanical Systems. Malabar, FL:

Krieger.Crawley E and Anderson E (1990) Detailed models for the

piezoceramic actuation of beams. Journal of Intelligent

Material Systems and Structures 1(4): 4–25.Daqaq M, Stabler C, Qaroush Y, et al. (2009) Investigation

of power harvesting via parametric excitations. Journal of

Intelligent Material Systems and Structures 20(5): 545–557.

Da Silva MC and Glynn C (1978) Nonlinear flexural-flex-ural-torsional dynamics of inextensional beams: I. Equa-

tions of motion. Mechanics Based Design of Structures and

Machines 6(4): 437–448.DuToit N, Wardle B and Kim S (2005) Design considerations

for MEMS-scale piezoelectric mechanical vibration energyharvesters. Integrated Ferroelectrics 71: 121–160.

Erturk A and Inman D (2008a) A distributed parameter

electromechanical model for cantilevered piezoelectric energyharvesters. Journal of Vibration and Acoustics 130: 1–15.

Erturk A and Inman D (2008b) Comment on modeling and

analysis of a bimorph piezoelectric cantilever beam forvoltage generation. Smart Material Structures 17(5),doi:10.1088/0964-1726/17/5/058001

Erturk A and Inman D (2008c) Issues in mathematical mod-

eling of piezoelectric energy harvesters. Smart Material

Structures 17(6): 1–14.Erturk A and Inman D (2009a) An experimentally validated

bimorph cantilever model for piezoelectric energy harvestingfrom base excitations. Smart Material Structures 18: 1–18.

Erturk A and Inman D (2009b) On mechanical modeling of

cantilevered piezoelectric vibration energy harvesters. Jour-nal of Intelligent Material Systems and Structures 19: 1311–

1325.Hagood N, Chung W and Von Flotow A (1990) Modelling of

piezoelectric actuator dynamics for active structural con-trol. Journal of Intelligent Material Systems and Structures

1(2): 327–354.Hu Y, Xue H and Yang J (2006) Nonlinear behavior of a

piezoelectric power harvester near resonance. IEEE Trans-

actions on Ultrasonics, Ferroelectrics, and Frequency Con-

trol 53(7): 1387–1391.Ikeda T (1996) Fundamentals of Piezoelectricity. Oxford:

Oxford University Press.Jiang S, Li X, Guo S, et al. (2005) Performance of a piezoelec-

tric bimorph for scavenging vibration energy. Smart Mate-

rials and Structures 14(4): 769–774.Kamel T, Elfrink R, Renaud M, et al. (2010) Modeling and

characterization of MEMS-based piezoelectric harvesting

devices. Journal of Micromechanics and Microengineering

20(105023): 14.Kim JDM, Hoegen M and Wardle BL (2010) Modeling and

experimental verification of proof mass effects on vibration

energy harvester performance. Smart Material Structures

19(4): 045023.Masana R and Daqaq M (2011) Electromechanical modeling

and nonlinear analysis of axially-loaded energy harvesters.Journal of Vibration and Acoustics 133(1): 011007.

Maugin G (1985) Nonlinear Electromechanical Effects and

their Application. World Scientific Publishing Co. Pte.Ltd., Singapore.

Nayfeh A and Mook D (1979) Nonlinear Oscillations. New

York: Wiley.Parashar S and Von Wagner U (2004) Nonlinear longitudinal

vibrations of transversally polarized piezoceramics: experi-ments and modeling. Nonlinear Dynamics 37: 51–73.

Priya S, Viehland D, Carazo A, et al. (2001) High-power reso-

nant measurements of piezoelectric materials: importanceof elastic nonlinearities. Journal of Applied Physics 90(3):

1469–1479.Sader J (1998) Frequency response of cantilever beams

immersed in viscous fluids with applications to the

198 Journal of Intelligent Material Systems and Structures 23(2)

atomic force microscope. Journal of Applied Physics 84(1):64–76.

Samal M, Seshu P, Parashar S, et al. (2006a) Nonlinear beha-viour of piezoceramics under weak electric fields: Part-i: 3-d finite element formulation. International Journal of Solidsand Structures 43(6): 1422–1436.

Samal M, Seshu P, Parashar S, et al. (2006b) Nonlinear beha-viour of piezoceramics under weak electric fields. Part-ii:numerical results and validation with experiments. Interna-tional Journal of Solids and Structures 43(6): 1437–1458.

Smits JG, Dalke SI and Cooney TK (1991) The constituentequations of piezoelectric bimorphs. Sensors and Actuators

A: Physical 28(1): 41–61.Sodano H, Park G and Inman D (2004) Estimation of electric

charge output for piezoelectric energy harvesting. Strain40: 49–58.

Stanton S and Mann B Nonlinear electromechanical

dynamics of piezoelectric inertial generators: Modeling,analysis, and experiment. Nonlinear Dynamics.

Stanton S, Erturk A, Mann B, et al. (2010a) Nonlinear piezo-electricity in electroelastic energy harvesters: modeling andexperimental identification. Journal of Applied Physics

108(1): 074903.Stanton S, Erturk A, Mann B, et al. (2010b) Resonant mani-

festation of intrinsic nonlinearity in electroelastic micro-power generators. Applied Physics Letters 97: 254101.

Tiersten H (1993) Electroelastic equations for electroded thinplates subject to large driving voltages. Journal of AppliedPhysics 74(5): 3389–3393.

Triplett A and Quinn D (2009) The effect of non-linear piezo-electric coupling on vibration-based energy harvesting.Journal of Intelligent Material Systems and Structures

20(16): 1959–1967.Von Wagner U (2003) Non-linear longitudinal vibrations of

piezoceramics excited by weak electric fields. InternationalJournal of Nonlinear Mechanics 38(4): 565–574.

Von Wagner U and Hagendorn P (2002) Piezo-beam systemssubjected to weak electric field: experiments and modelingof nonlinearities. Journal of Sound and Vibration 256(5):861–872.

Yang J (2005) Analysis of Piezoelectric Devices. New York:Springer.

Yu S, He S and Li W (2010) Theoretical and experimentalstudies of beam bimorph piezoelectric power harvesters.Journal of Mechanics and Materials of Structures 5(3):427–445.

Zhang W and Meng G (2005) Nonlinear dynamical system ofmicro-cantilever under combined parametric and forcing

excitations in MEMS. Sensors and Actuators A: Physical

119(2): 291–299.

Appendix 1

Definition of modal coefficients

M = b rbhb + 2rphp

� �b

ðL0

f2 dx + Mtf Lð Þ2 + Jof9 Lð Þ2

ðA:1Þ

K = cxx, bIb + 2cExxIp

� � ðL0

f002 dx ðA:2Þ

G =1

240cE

4 bhp 16h4p + 40h3

phs + 40h2ph2

s + 20hph3s + 5h4

s

� �

3

ðL0

f4 dx ðA:3Þ

Y=1

2~Y=

1

2be31 hp + hs

� �f9 Lð Þ ðA:4Þ

C = h�1p

~C = eSZZ

bL

2hp

ðA:5Þ

Fg = g rbhb + 2rphp

� �b

ðL0

f dx + Mtf Lð Þ

24

35 ðA:6Þ

Da = hbIb + 2hpIp

� � ðL0

f002 dx ðA:7Þ

Db =3

40zbhp 16h4

p + 40h3phs + 40h2

ph2s + 20hph3

s + 5h4s

� � ðL0

f004dx

ðA:8Þ

Dd =1

2~L =

1

2gb hp + hs

� �f9 Lð Þ ðA:9Þ

De =1

2x = n

bL

2hp

ðA:10Þ

Stanton et al. 199