Nonlinear vibration analysis of a circular composite plate harvester via harmonic balance

Tian-Chen Yuan·Jian Yang·Li-Qun Chen,3,4

Abstract A lumped parameter transverse vibration model of a composite plate harvester is analyzed via harmonic balance approaches.The harvester is mainly composed of a piezoelectric circular composite clamped by two steel rings and a proof mass on the plate.The lumped parameter model is a 1.5 degree- of-freedom strongly nonlinear system with a higher order polynomial stiffness.A harmonic balance approach is developed to analyze the system, and the resulting algebraic equations are numerically solved by adopting an arc-length continuation technique.An incremental harmonic balance approach is also developed for the lumped parameter model.The two approaches yield the same results.The amplitude-frequency responses produced by the harmonic balance approach are validated by the numerical integrations and the experimental data.The investigation reveals that there coexist hardening and softening characteristics in the amplitude-frequency response curves under sufficiently large excitations.The harvester with the coexistence of hardening and softening nonlinearities can outperform not only linear energy harvesters,but also typical hardening nonlinear energy harvesters.

Keywords Piezoelectric energy harvester·Circular composite plate·Transverse nonlinear vibration·Harmonic balance·Arc-length continuation

1 Introduction

Vibration-based energy harvesting can capture waste kinetic energy from the ambient environment to produce electricity.The power produced can be used to run low-power consumption appliances[1,2].Therefore vibratory energy harvesting is a promising research field of great interest.

Linear vibratory energy harvesters work in a narrow frequency range.It is a significant is sue to increase the working frequency range[3,4].A potential approach is to introduce appropriate nonlinearities into harvesting systems in order to broaden the response frequency b and width[5,6].Many researchers have explored the applications of nonlinearities in vibratory energy harvesting,for example,by adding magnets[7–9],establishing multi-stable systems[10–12] and designing special shapes[13–15]so that nonlinear restoring forces are produced.In addition to multi-stable systems,a typical nonlinear behavior,hardening or softening the nonlinear characteristic of the amplitude-frequency response curve,is widely adopted to enlarge harvester working frequency b and s.The study by Tang and Yang[8]designed a harvester with a magnetic oscillator and found numerically and experimentally a softening nonlinearity type.In their work,Leadenham and Erturk[13]proposed an M-shaped asymmetric nonlinear oscillator and predicted the bending of hardening nonlinearity via a harmonic balance method to with experimental validation.Furthermore,Chen et al.[16–18]proposed internal resonances as a novel possibility to improve energy harvesting.The double-bending of the amplitude-frequency was predicted by the method of multiple scales[16,18] and validated by time-domain simulations[16] and experimental data.Work by Tvedt et al.[19]investigated an electrostatic vibration-based energy harvester that exhibits rich nonlinear behaviors.Numerical simulations and experimental results showed that the softening nonlinearity and the hardening nonlinearity coexist.

Design of appropriate structures for energy harvesting is another way to enhance the performance of vibratory energy harvesters.The most used structure of piezoelectric energy harvesters is a cantilever beam with piezoelectric ceramic layers.Many investigations are devoted to the design,modeling and experiments of cantilever beam harvesters[20,21].Circular piezoelectric plate structures are possibly effective for energy harvesting as their coupling coefficient is high and structures are equally compact in all directions in plate planes.Since such axial symmetric structure has the advantages of simplicity,low cost and easy machining,it attracts increasing attention[22–29].Investigations by Yuan et al.[24] and Wang et al.[25]theoretically and experimentally studied drum-like harvesters consisting of two circular piezoelectric plates.Also,Kim et al.[26,27]studied a circular plate for piezoelectric harvesting energy from pressures.The work of Chen et al.[28]presented a circular piezoelectric plate with a center mass as a harvester in order to reduce the fundamental frequency of the circular plate structure.In their study,Zhao et al.[29]proposed a harvesting array with four different frequency linear piezoelectric circular diaphragms and implemented broad frequency vibratory energy harvesting to overcome the narrow b and width of linear systems.Because of the brittleness of isotropic piezoelectric plates,a composite plate produced by a metal base was adopted in the above designs.In addition,a center mass is needed to reduce natural frequencies of the plate harvester to resonate with environmental vibrations.To this end,Yang et al.[30]studied linear vibration based on a partial differential equation of a pro of-mass circular diaphragm harvester with a center mass.A study by Xue et al.[31]explored nonlinear behaviors of a circular isotropic piezoelectric plate without any center masses and found a hardening characteristic due to the von Karman geometric nonlinearity.Work by Hedrih and Simonovic[32]investigated a circular sandwich plate with an interconnecting viscoelastic layer and found one or more jumps with softening characteristics occurring in its amplitude-frequency response curves.Previous works[33,34]experimentally and numerically revealed softening and hardening characteristics of a center-massed circular composite plate harvester under larger excitations.Such phenomenon can serve as a mean to broaden the b and width more effectively than typical hardening nonlinear systems.

