Approximate solutions and their stability of a broadband piezoelectric energy harvester with a tunable potential function

Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Research paper

Approximate solutions and their stability of a broadband piezoelectric energy harvester with a tunable potential function Feng Qian a, Shengxi Zhou a,b, Lei Zuo a,∗ a b

Department of Mechanical Engineering, Virginia Tech, Blacksburg, VA 24061, USA School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China

a r t i c l e

i n f o

Article history: Received 10 October 2018 Revised 14 June 2019 Accepted 28 August 2019 Available online 30 August 2019 Keywords: Energy harvesting Nonlinear dynamics Monostable and bistable Harmonic balance Stability analysis

a b s t r a c t A broadband piezoelectric energy harvester (PEH) with a mechanically tunable potential function is modeled and analytically analyzed. The harvester consisting of a beam and a pre-compression spring at one end can be tuned to both monostable and bistable configurations. The axial motion of the beam resulting from the transverse vibration and spring load induces two coupled higher-order terms of displacement, velocity and acceleration into the governing equations. This significantly complicates the theoretical analysis, especially the stability analysis of solutions. Harmonic balance method is employed to investigate the dynamic characteristics of the nonlinear energy harvester. An effective approach is developed to solve the entries of the Jacobian matrix for determining the stability of analytical solutions. This approach offers a criterion for solution stability analysis of congeneric nonlinear systems with coupled higher-order terms. The energy harvesting performance and the nonlinear dynamic characteristics of the proposed PEH are explored for various base excitation levels, electrical resistive loads and pre-deformations of the spring. The approximate analytical solutions are validated by numerical simulations. Results demonstrate that the energy harvesting performance of the proposed PEH could be effectively tuned by the pre-deformation of the spring. The proposed PEH could harvest vibration energy in a wide frequency range of 0–91 Hz at the excitation level of 0.5 g. Published by Elsevier B.V.

1. Introduction Energy harvesting from ambient vibrations has received increasing attention in the last decade motivated by low-power electronic devices and wireless sensors [1–4]. The most widely studied energy harvesting device is the piezoelectric energy harvester (PEH) consisting of a beam, piezoelectric transducers and a tip mass [5,6]. PEHs are easily shaped, small in size and have high power densities. However, one of their inherent disadvantages is the narrow frequency bandwidth [7], which renders them only efficient at the resonant frequencies. Techniques including the introduction of nonlinear magnetic force, geometric nonlinearity and preload have been developed to address this issue. These methods effectively expend the frequency bandwidth of the harvesters by changing the system stiffness [8]. Nonlinear energy harvesters can be mainly categorized into monostable [9], bistable [10–14] and tristable [15–18] systems. These nonlinear energy harvesters have been theoretically and experimentally demonstrated to ∗

Corresponding author. E-mail addresses: [email protected] (F. Qian), [email protected] (S. Zhou), [email protected] (L. Zuo).

https://doi.org/10.1016/j.cnsns.2019.104984 1007-5704/Published by Elsevier B.V.

2

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

own a much wider bandwidth [19,20]. For instance, a broadband internally-resonant vibratory energy harvester was realized by using a magnetically tunable L-shape piezoelectric beam structure [21]. The knowledge achieved in nonlinear energy harvesting are also suitable for vibration control, shock absorption and mechanical signal processing [22]. Besides magnet based nonlinearity, geometrically nonlinearity have also been widely investigated in energy harvesting [9,23–27]. Chen et al. [28] designed a two DOF nonlinear energy harvester and explored the snap-through mechanism at internal resonance. A nonlinear piezoelectric micro-plate harvester with axial pre-load was analytically studied for harvesting energy from pulsating arterial pressure [29]. The effect of the axial load on the pre-buckling and post-buckling nonlinear dynamics and voltage outputs were studied. Buckling beam harvesters are receiving more attention for their simple structures and high-efficiency power outputs [30–35]. Leland and Wright [30] designed a resonance-changeable PEH with the adjustable axial preload achieved by using a differential micrometer drive. Xu et al. [31] systematically explored the influence of the axial compressive displacement on the energy harvesting performance, and found that the open-circuit RMS voltage could be improved by 300%. Ando et al. [32] studied a bistable beam with a proof mass at the middle, and experimentally tested the influence of external electrical load on the peak power. Liu et al. [33,34] developed a bistable piezoelectric stack energy harvester based on the buckled spring-mass architecture. The buckling bistable PEH can continuously snaps from one stable equilibrium to the other under proper initial conditions and excitations. This snap-through behavior (also referred to as interwell oscillation) exhibits versatile features like large-amplitude vibration, which is favorable to energy harvesting compared with the small-amplitude intrawell oscillation. Therefore, the investigation of the snap-through mechanism of the bistable PEHs is critical to their application in energy harvesting. However, modeling and analysis of buckling beam harvesters is usually a daunting task, especially the stability analysis of solutions. This is largely due to the appearance of some coupled higher-order nonlinear terms of displacement, velocity and acceleration in the governing equations. These coupled higher-order terms have appeared in many nonlinear systems [9,35–47]. Among these, numerical integration methods are widely used to simulate the nonlinear equations in the time domain [35–38]. Numerical simulation is inconvenient and time consuming for parameter analysis. Approximate analytical solutions of the congeneric systems are also obtained in the frequency domain by the method of multiple-scale analysis [39–42], which has been proven to be an effective way to system parameter analysis. The method of multiple-scale analysis has good accuracy for small excitation levels. He’s parameter expanding method (HPEM) [43,44] is used to obtain the first-order analytical response of a similar nonlinear system in the time domain. Besides, the residue harmonic balance and energy balance methods [45–47] are also successfully applied to approximate the periodic response of nonlinear systems with coupled higher-order terms. Although informative knowledge for this kind of a specific nonlinear system has been provided, the theoretical study of the snap-through behavior has not been adequately explored yet [22], especially the stability analysis of solutions. This paper proposes a novel broadband nonlinear PEH, studies the solution stability, dynamics and the energy harvesting performance. Specifically, a constrained spring is exerted to one moveable end of the PEH for the realization of the axial compression. Approximate analytical frequency responses are obtained via the harmonic balance method. The coupled higher-order nonlinear terms in the governing equations significantly complicate the stability analysis of the solutions. A new approach is developed to analytically derive the entries of the Jacobian matrix for the stability analysis of solutions. The approach can be generalized and used for similar nonlinear systems with coupled higher-order terms. The approximate analytical solutions are in good agreements with the numerical simulations. The performance of the monostable and bistable PEHs under different external resistive loads, excitation levels and the tuning pre-deformation of the spring are investigated. Results demonstrate that the proposed tunable nonlinear energy harvester can scavenge energy in a wide frequency range of 0–91 Hz at the excitation level of 0.5 g. 2. Theoretical model The proposed novel nonlinear PEH is constituted of a steel beam, two piezoelectric patches symmetrically bounded on the top and bottom of the beam, and a spring preloaded at one end, as illustrated in Fig. 1(a). Two guide rails constrain the vertical motion and free the axial motion of the clamped-clamped beam to achieve the tunable function and buckling mechanism, as shown in the Fig. 1(b). The PEH is subjected to a transverse base excitation y¨ (t ). To prevent charge cancellation, only part of the piezoelectric patches is covered by electrode layers. The substrate, piezoelectric and electrode layers are assumed to be perfectly bonded together so that no relative motions happen between them. The PEH can be modeled as a current source and a capacitor [3,6], which are connected with an external resistance, as shown in Fig. 1(d), where v(t) is the voltage across the external resistive load R. Let ρ s and ρ p denote the densities of the substrate and piezoelectric material, respectively. The length and width of the piezoelectric layers are identical to those of the substrate and denoted as L and b, respectively. The mass of per unit length of the composite beam is m = ρs bhs + 2ρ p bh p , where hs and hp are the height of the substrate and piezoelectric layers. The length of the electrode is le . The transverse and axial displacements of the middle plane of the harvester are w(x, t) and u(x, t), respectively. Ignoring the effect of the very thin electrode layers on the mass and stiffness of the system, the kinetic energy of the system can be written as:

T =

1 2



L 0



2



m[w˙ (x, t ) + y˙ (x, t )] + u˙ 2 (x, t ) dx