The previous works[33,34]presented only experimental and numerical results without any analytical predictions.It is a challenging task to predict analytically dynamic responses of nonlinear vibratory energy harvesters.A few approximate analytical methods have been used to analyze vibratory energy harvesters,for example,the method of multiple scales[16–18] and the method of harmonic balance[13,35,36].Although the method of harmonic balance has been successfully applied to determine steady-state responses of vibration energy harvesters[13,36],it is still difficult to deal analytically with higher order harmonic terms because of the strongnonlinearity[33].The method of incremental harmonic balance[37,38]is also a powerful approach to analyze a system with complex and strong nonlinearities[39].In practice,coefficients of Fourier series can be approximated by numerical methods,for example,the numerical integrations[40]or the fast Fourier transforms[41].The New ton–Raphs on iteration was widely used in harmonic balance(HB) and incremental harmonic balance(IHB)methods.Never the less,the iteration is failed at turning points where the Jacobian matrix .It is also difficult to pass the turning point by using the frequency continuation only.The arc-length continuation technique can be applied to overcome the defects[42].

The work of Yuan et al.[43]developed a complicated method to calculated multi-degrees- of-freedom piezoelectric mechanical systems.However,it is overkill to investigate a low dimension system,such as the circular composite plate harvester.In this work,the method of HB incorporated by the arc-length continuation technique is developed to analyze a piezoelectric energy harvester.The IHB method is also developed and the results via the two approaches are compared.The complete strategies of the proposed approach include an HB progress,an arc-length continuation and a stable analysis.The developed methods are more succinct and easier to be applied than that in Ref.[43].The resulting amplitude-frequency responses of the harvester are validated by numerical iterations and experimental data.

The manuscript consists of six sections.Section 2 describes a mathematical model of a circular composite plate harvester.Section 3 proposes HB and IHB methods implemented with the arc-length continuation technique.Section 4 determines the amplitude-frequency curves in different excitation levels and compares the approximate analytical solutions with experimental results and time-domain simulations.Section 5 highlights the effects of higher order harmonic and discusses the voltage responses under different parameters.The performance of the harvester in different resistance is also demonstrated.Section 6 presents the conclusions.

2 Lumped parameter model with identified parameters

Figure 1 shows a pro of-mass circular composite plate harvester.The dimensions are proportionally scaled.The composite plate consists of a brass base plate with radius Rb?40 mm and thickness hb?0.2 mm and two identical PZT-5H piezoelectric diaphragms with radius Rp?20mm and thickness hb?0.2 mm.An inverted cone shaped mass is fixed at the composite plate center.Two steel circular rings with radius Rs?22.5 mm are used to clamp the composite plate.The rings are excited by a shaker.The two piezoelectric elements are connected in series.The upper and the lower surfaces of the two piezoelectric elements are the silver electrodes.

Fig.1 Pro of-mass circular composite plate harvester

Fig.2 Lumped parameter model

If only transverse displacement is considered,one can establish a lumped parameter model as a forced mass damper-spring oscillator with an additional electric equation(Fig.2).The dynamic equations of the model are

where m,c and η denote,respectively,the equivalent mass,the linear damping coefficient and the electromechanical coupling coefficient.Cpis the total capacitance of the two piezoelectric elements.z?za?zbis the displacement of the center mass(with absolute displacement za)relative to the moving base(with displacement zb) and u is the output voltage of the harvester.These parameters can be experimentally identified.The governing equations can be obtained by discretizing vibrating plate equations into a set of ordinary differential equations.However,the expression of the nonlinear stiffness function kn(z)is hard to obtain because of the complicated nonlinear factors such as non-ideal gluing,physical uncertainties of materials and large deformations.However,the nonlinear restoring force fs(z)?kn(z)z can be determined in further experiments thought system identification progress.