(1)

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

3

Fig. 1. (a) Schematic of the proposed nonlinear PEH, red and blue dashed lines indicate the stable and unstable positions; (b) close-up view of the movable end; (c) cross-section of the piezoelectric harvester; (d) equivalent circuit [3,6].

·

where the dot (·) denotes the derivative with respect to time, and y (x, t ) is the velocity of the base. The total potential energy of the system includes mechanical bending, axial compression energy, spring potential energy and electromechanical coupling energy, which can be written as:

U=

1 2

 T11 S11 dV −

Vs

1 2

 Vp

D3 E3 dV −

1 k u(L, t )2 − kd u(L, t ) 2 d

(2)

where S11 and T11 are the strain and stress in the axial direction, respectively; D3 and E3 are the electric displacement and intensity along the polarization direction; kd and  are the stiffness and pre-deformation of the spring. Vs and Vp are the volume of the substrate and piezoelectric layers; Assume  is very small, such that the preload kd  is almost constant during the small movement of the beam. Therefore, the potential energy induced by the spring preload can be written as the last term in Eq. (2). In terms of the inextensibility [9,48], the axial displacement u (x, t) can be expressed by x u(x, t ) = − 21 0 [w (x, t )]2 dx. Assuming the damping coefficient of the system is c and the electric charge generated by the harvester is q(t). The work done by the generalized non-conservative forces can be written as:

 W =−

L 0

cw˙ (x, t )w(x, t )dx − q(t )Rq˙ (t )

(3)

where q˙ (t ) = − v(Rt ) is the current flowing in the circuit. p s =E S The stress in the substrate and piezoelectric layers are T11 s 11 and T11 = E p S11 − e31 E3 , where Es and Ep are the Young’s modulus of the substrate and piezoelectric material. The electric displacement in the piezoelectric layers can be expressed as s E where e s D3 = e31 S11 + 33 3 31 and 33 are the piezoelectric coupling constant and relative permittivity measured at constant strain. E3 is the electric intensity approximated by E3 = − vh(t ) . Substitution of the stresses, strains and displacement into p

Eq. (2), the total potential energy can be written as:

U =

1 2



L

0



2

E I w (x, t ) dx +

− kd  E bh3

1 2



L 0



EA 2



L

1 2

0

 2 2

w (x, t )

dx −

1 k 2 d

L  1 2

0

2

2

w (x, t ) dx

 le  s 2 b 2 1 L 33 w (x, t ) dx − 2 e31 bh p (hs + h p )w (x, t )v(t )dx − v (t )dx 0

2

0

hp

(4)

E

where EI = s12 s + 6p bh p (4h2p + 6h p hs + 3h2s ) and EA = Es bhs + 2E p bh p are the equivalent bending and axial stiffness of the composite beam, respectively. Generally, the vibration system is dominated by the fundamental mode [32]. Research has showed that the first modal solution had a fairly good agreement with the simulation of partial differential equation including multiple modes [49]. Therefore, only the first mode is considered and the mode shape function is [31–33]:

ϕ (x ) =





1 2π x 1 − cos 2 L

(5)

4

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

Using Galerkin method and the separation of variables, the displacement can be approximated by w(x, t ) = η (t )ϕ (x ), where η(t) is the modal displacement. Substitution of the approximated displacement into Eqs. (1), (3) and (4) yields:

m T = 2



L



0

 U =

η2

EI 2

m (x ) + 2η˙ y˙ ϕ (x ) + y˙ dx + (ηη˙ )2 2

where C p =

s bL 233 hp



L



x



2

2 ϕ  (x ) d x

η˙ ϕ



⎡   2 ⎤  L  L  L 2 4 2     k  k EA 2 ϕ  (x ) dx − d ϕ  (x ) d x + η 4 ⎣ ϕ  (x ) d x − d ϕ  (x ) d x ⎦

L

0



2

2

− ηbe31 v(t )(hs + h p ) W = −cηη˙



2

2



L 0



Le 0

0

0

8

0

0

dx

(6)

8

0

1 2

ϕ  (x )dx − Cp v2 (t )

(7)

ϕ 2 (x )dx + q(t )v(t )

(8)

is the total capacitance of the piezoelectric layers. Using Lagrange method, the equations of motion of the

system are obtained by

  η¨ + 2ξ ωn η˙ + β ηη˙ 2 + η2 η¨ + (k1 − ks )η + (k2 − k3 )η3 − θ v = −μy¨

(9a)

v˙ + τ v + θ¯ η˙ = 0

(9b)

where ξ is the damping ratio; ωn =



k1 ; k1 , k2 and k3 are the scaled linear and nonlinear stiffness induced by the bending, axial stiffness of the composite structure and spring stiffness, respectively. ks is the effective spring stiffness. θ and θ¯ are the scaled electromechanical coupling stiffness of the piezoelectric layers. The system parameters in Eqs. (9a) and (9b) are given in Appendix A. Eq. (9a) shows that the introduction of the nonlinearity results in two coupled higher-order terms, i.e. the product of the displacement and the square of velocity, as well as the product of the square of displacement and acceleration. It will be shown later that these coupled higher-order nonlinear terms bring great trouble to the stability analysis of solutions. The free vibration system will always come to one of the equilibria due to damping effect, where the velocity and acceleration of the system are zeroes. Therefore, by setting all the time derivatives in Eq. (9a) to zeros, the equilibria can be obtained as:

 η1 = 0 , η2 , 3 = ± −

k1 − ks k2 − k3

(10)

In general, the scaled axial stiffness k2 is much larger than the spring stiffness k3 , i.e.

k3 k2

=

kd 3EA/L

< 1. For such a case,

the stability of the system is essentially determined by the effective spring stiffness ks induced by the pre-deformation . For instance, when  is very small (ks < k1 ), only the trivial solution η1 = 0 exists and the harvester is monostable. As  increases to a certain value that makes ks = k1 , the axial force in the PEH induced by the compressive spring equals to the critical buckling load of the beam. The PEH begins to enter into the post buckling regime after that. According to Eq. (9a), the potential energy function of the system can be written as:

( k1 − ks ) 2 ( k2 − k3 ) 4 η + η V¯ (η ) = 2 4

(11)

which shows that the effective spring stiffness ks determines the unbuckling and buckling states of the system. When ks < k1 , the harvester is unbuckling and only oscillates around the absolute minimum of the potential function at the equilibrium η1 = 0. For the case  of ks > k1 , it is in the buckling state and the quartic energy function has two minimums at k −k

the two local equilibria η2,3 = ± − k 1−k s , therefore, the system is bistable. 2

3

Fig. 2 depicts the potential energy function for the tunable pre-deformation of the spring at a constant spring stiffness of kd = 30 0 0N/m. The potential energy function gradually turns from a flatted parabolic shape to a double-well shape as  increases. This intuitively illustrates that the PEH varies from the monostable state that only has one potential well at the single equilibrium to the bistable state which has two potential wells at the two equilibria separated by a saddle point. There is a critical state between the mostable and bistable regimes, at which the pre-deformation of the spring is denoted by c . To achieve the large-amplitude interwell oscillation, the bistable PEH needs to snap through one potential well to the other by overcoming the barrier (saddle point) separating the wells. Therefore, a larger excitation is required for activating the snap-through as the pre-deformation of the spring increases. This is because the potential wells become deeper and the barrier becomes higher as the pre-deformation increases. Alternatively, tuning the potential wells in terms of the excitations can also facilitate the occurrence of interwell oscillations. For instance, reducing the pre-deformation of the spring could shallow the potential wells and low the barrier so that the large voltage output is still achievable when the excitation is

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

5

Fig. 2. Potential energy of the PEH with different pre-deformations of the spring.

small. When the excitation level is higher, properly increasing the pre-deformation of the spring can enlarge the vibration amplitude of interwell oscillations and thus produce larger voltage. The tuning process could be realized by applying an actuator in the axial direction of the harvester. Based on the analysis, the PEH should be tuned to monostable configuration when subjected to a small excitation and to bistable configuration with deep potential in a large excitation environment. 3. Harmonic balance (HB) solutions The harmonic balance (HB) method [35,50–52] is adopted to solve Eqs. (9a) and (9b) for approximate analytical solutions. The PEH is subjected to a harmonic base excitation of y¨ (t ) = Y cos(ωt ), where Y and ω are the amplitude (level) and frequency. Considering the fundamental harmonic, the displacement and voltage responses are assumed as:

η (t ) = a0 (t ) + a1 (t )sin(ωt ) + b1 (t )cos(ωt )

(12a)

v(t ) = c1 (t )sin(ωt ) + d1 (t )cos(ωt )

(12b)

where a0 , a1 , b1 , c1 , and d1 are unknown coefficients to be determined, which vary slowly in time. The coefficient a0 is the offset of the displacement from the saddle point (zero position) when the buckling PEH is in intrawell oscillations. For interwell oscillations, a0 is zero because of the symmetrical motion. Substituting Eqs. (12a) and (12b) into Eq. (9b), and neglecting the time derivatives, one has

c1 =

ωθ¯ ( τ b1 − ω a1 ) ω2 + τ 2

(13a)

d1 =

−ωθ¯ ( τ a1 + ω b1 ) ω2 + τ 2

(13b)

Substituting Eqs. (12a), (12b) and (13a), (13b) into Eq. 9(a), ignoring the higher harmonic terms, and balancing the terms of constant, sin(ωt) and cos(ωt), three differential Eqs. (B1)–(B3) of the unknown coefficients a0 , a1 , b1 and their derivatives are obtained and given in Appendix B. Considering the steady state solutions only by neglecting the first and second time derivatives in Eqs. (B1)–(B3), the following Eqs. (14a)–(14c) can be derived



k¯ 2 a30 +





3k¯ 2 − βω2 2 k¯ 1 + r a0 = 0 2



r2 b1 k¯ 1 − ω2 + ω − βω2 a20 + 2

 a1



r2 k¯ 1 − ω2 + ω − βω2 a20 + 2

(14a)



+ k¯ 2 3a20 +





3r 2 4

3r 2 + k¯ 2 3a20 + 4

θ¯ θ and r = where k¯ 1 = k1 − ks , k¯ 2 = k2 − k3 , = ωω2 + τ2 

a0 has three solutions (a0 = 0 and a0 = ± ¯1

k2



βω2 −3k¯ 2

cillations around the two local equilibria. Only if

2

 + a1 ( τ + 2ξ ωn ω ) = −μY

(14b)

− b1 ( τ + 2ξ ωn ω ) = 0

(14c)



a21 + b21 is the amplitude of the displacement. Eq. (14a) shows that

r 2 − k¯ 1 ) which are corresponding to the interwell and intrawell os-

βω2 −3k¯ 2 2 r − k¯ 1 > 0, i.e. 2

ω>



2k¯ 1

β r2

3k¯

+ β2 , the solutions of the intrawell

6

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

oscillation exists. Squaring the both sides of Eqs. (14b) and (14c), and adding the resultant equations together yield:





r2 r 2 k¯ 1 − ω2 + ω − βω2 a20 + 2



+ k¯ 2 3a20 +

3r 2 4

2 + r 2 ( τ + 2 ξ ω n ω ) = μ2 Y 2 2

(15)

The displacement amplitude of the system can be solved from Eqs. (14a) and (15). The voltage amplitude is obtained by:

V =



ωθ¯ r ω2 + τ 2

c12 + d12 = √

(16)

The stability of the solutions is determined by the real part of the eigenvalues of the Jacobian matrix [50–52]:

⎡ ∂ a˙

0

∂ a0 ⎢ J (a0 , a1 , b1 ) = ⎣ ∂∂ aa˙ 0 1 ∂ a˙ 0 ∂ b1

∂ a˙ 1 ∂ a0 ∂ a˙ 1 ∂ a1 ∂ a˙ 1 ∂ b1

⎤

∂ b˙ 1 ∂ a0 ∂ b˙ 1 ⎥ ∂ a1 ⎦ ∂ b˙ 1 ∂ b1

(17)

For a general nonlinear system without the coupled higher-order terms, the entries of the Jacobian matrix can be directly derived as the explicit functions of the coefficients [51]. However, that is impossible for the presented system, since the resulting equations by ignoring the second derivatives of a0 , a1 and b1 from Eqs. (B1) to (B3) are second-order implicit functions of a˙ 0 , a˙ 1 and b˙ 1 , as shown in Eqs. (B4)–(B6) in Appendix B. Directly solving for a˙ 0 , a˙ 1 and b˙ 1 as usually done for general nonlinear systems is impossible. This is ascribed to the presence of the quadratic terms in Eqs. (B4)–(B6), which are induced by the coupled higher-order terms of the acceleration and the square displacement (η2 η¨ ), and the displacement and the square velocity (ηη˙ 2 ) in Eqs. (9a) and (9b). An alternative way to get the entries of the Jacobian matrix is proposed and elaborated in Appendix B. Taking the partial derivatives of Eqs. (B1)–(B3) with respect to each of the coefficients gives nine individual equations about the nine entries of the Jacobian matrix. By combining the resultant nine equations, all the entries of the Jacobian matrix can be solved as the functions of the coefficients. This approach can be applied to the congeneric systems with coupled higher-order nonlinear terms, e.g., see [9,36–47], for determining the stability of the harmonic balance solutions. The moveable end of the clamped-clamped beam is realized by two guide rails in this study. Many other ways, including different boundary conditions [25,26,53], have been reported that could introduce the similar coupled higher-order nonlinear terms in nonlinear systems. The proposed method for the stability analysis of the HB solutions can be extended to these systems. Substituting the solutions of the coefficients a0 , a1 and b1 obtained from Eqs. (14a), (14b) and (14c) into Eq. (17), one can further attain the corresponding eigenvalues of the Jacobian matrix. If the real part of all the eigenvalues are negative, the corresponding solution is stable. Otherwise, it is unstable. 4. Approximate analytical solutions and numerical verification Three different configurations (one monostable, one bistable with shallow potential wells, and one bistable with deep potential wells) are considered by tuning the pre-deformation of the spring to ࢞ = 3.0%L, ࢞ = 4.5%L, ࢞ = 6.0%L. Three base excitation levels of 0.015 g, 0.15 g and 0.5 g are considered as low, medium and high excitations. Excitation amplitude and resistance sweeps are performed at the three excitation frequencies of 5 Hz, 30 Hz, and 60 Hz to study their influence on the dynamic characteristics and voltage output. The dynamic characteristics of the PEH at different pre-deformations of the spring are studied. Simulations are conducted by using numerical integration method to validate the approximate analytical solutions. Both forward and backward frequency sweeps at a very slow rate are performed to each of the three configurations. Simulations are also carried out under the harmonic excitations at the discrete frequency points, from which the steady-state amplitudes of the voltage and displacement responses are collected. Unless explicitly stated otherwise, the initial displacement, velocity and voltage are assumed to be zero in simulation. The geometric and material properties of the PHE are summarized in Table 1. 4.1. Frequency responses A large external resistance of R = 250 MΩ is used to approximate the open circuit condition. The voltage and displacement frequency curves of the monostable PEH obtained from the HB method are plotted in Fig. 3 (a-1) and (a-2) under the low excitation level of 0.015 g, along with the frequency sweep simulation results. The hysteretic behavior can be clearly observed from the difference between the forward sweep (from low to high) and backward sweep (from high to low). The circles represent the steady-state amplitudes simulated under harmonic excitations at discrete frequency points. The analytical solutions show a good agreement with the frequency sweep simulations. The steady-state amplitude agrees well with the backward frequency sweep results attributed to the zero initial condition. Three branches, namely, the resonant (upper or higher) branch, non-resonant (lower) branch and unstable branch can be observed from the frequency curves.

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

7

Table 1 The geometric and material properties of the PEH. Substrate (steel) and spring

Piezoelectric layer

Parameter

Value

Parameter

Value

Length L Width b Thickness hs Density ρ s Young’s modulus Es Damping ratio ξ Spring stiffness kd

200 mm 30 mm 0.5 mm 7850 kg·m−3 203 Gpa 0.011 3000 N·m−1

Length Lp Width b Thickness hp Density ρ p Young’s modulus Ep Piezoelectric constant d31 Relative permittivity  33 Permittivity of free space  0