The electrical parameters η and Cpare identified from the experimental data under harmonic excitations via the leastsquares method[33].The parameters are identified as m?0.1209 kg,c?6.28 N·s/m,η?0.0108 N/V, and Cp?65.5 nF.The nonlinear restoring force can be approximated by[33]

with the coefficients identified as k1?3.61×104N/m,k2?1.9×107N/m2,k3??10.8×1010N/m3,k4??8.1×1012N/m4, and k5?3.5×1017N/m5.

3 Harmonic balance analysis implemented via an arc-length continuation technique

3.1 Harmonic balance process

Assume the harmonic base excitation takes the form

where A is the excitation amplitude and ω is the driving frequency.The vibration response can be expanded into Fourier series.Because the restoring force Eq.(3)has quadratic and quartic nonlinearity terms,a direct current(DC)component and even-order harmonics are added in the solutions[17].The displacement and voltage are the n

In order to examine the effects of higher order harmonics in the response with comparing with the experimental results and numerical integrations(ode45 in MATLAB),the higher harmonics are adopted in the further analysis.The resulting nonlinear algebraic equations can be cast into the form where x?[a10,a11,…,a1k,b11,…,b1k,a20,a21,…,a2k,b21,…,b2k],ω is a independent parameter,such as the driving frequency.As an exact solution is not available,a predictor—corrector method is used to solve numerically the unknowns.Differentiating Eq.(7)leads to

The New ton–Raphson method can yield the correction solution by using the predictive solution as the initial value.Choosing ω as a continuation parameter(“frequency continuation”)fails at turning points due to the singularity of

Fx(det[Fx]?0).Choosing one of the variables in x as a continuation parameter(“displacement continuation”)may be possible to pass the turning point.However,it may fail at another point when the new Jacobian matrix is singularity.A more general approach termed the arc-length continuation technique can overcome the limitations, and it will be presented in details in Sect.3.2.

3.2 Arc-length continuation technique

Unlike the frequency continuation or the displacement continuation,arc-length s along the solution curve c(s)is chosen as an independent parameter.Let Y?(x,ω) and ignore the difference between x and ω,Eq.(7)can be rewritten as

Suppose that J?[J1,J2,…,Jn+1]Tis the tangent vector of the solution curve.The element in vector J can be calculated by

where D?[Fx,Fω].The tangent vector is normalized by τ?J/||J||, and τ is the unit tangent vector.In practice,the forward Euler method is used to give predictive solution, and the calculation formulation is

where Ypis the predictive value,Y0?[x0,ω0]Tis the initial value, and ?s is variable increment of the arc-length.The predictive value can be corrected by the New ton–Raphson method

The magnitude of increment length?s can be chosen with the curvature of the solution curves.The larger the curvature is,the smaller the increment length should be.Actually,a rather simple formula can be adopted[37],?sj??sj?1Nd/Ij?1,where?sj?1is the increment length at the (j?1)-th prediction step,Ij?1is the iteration number at the (j?1)-th correction step, and Ndis the desired iteration number.As demonstrated in the next section,the technique works in passing turning points and following a solution curve.

3.3 Stability analysis

The harmonic balance method yields both stable and unstable periodic solutions.Therefore,it is necessary to analyze the stability of the harmonic balance solutions.The Floquet theory can beused to determine the stability of periodic solutions of the nonlinear energy harvester system.Introduce U(t)?[z,,u],transform governing Eqs.(1) and (2)into

and superpose?U to perturb the assumed periodic solution U*(t) of the system.The stability can be examined by the linear stability of

Hsu’s method[44]approximates the monodromy matrix as

where I denotes an identity matrix.Solution period T is divided into Nssubinterval?T.Njdenotes the number of terms in the approximation of Hnexponential.Constant matrix Hn?H[U*(tn)]substitutes the time-varying matrix H[U*(t)]in the n-th time interval,where tn?n?T/Ns.

If all the eigenvalues of H(Floquet multipliers)are in a unit circle on the complex plane,the periodic solutions are stable.Otherwise,the periodic solutions are unstable.The solving steps for the method of harmonic balance with the arc-length continuation technique are surveyed in Fig.3.

Fig.3 flow chart of the method of harmonic balance with the arc-length continuation

Step 1 Prepare matrix F and augmented Jacobian matrix D;

Step 2 Suppose a initial value Y0?[x0,ω0]Tand determine the algorithm parameters Nd,?s and ing;

Step 3 Use Eq.(12)to calculate predictive solution Yp;

Step 4 Use Eq.(13)to calculate the correction solution