200 mm 30 mm 0.08 mm 4000 kg·m−3 40 Gpa −10 pC·N−1 100 8.854 × 10−12 F·m−1

Fig. 3. Voltage and displacement responses of the monostable PEH at different excitation levels; (a-1), (b-1) and (c-1): voltage responses at Y = 0.015 g, Y = 0.15 g and Y = 0.5 g; (a-2), (b-2) and (c-2): displacement responses at Y = 0.015 g, Y = 0.15 g and Y = 0.5 g.

8

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

The frequency curves bend to the right hand due to the nonlinearity which consequently broadens the bandwidth. The results of the system subjected to medium and high excitation levels are presented in Fig. 3 (b-1) and (b-2), (c-1) and (c-2), respectively. The monostable PEH experiences a unique higher branch oscillation at the low frequency range of 0–17.6 Hz, as marked in Fig. 3 (a-2). As the excitation level increases, this frequency range is further enlarged to 0–23.5 Hz at Y = 0.15 g in Fig. 3 (b-1) and (b-2), and to 0–30.5 Hz at Y = 0.5 g in Fig. 3 (c-1) and (c-2). From the perspective of energy harvesting, the system is expected to oscillate at the higher branch. However, as the excitation frequency increases, the system enters into the multi-solution range (two stable solutions and one unstable solution), as marked in Fig. 3 (a-2). Fig. 3(a)–(c) show that the multi-solution range is enlarged from 17.6–20.7 Hz to 23.5–52 Hz, and 30.5–92.5 Hz as the excitation level increases from 0.015 g to 0.15 g and 0.5 g. This frequency range is favorably exploited by researchers for the large voltage output by activating the higher branch vibration via appropriate initial conditions and excitation levels [19]. As the excitation frequency further increases, the system suddenly jumps to the undesirable lower branch oscillation. The bistable PEH exhibits more rich dynamic characteristics including the global snap through dynamics and local oscillations in one of the potential wells. The approximate analytical solutions and numerical frequency sweep results for the bistable PEH with shallow wells are presented in Fig. 4(a)–(c) for the three excitation levels. The analytical results in Fig. 4 (a-1) and (a-2) show that the desired large-amplitude vibration of the stable interwell oscillation only takes place in a very narrow frequency range of 0–11 Hz at Y = 0.015 g and always coexisted with the stable intrawell oscillations. The smallamplitudes stable intrawell oscillation occurs in the entire frequency range. The voltage output is very small because the system is confined to the local potential wells. The steady-state amplitudes simulated at discrete frequency points also show the local vibrations and agree with the stable intrawell solutions. The local intrawell oscillation could be clearly observed from the displacement response around the equilibrium (−3.7 mm) in Fig. 4 (a-2), where the dashed black lines are the two local equilibria. The frequency bandwidth of the large-amplitude interwell oscillation is extended to 0–49.5 Hz at the excitation level of 0.15 g, as shown in Fig. 4 (b-1) and (b-2), and the voltage output is also enhanced by the snap-through dynamics. The forward frequency sweep follows the large-amplitude vibration and the amplitudes agree with the interwell frequency responses. The steady-state amplitudes follow the backward frequency sweep simulation. Fig. 4 (b-2) demonstrates that the system maintains the large-amplitude vibration in the forward frequency sweep by alternatively snapping through the two local equilibria in a wide range of the frequency. This suggests that the large energy input could maintain the desirable interwell oscillation. As the excitation level is further increased to 0.5 g, the large voltage output generated by the large-amplitude interwell oscillation occurs in a much larger frequency range of 0–91 Hz in the forward frequency sweep, as shown in Fig. 4 (c1). Fig. 4 (c-2) shows that the snap through phenomenon is successfully achieved over a wider range of frequency. It is noteworthy that a special frequency gap of 0–20 Hz appears, as marked in Fig. 4 (c-1), in which only the stable interwell oscillation exists. The steady-state amplitudes simulated at discrete frequency points follow the backward frequency sweep results. The forward frequency sweep simulation is in good accordance with the HB solutions. The results of the bistable PEH with deep potential wells are plotted in Fig. 5(a)–(c) for the three excitation levels, respectively. The analytical solutions in Fig. 5 (a-1) and (a-2) suggest that the frequency range of the stable interwell oscillation is 0–7 Hz, which is narrowed down due to the deeper potential wells compared with this in Fig. 4 (a-1) and (a-2). Nevertheless, the voltage amplitudes becomes larger. No snap through occurred in either forward or backward frequency sweeps because of the low excitation. The steady-state amplitudes also only exhibit local vibrations. Fig. 5 (b-1) and (b-2) indicate that the system still experiences local vibrations over most of the frequency range under the medium excitation level. This is quite different from the results presented in Fig. 4 (b-1) and (b-2), which exhibit the large amplitude interwell oscillations under the same excitation level. This suggests that the bistable PEH with shallow potential wells outperforms the one with deep potential wells at the medium excitation. The steady-state amplitude responses simulated at discrete frequency points follow frequency sweep results but display the large-amplitude snap through dynamics at around 20 Hz and super-harmonic vibrations in the frequency range of 60–80 Hz. Fig. 5 (c-1) and (c-2) show that the large-amplitude interwell oscillations are achieved in a wide frequency range of 0–90 Hz, when the bistable PEH with the deep potential wells is subjected to the excitation level of 0.5 g. The forward frequency sweep agrees well with the stable interwell solutions. The steady-state amplitudes follow the backward frequency sweep. The super-harmonic vibration is observed at 69 Hz from both the backward frequency sweep and the steady-state amplitudes, as shown in Fig. 5 (c-2). The super-harmonic vibration is not captured in the analytical solutions since only the steady-state solutions of the fundamental harmonic are considered in the HB analysis. The voltage frequency responses of the mono and bistable PEHs under the three excitation levels are plotted together in Fig. 6 for comparison. The ‘Bistable-S’ and ‘Bistable-D’ in the legend stand for the bistable systems with shallow and deep potential wells, respectively. Only the higher branch solutions of the monostable PEH and the stable interwell solutions of the bistable PEHs are presented because the associated large voltage output is concerned in energy harvesting. It is shown that the frequency range of the large voltage outputs increases along with the excitation levels. At the same excitation level, the bistable PEH with the deep potential wells has the largest voltage output, while the monostable system has the smallest voltage output. It should be noted that only when the interwell oscillation occurs the bistable PEH with the deep potential wells could outperform the one with the low potential wells and monostable PEH. However, the results demonstrate that

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

9

Fig. 4. Voltage and displacement responses of the bistable PEH with shallow wells at different excitation levels; (a-1), (b-1) and (c-1): voltage responses at Y = 0.015 g, Y = 0.15 g and Y = 0.5 g; (a-2), (b-2) and (c-2): displacement responses at Y = 0.015 g, Y = 0.15 g and Y = 0.5 g.

the bistable PEH with deep potential wells is more difficult to achieve the large-amplitude vibration at low and medium excitation levels than the one with low potential wells and monostable system. In summary, the system should be tuned to the monostable configuration when the excitation is too weak to activate the interwell oscillation. Once the excitation level is high enough to realize the snap-through oscillation, the system should be tuned to the bistable configuration for a larger voltage output. The bistable PEH with shallow potential wells is more preferable for the energy harvesting at medium excitation level. Since the shallow potential wells separated by the flatter barrier is much easier to overcome under low-level excitation. 4.2. Influence of the excitation levels To investigate the influence of the excitation level on the energy harvesting performance, the excitation amplitude sweep is carried out at three frequencies of 5 Hz, 30 Hz and 60 Hz, which represent lower, medium and higher frequencies. Only stable solutions are presented for brevity. The voltage outputs of the monostable and bistable PEHs at the three excitation frequencies are given in Fig. 7(a)–(c), respectively. The voltage outputs increase along with the excitation level at the fixed

10

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

Fig. 5. Voltage and displacement responses of the bistable PEH with deep wells at different excitation levels; (a-1), (b-1) and (c-1): voltage responses at Y = 0.015 g, Y = 0.15 g and Y = 0.5 g; (a-2), (b-2) and (c-2): displacement responses at Y = 0.015 g, Y = 0.15 g and Y = 0.5 g.

Fig. 6. Voltage frequency responses of the monstable and bistable PEHs under different excitation levels.

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

11

Fig. 7. Analytical voltage responses of the monostable and bistable PEHs at varying excitation levels: (a) 5 Hz; (b) 30 Hz; (c) 65 Hz.

frequency. The interwell oscillation of the bistable PEH with deep potential wells produces largest voltage. Both the bistable PEHs with shallow and deep wells outperform the monostable one when the interwell oscillation take places. Fig. 7(a) shows that only the higher branch solution exists over all the considered excitation levels for the monostable PEH at 5 Hz. The interwell solution coexists with two intrawell solutions over the most of the excitation levels for the bistable PEH with deep potential wells. The bistable PEH with shallow wells has one interwell and two intrawell solutions at low excitation levels, specifically in the range of 0–0.11 g, and only has the interwell solution when the excitation level exceeds 0.11 g. Fig. 7(b) shows that both the higher and lower branch solutions appear for the monostable PEH in the low excitation levels of 0–0.45 g. One intrawell solution appears along with the interwell solutions for the bistable PEH with shallow wells. But two inrawell solutions at range of 0–0.25 g appear for the bistable PEH with deep wells. As the excitation level increases, only one intrawell solution exists, which disappears after 0.85 g. As the frequency is increased to 60 Hz, the voltage outputs of all the systems increase along with the excitation level, as shown in Fig. 7(c). Only the lower branch solution of the monostable PEH and the intrawell solution of the bistable PEHs exist when the excitation level is less than 0.21 g. After that, the interwell and intrawell solutions always coexist.

4.3. Influence of the external electrical resistive loads The resistance sweep is performed to the three nonlinear PEHs subjected to the excitation level of Y = 0.5 g at the three frequencies of 5 Hz, 30 Hz and 60 Hz. Only higher branch solution of the monostable PEH and interwell solutions of the bistable PEHs are considered. The analytical power outputs of the harvesters over varying resistance are plotted in Fig. 8(a)– (c), respectively. The power of all the PEHs firstly increases along with the growth of the resistance and then comes to the maximum, after which it begins to decay. For the three different excitation frequencies, both the bistable PEHs produce larger power than the monostable one. While the bistable PEH with deep potential wells has a larger power output than the one with shallow potential wells. The three different PEHs have the same optimal resistance at each excitation frequency because they have the same capacitance. The values of the optimal resistance at the three excitation frequencies are identified as 0.48 MΩ, 0.08 MΩ and 0.04 MΩ, which are found to be in accordance with the analytically calculated values from the formula ω1C p [54].

12

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

Fig. 8. The power amplitudes of the monostable and bistable PEHs at varying resistive loads: (a) 5 Hz; (b) 30 Hz; (c) 60 Hz.

Fig. 9. The voltage response of the harvester with respect to varying pre-deformation of the spring at: (a) different excitation levels of 0.015 g, 0.15 g and 0.5 g; (b) different excitation frequencies of 5 Hz, 30 Hz and 60 Hz.

4.4. Influence of the pre-deformation of the spring The harvester can be tuned from monostable to bistable with varying potential wells by adjusting the pre-deformation of the spring. The above results demonstrate that the dynamic characteristics as well as the energy harvesting performance depends on the tuning pre-deformation of the spring. It is necessary to explore the influence of the tuning pre-deformation of the spring on the voltage output of the harvester subjected to different excitation amplitudes and frequencies. Fig. 9(a) presents the voltage amplitudes of the harvester versus the tuning pre-deformation of the spring at three excitation levels of 0.015 g, 0.15 g and 0.5 g and the frequency of 5 Hz. When the ratio of the pre-deformation of the spring over the length of the harvester is tuned to less than the critical value Lc × 100% = 3.7%, the system is monostable and only higher branch solutions are observed in Fig. 9(a). These higher branch vibrations are desirable under lower excitation level not enough to activate the snap through of the bistable harvester. When the ratio is tuned to greater than the critical value, the system enters in the bistable regime with two intrawell and one interwell solutions. The voltage amplitude of the

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

13

interwell oscillation has an increasing trend along with the increasing pre-deformation of the spring. One of the intrawell oscillations exhibits a decaying voltage output tendency while the other one has a rising tendency in the voltage output. This is because the distance at the two potential wells becomes larger but narrower as the pre-deformation increases. An evident region of the pre-deformation of the spring can be identified from Fig. 9(a), in which only the favorable interwell oscillation happens. Theoretically, it’s the best choice to tune the pre-deformation of the spring in this region, since the preferable large amplitude vibration will always happen regardless of the initial conditions. The region is enlarged along with the increase of the excitation level. Fig. 9(b) presents the voltage response over the tuning pre-deformation of the spring at the excitation level of 0.5 g and three different excitation frequencies of 5 Hz, 30 Hz, and 60 Hz. Similarly, the voltage amplitudes of the interwell oscillation increases as the pre-deformation of the spring becomes large. In addition, it is suggested that the higher excitation frequency produces larger voltage for the interwell motion. The monostable harvester ( Lc × 100% < 3.7%) has two stable solutions respectively belonging to the upper and lower branches at the high excitation frequencies of 30 Hz and 60 Hz, and a single upper branch solution at the low excitation frequency of 5 Hz. This observation follows the results in Fig. 3(a)–(c). With the augment of the spring pre-deformation ( Lc × 100% > 3.7%), the stable intrawell solutions at the excitation frequencies of 30 Hz, and 60 Hz are observed. An interval of the spring pre-deformation can be observed where only interwell oscillations occur. This interval gradually shifts to the large spring pre-deformation area as the excitation frequency increases. These findings offer a practical reference for the design of the nonlinear beam harvester subjected to different ambient excitation levels. The tuning of the potential wells could be achieved through active control strategies in real applications, which could be a potential research direction in the future. The results indicate that the pre-deformation of the spring should be tuned according to the excitation amplitude and frequency to achieve the maximum voltage output. Therefore, a close-loop feedback control system including power tracking and actuating elements is necessary for the proposed PEH with tuning potential function. The input of the control system may include the excitation amplitude, frequency, initial condition and voltage output of the system. Proper shape memory and piezoelectric actuators could be designed and integrated into the system as actuators to realize the tuning of the pre-deformation of the spring. 5. Conclusions A novel broadband piezoelectric energy harvester with mechanically tunable potential energy function is proposed for enhancing vibration energy harvesting. The theoretical model is derived based on the energy method and Euler–Lagrange equation and verified by numerical simulations. Harmonic balance method is used to derived the approximate analytical solutions of the PEH with the coupled higher-order nonlinear terms. An approach is proposed to solve for the entries of Jacobian matrix for determining the stability of the approximate solutions. This approach offers a reference for solution stability analysis of congeneric nonlinear systems with coupled higher-order terms. The energy harvesting performance, as well as the nonlinear dynamic characteristics of the harvesters (both monostable and bistable configurations) subjected to different excitation levels are investigated. The results demonstrate that the bistable PEH can produce large voltage output due to the large amplitude interwell oscillation. The harvester should be tuned to the monostable state for a large voltage output under the low-level excitation. The snap-through behavior of the bistable PEH depends on the depth of the tunable potential energy wells apart from the initial condition and the excitation condition. The harvester should be tuned to the bistable state with deep potential wells when the excitation could activate the interwell oscillation. The pre-deformation of the spring has a significant effect on the nonlinear dynamic characteristics and the energy harvesting performance. A region of the pre-deformation of the spring is identified that only the interwell oscillations happen. A close-loop control strategy may be integrated in the harvester to realize self-tuning function in terms of the excitation level, frequency, and initial conditions in future study. These findings have good value for the design of the nonlinear harvesters for enhancing energy harvesting performance. Acknowledgments This research was supported by National Science Foundation under grant no. 1529842. Compliance with ethical standards The authors declare that there is no conflict of interest concerning the publication of this manuscript. Apendix A: System parameters in Eq. (9)

 

ξ

c = , 2mω n

k3 =

β=

2 ∫L0 ∫x0

m ∫L0 ϕ 2 (x )dx

   2 2 ϕ  (x ) dx