Step 5 Analyze stability by Eq.(16) and record the solution

Step 6 Update the ?s and Y0and repeat Step 2 to Step 5 until all the iteration is finished.

3.4 Incremental harmonic balance process

Replace time variable t by a new variables τ? ωt.Equations(1) and (2)become

with

Supposing z0and u0are the approximations of a periodic solution at frequency ω0,the neighboring state can be obtained by adding the m the corresponding increments

Substituting Eq.(20)into Eqs.(17) and (18) and neglecting higher order terms,one gets the linearized incremental equation

with

The sums of harmonic terms are assumed for the harvester’s responses.The displacement and voltage are the n

with

The vectors of harmonic coefficients and their increments can be expressed as A?[A1,A2]Tand ?A?[?A1,?A2]T.Substituting Eqs.(25)–(31)into Eqs.(21) and (22) and performing the Galerk in procedure,a set of linear equations in terms of?A and ?ω can be obtained readily

with

In Eq.(32),Kmcknown as the Jacobi matrix and its has the same form as Fx.So the augmented Jacobian matrix D also can be expressed as D?[Kmc,Rmc].In Eq.(37),C p is the capacitance matrix.In Eq.(38),RLis the resistance matrix.The solutions of HB and IHB are compared in Fig.4.The external excitation is2g(g is gravitational acceleration) and driven frequency is 75 Hz.It is worth emphasizing that the matrix F in HB method has the same expression as the matrix R in IHB method.So the curve and scatters have good agreement with each other in Fig.4.

Fig.4 Comparisons of HB solutions and IHB solutions

Both the incremental harmonic balance and harmonic balance with the arc-length approach have been developed for a piezoelectric energy harvesting system.The novelty lies in the fact that dynamics equations of such systems include a first-order differential equation.That is,those systems have an additional half degree- of-freedom.The methodological frameworks of these two methods have been expressed in detail.The y provide two possible approaches to analyze electromechanical systems,including but not limited to piezoelectric energy harvesting systems with strong nonlinearity.

The incremental harmonic balance method and the harmonic balance Newton–Raphson method are equivalent, and it was proofed by Ferri[45].The incremental process and the harmonic balance process can beexchanged.So,both of the methods have same results as shown in Fig.4.The incremental harmonic balance method may have the advantage of a multi-degrees- of-freedom system because of its standardized steps.The present work develops the incremental harmonic balance method to a lumped parameter model of the energy harvester with a so-called half-degree- of-freedom.It is also found that two approaches yield same results as expected.In the forthcoming treatments,the harmonic balance method with the arc-length approach will be utilized to reveal the characteristic of the harvester.

4 Confirmations of harmonic balance solutions

4.1 Numerical verification

In order to verify the harmonic balance solutions,both numerical integrations and harmonic balance solutions(the first-order term only and up to the third-order term)are pre-sented at four harmonic base excitation amplitudes(0.3g,1g,1.5g and 2g).The load resistance is fixed at 25 k?.The voltage amplitude-frequency response curves can be calculated by the HB method and the numerical integrations,which are compared in Fig.5.The HB solutions are represented by orange thick line(N?1) and green thick line(N?3)while numerical results for the frequency sweep are plotted as the triangular scatters(upward) and the circular blue scatters(downward).The Fourier series solutions with even-order harmonics and constant term can accurately approximate the responses of the asymmetric circular composite plate harvester.

Fig.5 Comparisons of time-domain simulations and HB solutions under different excitation levels

The system exhibits linear behavior in the 0.3g excitation case and the black dashed lines in Fig.5a reflect the linear natural frequency(86.97 Hz).The simulations and HB solutions(N?1 and N?3)agree well at such low excitation.Under the 1g and the 1.5g excitations,voltage amplitude-frequency response curves have softening characteristic(bending to the left)with a slight jump phenomenon both in the HB and numerical solutions(Fig.5b,c).The resonant frequency is smaller than the natural frequency of the linear system.

Both the softening and the hardening characteristics are demonstrated in Fig.5d.That is,a slight jump due to the softening is observed on the left side of the resonance peak, and a jump due to the hardening on the right.It should be noted that the resonant frequency is also smaller than the linear one,which is different from the typical Duffing oscillator.Overall,the solution up to order 3(N?3)agrees with the numerical solution much better than the first order solution(N?1).