kd ∫L0

8m ∫L0

ϕ 2 (x )dx

 2 2 ϕ  (x ) dx dx

,

θ=

, k1 =

   2 2 2  2 3EA ∫L0 ϕ  (x ) dx ϕ  (x ) dx kd  ∫L0 ϕ  (x ) dx , ks = , k2 = m ∫L0 ϕ 2 (x )dx m ∫L0 ϕ 2 (x )dx 8Lm ∫L0 ϕ 2 (x )dx

EI ∫L0



e31 b(hs + h p ) ∫l0e ϕ  (x )dx 4m ∫L0

ϕ 2 (x )dx

,

μ=

∫L0 ϕ (x )dx ∫L0

ϕ 2 (x )dx

,

τ=

e31 bh p (hs + h p ) 1 , θ¯ = Cp R 4C p



le 0

ϕ  (x )dx.

14

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

Appendix B: Derivation of the entries of the Jacobian matrix for the stability analysis of solutions Eqs. (B1)–(B3) are obtained after substituting Eqs. (12) and (13) into (9a) by ignoring the higher harmonic terms and balancing the terms of constant, sin(ωt) and cos(ωt).



a¨ 0 1 + β a0 + 2

β a0 

+

2

β 2

a21

+

b21

 





 3 2 + β a0 a1 a¨ 1 + β a0 b1 b¨ 1 + k¯ 1 a0 + k¯ 2 a0 a20 + a + b21 − 2 1







βω2 a0  2

a21 + b21



a˙ 21 + b˙ 21 + βωa0 b1 a˙ 1 − a1 b˙ 1 + β a˙ 0 a1 a˙ 1 + b1 b˙ 1 + a0 a˙ 0 + 2ξ ωn a˙ 0 = 0



+ a1

3k¯ 2 − 4

 βω2  2

(B1)



  3 2 1 2 1 β a1 + β b1 + β a1 b1 b¨ 1 + (a1 ω − b1 τ ) + k¯ 1 − ω2 a1 4 4 2

2β aa0 a¨ 0 + a¨ 1 1 + β a20 +











a21 + b21 + 3k¯ 2 − βω2 a1 a20 + 2ξ ωn (a˙ 1 − ωb1 ) +







− 2βωa0 b1 a˙ 0 + a0 b˙ 1 + β a˙ 1 2a0 a˙ 0 +



βa  4

4a˙ 20 + a˙ 21 + b˙ 21



  1 ˙ 1 b1 b1 + a1 a˙ 1 − ωb˙ 1 2 + β a21 + β b21 = 0 2 2



(B2)



  1 1 3 β a1 b1 a¨ 1 + b¨ 1 1 + β a20 + β a21 + β b21 + (a1 τ + b1 ω ) + k¯ 1 − ω2 b1 2 4 4

      βb  2  3k¯ 2 βω2  2 + b1 − a1 + b21 + 3k¯ 2 − βω2 b1 a20 + 2ξ ωn b˙ 1 + ωa1 + 4a˙ 0 + a˙ 21 + b˙ 21 4 2 4

2β ba0 a¨ 0 +





  1 1 + 2βωa0 (a1 a˙ 0 + a0 a˙ 1 ) + β b˙ 1 2a0 a˙ 0 + b1 b˙ 1 + a1 a˙ 1 + ωa˙ 1 2 + β a21 + β b21 = −Y μ 2 2

(B3)

Ignoring the second derivatives of the coefficients in Eqs. (B1)–(B3) gives the following equations that only include the first derivatives of the coefficients and themselves.



  β a0  2    3 2 βω2 a0  2 k¯ 1 a0 + k¯ 2 a0 a20 + a + b21 − a1 + b21 + a˙ 1 + b˙ 21 + βωa0 b1 a˙ 1 − a1 b˙ 1 2 1 2 2   + β a˙ 0 a1 a˙ 1 + b1 b˙ 1 + a0 a˙ 0 + 2ξ ωn a˙ 0 = 0

(B4)

      3k¯ 2 βω2  2

(a1 ω − b1 τ ) + k¯ 1 − ω2 a1 + a1 − a1 + b21 + 3k¯ 2 − βω2 a1 a20 + 2ξ ωn (a˙ 1 − ωb1 ) 4

+

( a1 τ

βa 

 2

2





   1 1 4a˙ 20 + a˙ 21 + b˙ 1 − 2βωa0 b1 a˙ 0 + a0 b˙ 1 + β a˙ 1 2a0 a˙ 0 + b1 b˙ 1 + a1 a˙ 1 − ωb˙ 1 2 + β a21 + β b21 = 0 4 2 2 (B5) 







+ b1 ω ) + k¯ 1 − ω2 b1 + b1





+ 2ξ ωn b˙ 1 + ωa1 +



+ β b˙ 1 2a0 a˙ 0 +

βb  4

3k¯ 2 − 4

 βω2  2







a21 + b21 + 3k¯ 2 − βω2 b1 a20



4a˙ 20 + a˙ 21 + b˙ 21 + 2βωa0 (a1 a˙ 0 + a0 a˙ 1 )



  1 ˙ 1 b1 b1 + a1 a˙ 1 + ωa˙ 1 2 + β a21 + β b21 = −Y μ 2 2

(B6)

The goal is to attain the entries of the Jacobian matrix to determine the stability of the HB solutions from the above Eqs. (B4)–(B6). However, it is impossible to directly solve for the entries of the Jacobian matrix from the above equations as a usual way to most general nonlinear problems. Thus, an alternative approach is presented to obtain the Jacobian matrix for the current problem. Taking the partial differentials of the above equations with respect to the coefficient a0 and neglecting ∂ a˙

∂ a˙

∂ b˙

the time derivatives result in the following equation set about ∂ a0 , ∂ a1 , and ∂ a1 . 0 0 0

2ξ ωn

 βω2  2  ∂ a˙ 0 ∂ a˙ 1 ∂ b˙ 1 3 2 + βωa0 b1 − βωa0 a1 = −k¯ 1 − k¯ 2 3a20 + a1 + b21 + a1 + b21 ∂ a0 ∂ a0 ∂ a0 2 2

(B7)

   ∂ a˙ 0 ∂ a˙ 1 ∂ b˙ 1  + 2ξ ωn −ω 2 + β a21 + β b21 + 2β a20 = −2 3k¯ 2 − βω2 a1 a0 ∂ a0 ∂ a0 ∂ a0

(B8)

−2βωa0 b1

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

2βωa0 a1

   ∂ a˙ 0 ∂ a˙ 1  ∂ b˙ 1 +ω 2 + β a21 + β b21 + 2β a20 + 2ξ ωn = −2 3k¯ 2 − βω2 b1 a0 ∂ a0 ∂ a0 ∂ a0 ∂ a˙

15

(B9)

∂ a˙

∂ b˙

Combining the above equations and solving for ∂ a1 , ∂ a1 , and ∂ a0 , one has, 0 0 0

∂ a˙ 1 B2C1 − B1C2 = ∂ a0 A1 B2 − A2 B1

(B10)

∂ b˙ 1 A2C1 − A1C2 = ∂ a0 A1 B2 − A2 B1 

  βω2  2  ∂ a˙ 0 1 B2C1 − B1C2 A2C1 − A1C2 3 2 = −βωa0 b1 − a1 − k¯ 1 − k¯ 2 3a20 + a1 + b21 + a1 + b21 ∂ a0 2ξ ωn A1 B2 − A2 B1 A1 B2 − A2 B1 2 2

(B11)

(B12)

where

A1 = 2(ξ ωn ) + (βωa0 b1 ) 2

2





(B13)





B1 = ξ ωn ω 2 + β a21 + b21 + 2β a20 + (βωa0 ) a1 b1 2





(B14)



  3 2 C1 = −2 3k¯ 2 − βω2 a1 a0 ξ ωn + βωa0 b1 −k¯ 1 − k¯ 2 3a20 + a + b21 + 2 1 





A2 = 2a1 ξ ωn + ωb1 2 + β a21 + b21 + 2β a20











a21 + b21

2

a21 + b21





B2 = −2b1 ξ ωn + ωa1 2 + β a21 + b21 + 2β a20 C2 = −2a0 3k¯ 2 − βω2