The possible reason for the softening nonlinearity is that the brass base plate and piezoelectric diaphragm are assembled together by the glue.The Young’s modulus of the glue is much lower than the brass and the piezoelectric ceramics.Under small excitations,the negative cubic term of Eq.(3)dominates the nonlinear effect and the softening occurs.For sufficiently large excitations,both the cubic and the positive quintic terms play apart, and thus both the softening and the hardening nonlinearities become visible.

The solutions up to order 5(N?5) and the solutions up to order3(N?3)are compared in Fig.6.The external excitation amplitude is 2g.The comparisons of turning points A,B,C and D in Fig.6 are listed in Table1.Two solutions are close to each other and agree well with the numerical results.So the solution up to order 3(N?3)is accuracy enough to meet the needs of practice.

4.2 Experimental validations

Upward and downward frequency sweep experiments are carried out.Figure 7 presents the voltage amplitude frequency response curves of both the experimental frequency sweep and the HB solutions(N?3).The experimental results are shown in the left column and HB results are represented in the right column.The HB solutions agree well with the experimental results which display both the softening and the hardening characteristics.The resonant frequencies of HB solutions under different excitation levels are very close to the experimental ones.Overall,the model-based HB solutions match well with the experimental results.

Fig.6 Comparisons of higher order harmonic solutions

Table1 Turning points ofhigher order harmonic solutions

Fig.7 Comparisons of experimental results and HB solutions under different excitation levels

5 Discussions based on harmonic balance solutions

5.1 Solution branch of higher harmonics

The effect of the higher order harmonics and the solution curves of the harmonic terms will be discussed in detail.Figures,,display the solution curves for harmonic coefficients under different excitations with the 25 k?load resistance via the method of harmonic balance with the arc-length continuation technique,which was presented in Sect.4.

Fig.8 Solution curves of the first order harmonic.a Voltage and b displacement

Fig.9 Solution curves of the second order harmonics.a Voltage,b displacement

Figures 8,9, and 10 are the first order,the second order, and the third order harmonic solution curves,respectively.Figures 8a,9a, and 10a are the voltage harmonic solutions curves, and Figs.8b,9b, and 10b are the displacement solution curves.The solid and dash space curves are the stable and the unstable parts of the solutions,respectively.The space solution curves are smooth and continuous without sharp point and knot phenomenon.However,the harmonic coefficient response curves(black lines)have sharp points and knots at the switch between the stable and the unstable solutions.The gray lines show the relationships between the two harmonic coefficients,aijand bij,for each order harmonic,which are the Nyquist diagrams.The Nyquist diagrams are a combination of single or several incomplete circles.It can be found that the higher order harmonic solutions have lower value and more complex form.Because of the coupling between the harmonics,ignoring these high harmonics may lead to the error shown in Fig.5.The method of harmonic balance with the arc-length continuation technique can be applied to pass the sharp points and knots and follow the solution curve of the higher order harmonics.

Figure 11 demonstrates the contributions of the first,the second, and the third harmonics to the displacement,the acceleration,the voltage and the power at2g amplitude acceleration level.The magnitude order of the first order harmonic is higher than all other orders components least one.However,these different harmonics are coupled in the nonlinear algebraic equations.Thus ignoring higher order harmonics may yield significant errors in the prediction of the response,such as the frequency and the magnitude in Fig.5.The DC component of the displacement has a higher magnitude than the second and the third order harmonics near the main resonance, and it may cause the DC shift of the displacement.Figure11expounded that the second or the third order superharmonic resonance arises when the excitation frequency approaches one half or one third of its liberalized natural frequency, and the y are caused by the quadratic or the cubic components in the restoring force.These super-harmonic resonances are more obvious in the acceleration-frequency curve(Fig.11c).

Fig.10 Solution curves of the third order harmonics.a Voltage,b displacement

Fig.11 Magnitudes of the harmonics at 2g excitation for different response.a Displacement,b voltage,c acceleration,d power

Fig.12 Comparisons for the voltage responses subjected to the restoring force parameter disturbance.All stiffness coefficients a increases5% and b decreases 5%

5.2 Robustness of the co-existence of the hardening and softening characteristics

To examine the robustness of the co-existence of hardening and softening characteristics,the voltage-frequency response will be calculated with the stiffness coefficients slightly different from original parameters.Figure 12a displays the voltage response for the stiffness coefficients increased by 5%.The linear resonance frequency also increases from 86.97 to 89.12 Hz.Figure12b displays the voltage response in the stiffness coefficients decreased by 5%.The linear resonance frequency also decreases from 86.97 to 84.77 Hz.The voltage-response curves have no qualitative differences in the parameter disturbance of the restoring force.