βω2 

(B15) (B16)



(B17)



(B18)

Taking the partial differentials of the above Eqs. (B4)–(B6) with respect to the coefficient a1 and neglecting the time ∂ a˙

∂ a˙

∂ b˙

derivatives result in the following equation set about ∂ a0 , ∂ a1 , and ∂ a1 . 1 1 1

2ξ ωn

  ∂ a˙ 0 ∂ a˙ 1 ∂ b˙ 1 + βωa0 b1 − βωa0 a1 = − 3k¯ 2 − βω2 a0 a1 ∂ a1 ∂ a1 ∂ a1

(B19)

     ∂ a˙ 0 ∂ a˙ 1 ∂ b˙ 1  + 2ξ ωn − ω 2 + β a21 + b21 + 2β a20 = − ω − k¯ 1 − ω2 ∂ a1 ∂ a1 ∂ a1

    3k¯ 2 βω2  2 − − 3a1 + b21 − 3k¯ 2 − βω2 a20

− 2βωa0 b1

4

2

  2   ∂ a˙ 0 ∂ a˙ 1  ∂ b˙ 1 3k¯ 2 βω2 2 2 2βωa0 a1 + ω 2 + β a1 + b1 + 2 β a0 + 2 ξ ωn = − τ − 2a1 b1 − − 2ξ ωn ω ∂ a1 ∂ a1 ∂ a1 4 2 ∂ a˙

(B20)

(B21)

∂ a˙

∂ b˙

Similarly, one can solve for ∂ a1 , ∂ a1 , and ∂ a0 as following, 1 1 1

∂ a˙ 1 E2 F1 − E1 F2 = ∂ a1 D1 E2 − D2 E1

(B22)

∂ b˙ 1 D2 F1 − D1 F2 = ∂ a1 D1 E2 − D2 E1  ∂ a˙ 1 0

∂ a1

=

(B23)





− 3k¯ 2 − βω2 a0 a1 − βωa0 b1

2ξ ωn

where







E2 F1 − E1 F2 D2 F1 − D1 F2 − a1 D1 E2 − D2 E1 D1 E2 − D2 E1

D1 = (βωa0 ) a1 b1 − ξ ωn ω 2 + β a21 + b21 + 2β a20 2

E1 = (βωa0 a1 ) + 2(ξ ωn ) 2

2



(B24)

(B25) (B26)

16

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984







F1 = −βω (a0 a1 )2 3k¯ 2 − βω2 + ξ ωn τ + 2ξ ωn a1 b1







D2 = 2a1 ξ ωn + b1 ω 2 + β a21 + b21 + 2β a20







F2 = − (τ b1 + ωa1 ) − 3a1

3k¯ 2 − 4

 βω2  2

βω2



2

+ 2 ( ξ ωn )

2

ω

(B27)



E2 = −2b1 ξ ωn + a1 ω 2 + β a21 + b21 + 2β a20



3k¯ 2 − 4

(B28)



(B29)











a21 + b21 − 2b1 ξ ωn ω − a1 k¯ 1 − ω2 − a1 3k¯ 2 − βω2 a20

(B30)

Taking the partial differentials of the above Eqs. (B4)–(B6) with respect to the coefficient b1 and neglecting the time ∂ a˙

∂ a˙

∂ b˙

derivatives result in the following equation set about ∂ b0 , ∂ b1 , and ∂ b1 . 1 1 1

2ξ ωn

 ∂ a˙ 0 ∂ a˙ 1 ∂ b˙ 1  2 + βωa0 b1 − βωa0 a1 = βω − 3k¯ 2 a0 b1 ∂ b1 ∂ b1 ∂ b1

(B31)

  2   ∂ a˙ 0 ∂ a˙ 1 ∂ b˙ 1  3k¯ 2 2 2 2 −2βωa0 b1 + 2ξ ωn − ω 2 + β a1 + b1 + 2 β a0 = τ − a1 b1 − βω + 2ξ ωn ω ∂ b1 ∂ b1 ∂ b1 2 2βωa0 a1

     ∂ a˙ 0 ∂ a˙ 1  ∂ b˙ + ω 2 + β a21 + b21 + 2β a20 + 2ξ ωn 1 = −ω − k¯ 1 − ω2 ∂ b1 ∂ b1 ∂ b1

    3k¯ 2 βω2  2 − − a + 3b2 − 3k¯ 2 − βω2 a20 4

∂ a˙

∂ a˙

(B32)

(B33)

2

∂ b˙

Solving for ∂ b0 , ∂ b1 , and ∂ b1 , one has, 1 1 1

∂ a˙ 1 H2 M1 − H1 M2 = ∂ b1 G1 H2 − G2 H1

(B34)

∂ b˙ 1 G2 H2 M1 − G1H1 = ∂ b1 G1H2 − G2 H1

(B35)

  H M −H M  ∂ a˙ 0 1 G2 H2 M1 − G1H1 2 1 1 2 = − 3k¯ 2 − βω2 a0 b1 − βωa0 b1 − a1 ∂ b1 2ξ ωn G1 H2 − G2 H1 G1H2 − G2 H1

(B36)

G1 = (βωa0 b1 ) + 2(ξ ωn )

(B37)

where 2

2







H1 = (βωa0 ) a1 b1 + ξ ωn ω 2 + β a21 + b21 + 2β a20 2

M1 = βω (a0 b1 )

2



βω

2





− 3k¯ 2 + ξ ωn





τ − a1 b1



G2 = 2a1 ξ ωn + b1 ω 2 + β a21 + b21 + 2β a20









M2 = (τ a1 − b1 ω ) − b1 k¯ 1 − ω2 − 3b1



(B38)

3k¯ 2 − βω2 2



 + 2ξ ωn ω

(B39)



H2 = −2b1 ξ ωn + a1 ω 2 + β a21 + b21 + 2β a20





(B40)



3k¯ 2 − 4

(B41)

 βω2  2



a21 + b21 −





3k¯ 2 − βω2 a20 b1 + 2a1 ξ ωn ω

(B42)

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

17

References [1] Roundy S, Wright PK, Rabaey J. A study of low level vibrations as a power source for wireless sensor nodes. Comput Commun 2003;26:1131–44. [2] Tang L, Wang J. Size effect of tip mass on performance of cantilevered piezoelectric energy harvester with a dynamic magnifier. Acta Mech 2017;228(11):3997–4015. [3] Liu W, Qin G, Zhu Q, Hu G. Synchronous extraction circuit with self-adaptive peak-detection mechanical switches design for piezoelectric energy harvesting. Appl Energy 2018;230:1292–303. [4] Tang L, Wang J. Modeling and analysis of cantilever piezoelectric energy harvester with a new-type dynamic magnifier. Acta Mech 2018;229(11):4643–62. [5] Zhou W, Penamalli GR, Zuo L. An efficient vibration energy harvester with a multi-mode dynamic magnifier. Smart Mater Struct 2011;21(1):015014. [6] Yan Z, Sun W, Tan T, Huang W. Nonlinear analysis of galloping piezoelectric energy harvesters with inductive-resistive circuits for boundaries of analytical solutions. Commun Nonlinear Sci Numer Simul 2018;62:90–116. [7] Liu J, Fang H, Xu Z, Mao X, Shen X, Chen D, Liao H, Cai B. A MEMS-based piezoelectric power generator array for vibration energy harvesting. Microelectron J 2008;39(5):802–6. [8] Wei C, Jing X. A comprehensive review on vibration energy harvesting: modelling and realization. Renew Sust Eneryg Rev 2017;74:1–18. [9] Fang F, Xia G, Wang J. Nonlinear dynamic analysis of cantilevered piezoelectric energy harvesters under simultaneous parametric and external excitations. Acta Mech Sin 2018;34(3):561–77. [10] Litak G, Friswell MI, Adhikari S. Magnetopiezoelastic energy harvesting driven by random excitations. Appl Phys Lett 2010;96(21):214103. [11] Liu W, Formosa F, Badel A, Hu G. A simplified lumped model for the optimization of post-buckled beam architecture wideband generator. J Sound Vib 2017;24(409):165–79. [12] Tchoukuegno R, Woafo P. Dynamics and active control of motion of a particle in a ϕ ;6 potential with a parametric forcing. Physica D 20 02;167(1):86–10 0. [13] Stanton SC, McGehee CC, Mann BP. Nonlinear dynamics for broadband energy harvesting: investigation of a bistable piezoelectric inertial generator. Physica D 2010;239(10):640–53. [14] Zhou S, Cao J, Erturk A, Lin J. Enhanced broadband piezoelectric energy harvesting using rotatable magnets. App Phys Lett 2013;102(17):173901. [15] Zhou S, Zuo L. Nonlinear dynamic analysis of asymmetric tristable energy harvesters for enhanced energy harvesting. Commun Nonlinear Sci Numer Simul 2018;61:271–84. [16] Tékam GTO, Kwuimy CAK, Woafo P. Analysis of tristable energy harvesting system having fractional order viscoelastic material. Chaos 2015;25(1):013112. [17] Zhou S, Cao J, Inman DJ, Lin J, Liu S, Wang Z. Broadband tristable energy harvester: modeling and experiment verification. Appl Energy 2014;133:33–9. [18] Li H, Qin W, Lan C, Deng W, Zhou Z. Dynamics and coherence resonance of tri-stable energy harvesting system. Smart Mater Struct 2015;25(1):015001. [19] Daqaq MF, Masana R, Erturk A, Quinn DD. On the role of nonlinearities in vibratory energy harvesting: a critical review and discussion. Appl Mech Rev 2014;66(4):040801. [20] Yuan TC, Yang J, Chen LQ. Nonlinear dynamics of a circular piezoelectric plate for vibratory energy harvesting. Commun. Nonlinear Sci Numer Simul 2018;59:651–6. [21] Chen LQ, Jiang WA, Panyam M, Daqaq MF. A broadband internally resonant vibratory energy harvester. J Vib Acoust 2016;138(6):061007. [22] Harne RL, Wang KW. Harnessing bistable structural dynamics: for vibration control, energy harvesting and sensing. Chichester: John Wiley & Sons Ltd.; 2017. [23] Gu L, Livermore C. Impact-driven, frequency up-converting coupled vibration energy harvesting device for low frequency operation. Smart Mater Struct 2011;20(4):045004. [24] Kim IH, Jung HJ, Lee BM, Jang SJ. Broadband energy-harvesting using a two degree-of-freedom vibrating body. App Phys Lett 2011;98(21):214102. [25] Rezaei M, Khadem SE, Firoozy P. Broadband and tunable PZT energy harvesting utilizing local nonlinearity and tip mass effects. Int J Eng Sci 2017;118:1–15. [26] Firoozy P, Khadem SE, Pourkiaee SM. Power enhancement of broadband piezoelectric energy harvesting using a proof mass and nonlinearities in curvature and inertia. Int J Mech Sci 2017;133:227–39. [27] Priya S, Song HC, Zhou Y, Varghese R, Chopra A, Kim SG, Kanno I, Wu L, Ha DS, Ryu J, Polcawich RG. A review on piezoelectric energy harvesting: materials, methods, and circuits. Energy Harvest Syst 2017;4(1):3–39. [28] Chen LQ, Jiang WA. Internal resonance energy harvesting. J Appl Mech 2015;82(3):031004. [29] Nwagoum Tuwa PR, Woafo P. Micro-plate piezoelectric energy harvester for pulsating arterial pressure. J Mech Med Biol 2016;16(05):1650073. [30] Leland ES, Wright PK. Resonance tuning of piezoelectric vibration energy scavenging generators using compressive axial preload. Smart Mater Struct 2006;15(5):1413–20. [31] Xu C, Liang Z, Ren B, Di W, Luo H, Wang D, Wang K, Chen Z. Bi-stable energy harvesting based on a simply supported piezoelectric buckled beam. J Appl Phys 2013;114(11):114507. [32] Andò B, Baglio S, Bulsara AR, Marletta V. A bistable buckled beam based approach for vibrational energy harvesting. Sens Actuators A Phys 2014;211:153–61. [33] Liu W, Badel A, Formosa F, Wu YP, Agbossou A. Novel piezoelectric bistable oscillator architecture for wideband vibration energy harvesting. Smart Mater Struct 2013;22(3):035013. [34] Liu W, Formosa F, Badel A. Optimization study of a piezoelectric bistable generator with doubled voltage frequency using harmonic balance method. J Intell Mater Syst Struct 2017;28(5):671–86. [35] Bayat M, Pakar I. On the approximate analytical solution to non-linear oscillation systems. Shock Vib 2013;20(1):43–52. [36] Tan T, Yan Z, Lei H, Sun W. Geometric nonlinear distributed parameter model for cantilever-beam piezoelectric energy harvesters and structural dimension analysis for galloping mode. J Intell Mater Syst Struct 2017;28(20) 1045389X17704922. [37] Yan Z, Hajj MR. Energy harvesting from an autoparametric vibration absorber. Smart Mater Struct 2015;24(11):115012. [38] Ghayesh MH, Farokhi H, Gholipour A, Hussain S. On the nonlinear mechanics of layered microcantilevers. Int J Eng Sci 2017;120:1–4. [39] Pai PF, Rommel B, Schulz MJ. Non-linear vibration absorbers using higher order internal resonances. J Sound Vib 20 0 0;234(5):799–817. [40] Mahmoodi SN, Jalili N. Non-linear vibrations and frequency response analysis of piezoelectrically driven microcantilevers. Int J Non Linear Mech 2007;42(4):577–87. [41] Yan Z, Taha HE, Tan T. Nonlinear characteristics of an autoparametric vibration system. J Sound Vib 2017;390:1–22. [42] Yan Z, Hajj MR. Nonlinear performances of an autoparametric vibration-based piezoelastic energy harvester. J Intell Mater Syst Struct 2017;28(2):254–71. [43] He JH. Recent development of the homotopy perturbation method. Topol Methods Nonlinear Anal 2008;31(2):205–9. [44] Sedighi HM, Shirazi KH. A new approach to analytical solution of cantilever beam vibration with nonlinear boundary condition. J Comput Nonlinear Dyn 2012;7(3):034502. [45] Guo Z, Zhang W. The spreading residue harmonic balance study on the vibration frequencies of tapered beams. App Math Model 2016;40(15):7195–203. [46] Mohammadian M, Akbarzade M. Higher-order approximate analytical solutions to nonlinear oscillatory systems arising in engineering problems. Arch Appl Mech 2017;87(8):1–6. [47] Akbarzade M, Farshidianfar A. Nonlinear transversely vibrating beams by the improved energy balance method and the global residue harmonic balance method. App Math Model 2017;45:393–404.

18

F. Qian, S. Zhou and L. Zuo / Commun Nonlinear Sci Numer Simulat 80 (2020) 104984

[48] Mahmoodi SN, Afshari M, Jalili N. Nonlinear vibrations of piezoelectric microcantilevers for biologically-induced surface stress sensing. Commun Nonlinear Sci Numer Simul 2008;13(9):1964–77. [49] Tuwa PN, Woafo P. Electromechanical control of the dynamics of a thin elasticplate: analytical method and finite differences simulation. Mech Res Commun 2014;61:19–26. [50] Stanton SC, Owens BA, Mann BP. Harmonic balance analysis of the bistable piezoelectric inertial generator. J Sound Vib 2012;331(15):3617–27. [51] Zhou S, Cao J, Inman DJ, Lin J, Li D. Harmonic balance analysis of nonlinear tristable energy harvesters for performance enhancement. J Sound Vib 2016;373:223–35. [52] Kovacic I, Brennan MJ. The Duffing equation: nonlinear oscillators and their behavior. Chichester: John Wiley & Sons Ltd.; 2011. [53] Xu C, Liang Z, Ren B, Di W, Luo H, Wang D, Wang K, Chen Z. Bi-stable energy harvesting based on a simply supported piezoelectric buckled beam. J Appl Phys 2013;114(11):114507. [54] Qian F, Zhou W, Kaluvan S, Zhang H, Zuo L. Theoretical modeling and experimental validation of a torsional piezoelectric vibration energy harvesting system. Smart Mater Struct 2018;27(4):045018.