Figure 13 shows the voltage-frequency response for the damping coefficient and the mass different from Table 1.Figure 13a displays voltage response when the damping coefficient changes from 6.2 to 2.54 N·s/m.The system responses show stronger nonlinear phenomenon with the decreasing damping coefficient.Figure13b displays the voltage response when the mass changes from 174.1 to 77.4 g.The system responses show stronger nonlinear phenomenon with the increasing mass.

Fig.14 The comparisons for the output power under different excitation levels

Table2 Comparison of b and width for different systems

5.3 Performance of the energy harvesting

Figure14 shows the output power in different excitation levels from 0.3g to 2.5g.The solid and the dash lines are the stable and the unstable solutions,respectively.The black dash line shows the linear natural frequency.As shown in Fig.14,the harvester has the liner characteristic in the 0.3g excitation,the soft characteristic in 1g and 1.5g excitations and both the s oft and the hard characteristic in2g and 2.5g excitations.The b and width(frequency region between half power points)is exp and ed obviously with the increasing excitation level and the data is listed in Table 2.The b and width is21.3 Hz of the harvester here,15.8 Hz of the typical hardening nonlinearity and 10.5 Hz of the equivalent linear system in 2.5g excitation,which may expand the bandwidth 50.48%to the typical hardening nonlinearities and 102.86%to the equivalent linear system,respectively.The typical hardening nonlinear system has the same parameters(m,c,k1,η and Cp)with the proposed system, and the hardening nonlinear stiffness kh3is calculated by kh3??1.13×1011(Ep?74 GPa is the Young’s modulus of piezoelectric ceramics)base on the von Karman plate theory.The results demonstrate that the coexistence of hardening and softening stiffness can broaden the b and width of the harvester.However,a rather high excitation is needed to induce the nonlinearity of the harvester.

Figure15 shows the output surface of the voltage and the power as a function of frequency for the changing resistance values.For each of these curves,the point at which the voltage or the power is at a maximum for a given resistance is identified.The gray curves are the locus of these points and their projection curves are shown in the resistance–voltage,the frequency-voltage and the resistance-frequency planes in Fig.15a and in the resistance-power,the frequency-power and the resistance-frequency planes in Fig.15b,respectively.In Fig.15a,the projection curves of the resistance–voltage show that the maximum voltage increases with the resistance and tends to a constant for the resistance beyond 100 k?.The projection curves of the frequency–voltage show the backbone line of the voltage-frequency response and the line bends to the right with the increasing resistance.In Fig.15b,the project curve of the resistance-power shows that the maximum power is not monotone with the change of the resistance and it has a peak value when the resistance approaches 30 k?.The project curve of the frequency-power shows the backbone line of the power–frequency response and it is different from the backbone of the voltage-frequency response.However,the frequency of the peak point also increases with the increasing resistance.

Fig.15 a Voltage surface and b power surface over a range of the resistance values at 2g excitation

6 Conclusions

This investigation focuses on nonlinear vibration of a circular composite plate in a piezoelectric energy harvester.The harvester consists of a circular composite piezoelectric plate and a pro of-mass at the plate center.The method of harmonic balance enhanced by the arc-length continuation technique is proposed to seek the response of the harvester under different excitations.The method of incremental harmonic balance is developed to compare with the method of the harmonic balance.The nonlinear characteristics of the amplitude-frequency response are analytically,numerically,an experimentally investigated.The investigation yields the following conclusions.

1.The method of harmonic balance can solve the strongly nonlinear harvester governing equations and pass through turning points with the assistance of the arc-length continuation technique.

2.The stable approximate analytical solutions are confirmed by the numerical integrations and the experimental data.

3.The harvester has softening characteristic under small excitation and both hardening and softening characteristics under sufficiently large excitations.

4.The co-existence of the softening and hardening the characteristics broadens the working b and width of the harvester.

5.The backbone of the voltage-frequency response bends to the right with the increasing resistance and the peak frequency in the power–frequency response increases with the increasing resistance.

AcknowledgementsThis work was supported by the National Natural Science Foundation of China(Grants51575334 and 11802170),the State Key Program of National Natural Science Foundation of China(Grant 11232009),the Key Research Projects of Shanghai Science and Technology Commission(Grant 18010500100), and the Innovation Program of Shanghai Municipal Education Commission(Grant 2017-01-07-00-09-E00019).