Global Dynamics of the Vibrating System of a Tristable Piezoelectric Energy Harvester

Abstract

Global dynamics of a piezoelectric energy harvester with tristable potential is investigated. The dynamical model of a cantilever beam energy harvester is considered; its static bifurcation is also discussed. Multiple intra-well attractors and their basins of attraction are presented to discuss the mechanism of multistability and its initial sensitivity. Moreover, the Melnikov method is applied to present the conditions for global bifurcations and the induced complex dynamics. The results show that the variation of coefficients of the polynomial may affect the number and shapes of potential wells, while the increase of the excitation amplitude may trigger multistability around one equilibrium, initial-sensitive jump, inter-well attractor and chaos. The results may provide some theoretical reference for increasing the working performance of energy harvesters.

1. Introduction

Piezoelectricity is an attractive physical conversion generated by converting mechanical strain to electrical potential [1,2]. Most piezoelectric materials have fast piezoelectric reactions, meaning that they can do mechanical–electrical energy conversion thus harvesting energy from dynamical structures conveniently [3,4]. Owing to this advantage, the profile design of piezoelectric energy harvesters has been paid great attention during these decades. Various structures have been designed for piezoelectric energy harvesters, possessing monostable [5], bistable [6], tristable [7], qua-stable [8] or even quinstable [9] characteristics. Zhou et al. [10] proposed a broadband piezoelectric-based vibration energy harvester with a triple-well potential induced by a magnetic field and presented an experimental investigation, showing that the tristable configuration easily attained higher energy intra-well oscillations. Li et al. [11] considered a multistable piezoelectric energy harvester with a nonlinear spring subjected to wake-galloping and observed that intra-well motion and chaos occurred within a certain range of fluid velocity. Naseer et al. [12] constructed a piezoelectric cantilever-cylinder structure for the sake of energy harvesting from vortex-induced vibration (VIV) and found the energy harvester in the monostable configuration displayed a hardening behavior with higher amplitudes thus a larger output voltage, while in the bistable configuration, it had a wider synchronization region with period or non-period responses but produced a lower output power. Yang et al. [13] designed a magnetic levitation-based hybrid energy harvester and found via quantitative investigation that the tristable system required less kinetic energy to excite a large displacement motion, compared with monostable systems. Wang et al. [14] proposed an ultra-low-frequency energy harvester to harness structural vibration energy and displayed the benefits of multistability for energy harvesting.

It is evident that multistability and chaos show great potential in vibration energy harvesting techniques. The idea is that these two initial-sensitive phenomena, i.e., chaos and jump among intra-well motions and inter-well ones, can easily make a big deformation of the piezoelectric structure, thus collecting much electrical potential. Hence, many researchers have been devoted to investigating these complex phenomena of energy harvesters and their mechanisms. Rezaei et al. [15] exploited the Method of Multiple Scales to provide an approximate-analytical solution for the vibrating system of a nonlinear piezoelectric energy harvester under a hard harmonic excitation. Chen et al. [16] considered an arch-linear-composed-beam piezoelectric energy harvester with magnetic coupling and investigated numerically and experimentally the large-amplitude intra-well motion and chaos. Ju et al. [17] presented numerical results and experiments for a multistable piezoelectric vibration energy harvester with four potential wells, showing that the jump from the inter-well motion to the intra-well motion can be easily triggered under a low acceleration. Cao et al. [18] introduced the fractional model for magnetically coupling broadband energy harvesters under low-frequency excitation and presented the chaotic behavior clearly via numerical simulations and experiments. Lallart et al. [19] exposed the analytical results and numerical ones to provide conditions for the occurrence of multistability in the framework of energy harvesting. Tékam et al. [20] focused on a tristable energy harvesting system having fractional order viscoelastic material and computed its periodic responses by the Krylov–Bogoliubov averaging method. In the dynamical system of a parametrically amplified Mathieu–Duffing nonlinear energy harvester, Karličić et al. [21] obtained an approximation of the periodic response by using the incremental harmonic balance method and exhibited the coexistence of bistable periodic attractors via numerical results. Considering a bistable piezo-magnetoelastic structure for energy harvesting, Barbosa et al. [22] proposed a semi-continuous method to control chaos and presented the control effect by numerical results. Chen et al. [23] proposed Melnikov function-based necessary conditions of chaos in a bistable piezoelectric vibration energy harvesting system and verified their validations numerically. For a novel electromagnetic bistable vibration energy harvester with an elastic boundary, Zhang et al. [24] classified basins of attraction of different attractors and found that multistability will increase the occurring probability of the large-amplitude intra-well responses. Fu et al. [25] modeled a new sliding-mode triboelectric energy harvester in the form of a cantilever beam with a tip mass loaded by both magnetic and friction forces and found three types of multistability in its dynamical system. Sufficient works have studied the complex dynamics of energy harvesting systems in detail, but the mechanism behind multistability, jump and chaos is still not clear yet.

To this end, we consider a type of tristable piezoelectric energy harvester and study the mechanism behind its complex dynamics. The remaining contents are organized as follows: In the next section, the dynamical model is constructed, and its equilibria are discussed. In Section 3, the coexistence of multiple attractors and their mechanisms are discussed in detail. In Section 4, necessary conditions for global bifurcations are proposed and verified by numerical results. Finally, conclusions are discussed in Section 5.

2. Dynamical Model and Its Static Bifurcation

The simplified diagram of the considered tristable piezoelectric energy harvester [10] is shown in Figure 1 where X is the horizontal displacement of the end of the substrate layer at moment t, L the length of the substrate layer. In Figure 1, the magnet at the end of the substrate layer provides nonlinear forces via interacting with the other two magnets; B and D are unstable positions of the vibrating system, A, C and E are stable ones. The energy is stored by the energy harvesting circuits based on the piezoelectric effect of the piezoelectric layers driven by the substrate layer. According to the Second Law of Newton and Kirchhoff’s law, the vibrating system of the energy harvester can be expressed as

$$m\ddot{X}(t)+c\dot{X}(t)+k(X)-k_{em}V(t)=F\cos (\mathsf{\Omega }t),\hspace{1em}{c}_p\dot{V}(t)+\frac{V(t)}{r_l}+k_{em}\dot{X}(t)=0$$

(1)

where m, c, k_em are the equivalent mass, the equivalent damping and equivalent electromechanical coupling coefficient, respectively; c_p is the equivalent capacitance of the piezoelectric materials, $r_l$ the load resistance, V(t) the voltage across the electrical load, X(t) is the tip displacement of the harvester in the transverse direction, k(X) the magnetic force whose nonlinearity due to the effect of magnetic force, and F and Ω are the amplitude and frequency of the external excitation, respectively. The polynomial form for nonlinear magnetic force [15,17] is introduced in order to characterize the relationship between the tip displacement of the cantilever and it below:

$$k(X)=a_{1}X-a_2{X}^3+a_3{X}^5$$

(2)

where $a_1$, $a_2$ and $a_3$ are positive and the coefficients of the polynomial. By introducing the dimensionless time $T={\omega }_{0}t$ where ${\omega }_0=\sqrt{\frac{a_1}{m}}$, and the variables $x(T)=\frac{X(t)}{L},v(T)=\frac{V(t)}{V_0}$, the dimensionless form of Equation (1) can be expressed by

$$x^{″}(T)+\mu x^{\prime }(T)+x(T)-k_1{x}^3(T)+k_2{x}^5(T)-\eta v(T)=f_0\cos \omega T,\hspace{1em}{v}^{\prime }(T)+\gamma v(T)+\beta x^{\prime }(T)=0,$$

(3)

where

$$\omega =\frac{\mathsf{\Omega }}{{\omega }_0},k_1=\frac{a_2}{a_1}{L}^2,k_2=\frac{a_3}{a_1}{L}^4,\mu =\frac{c}{m{\omega }_0},\eta =\frac{k_{em}{V}_0}{L{a}_1},f_0=\frac{F}{L{a}_1},\gamma =\frac{1}{c_p{r}_l{\omega }_0},\beta =\frac{k_{em}L}{c_p{V}_0}.$$

(4)

By denoting $x(T)\triangleq x, v(T)\triangleq v, x^{\prime }(T)\triangleq \dot{x}, x^{″}(T)\triangleq \ddot{x}, v^{\prime }(T)\triangleq \dot{v}$ in Equation (3), the dimensionless state-space model of the piezoelectric vibration energy harvesting system can be obtained as follows:

$$\ddot{x}+\mu \dot{x}+x-k_1{x}^3+k_2{x}^5-\eta v=f_0\cos (\omega T),\hspace{1em}\dot{v}+\gamma v+\beta \dot{x}=0.$$

(5)

Note that all non-dimensional parameters in the above equation are positive. The nomenclatures of the system parameters are presented in

Table 1. Parameters of systems (1) and (5).

Parameter	Symbol
Equivalent mass of the proof mass (kg)	m
Equivalent damping of the piezoelectric beam (N·s/m)	c
Linear stiffness (N/m)	$a_1$
Cubic stiffness term (N/m3)	$a_2$
Penta power stiffness term (N/m5)	$a_3$
Electromechanical coupling coefficient (N/V)	kem
Amplitude of the external excitation (N)	F
Frequency of the external excitation (HZ)	Ω
Equivalent capacitance of piezoelectric layers (F)	cp
Load resistance (kΩ)	$r_l$
Length of substrate layer (mm)	L
Initial voltage of energy harvesting circuits (V)	$V_0$
Time	t
Tip displacement of the harvester at time t	$X(t)$
Voltage at time t	$V(t)$
Natural frequency of the dynamical system	${\omega }_0$
Dimensionless linear damping term	$\mu$
Dimensionless cubic stiffness term	$k_1$
Dimensionless penta power stiffness term	$k_2$
Dimensionless electromechanical coupling term	$\eta$
Dimensionless excitation amplitude	$f_0$
Dimensionless excitation frequency	$\omega$
Dimensionless stiffness term of the coil current system	$\gamma$
Dimensionless electromechanical damping coefficient term	$\beta$
Dimensionless time	T
Dimensionless displacement at time T	x
Dimensionless voltage at time T	v

.

Assuming $\mu$ = 0, η = 0, and f₀ = 0 in Equation (5) yields its unperturbed system

$$\dot{x}=y,\hspace{1em}\dot{y}=-x+k_1{x}^3-k_2{x}^5.$$

(6)

which is a Hamilton system. Letting the right side of Equation (6) be zero, one obtains its equilibria whose vertical coordinate y is zero and horizontal coordinates satisfy

$$x(1-k_1{x}^2+k_2{x}^4)=0.$$

(7)

Their stability is determined by the roots of the following characteristic equation

$${\lambda }^2+(1-3k_1{x}^2+5k_2{x}^4)=0.$$

(8)

Obviously, the origin O(0, 0) is the equilibrium of the system (6) whose eigenvalue is $\lambda =\pm i$, implying that O(0, 0) is a center. The number of the nontrivial equilibria and the shapes of possible potential wells of the unperturbed system (6) depend on the non-dimensional parameters k₁ and k₂. For example, for $k_1{}^2<4k_2$, there is no nontrivial equilibria in Equation (7), because according to Weda’s Theorem, apart from zero, there are no real roots of $x^2$ in Equation (7). For $k_1{}^2>4k_2$, We have the following theorem.

Theorem 1.

If $k_1{}^2>4k_2$, there will be four nontrivial equilibria of Equation (7), i.e., two centers and two saddles.

Proof of Theorem 1.

If $k_1{}^2>4k_2$, according to Weda’s Theorem, there are two pairs of positive roots of $x^2$ in Equation (7) expressed by $\frac{k_1\pm \sqrt{\Delta }}{2k_2}$ where $\Delta =k_1{}^2-4k_2>0$; thus, there are four real solutions of x in Equation (7), namely, $\pm \sqrt{\frac{k_1+\sqrt{\Delta }}{2k_2}}$ and $\pm \sqrt{\frac{k_1-\sqrt{\Delta }}{2k_2}}$.

In the neighborhood of the two equilibria $S_{\pm }(\pm \sqrt{\frac{k_1-\sqrt{\Delta }}{2k_2}},0)$, the two eigenvalues solved from Equation (8) are a positive one and a negative one, expressed by $\lambda =\pm \sqrt{\frac{\sqrt{\Delta }(k_1-\sqrt{\Delta })}{k_2}}$, respectively. It means that $S_{\pm }(\pm \sqrt{\frac{k_1-\sqrt{\Delta }}{2k_2}},0)$ are saddles. In the neighborhood of the other two nontrivial equilibria $C_{\pm }(\pm \sqrt{\frac{k_1+\sqrt{\Delta }}{2k_2}},0)$, the characteristic equation becomes

$${\lambda }^2+\frac{\sqrt{\Delta }(k_1+\sqrt{\Delta })}{k_2}=0.$$

(9)

As $\frac{\sqrt{\Delta }(k_1+\sqrt{\Delta })}{k_2}>0$, there are two pure imaginary roots in the above equation, illustrating that the equilibria $C_{\pm }(\pm \sqrt{\frac{k_1+\sqrt{\Delta }}{2k_2}},0)$ are centers. □

Accordingly, the critical condition for static bifurcation of the equilibria is $k_1{}^2=4k_2$. The different types of orbits of the unperturbed system (6) are classified in k₁-k₂ plane, as shown in Figure 2. It follows that for a given value of the dimensionless cubic stiffness term k₁, the lower the dimensionless penta power stiffness term k₂ is, the higher the probability of a triple well the system will have. Based on the relationship between the dimensionless parameters and the original parameters shown in Equation (4), it indicates that the lower the penta power stiffness coefficient of polynomial magnetic force $a_3$ is, the higher possibility for the occurrence of a triple potential well will be. For the case of multiple equilibria, two nontrivial equilibria $C_{\pm }(\pm \sqrt{\frac{k_1+\sqrt{\Delta }}{2k_2}},0)$ are the centers of two potential wells surrounded by homoclinic orbits, the origin is the center of a well surrounded by heteroclinic orbits crossing the other two nontrivial saddles $S_{\pm }(\pm \sqrt{\frac{k_1-\sqrt{\Delta }}{2k_2}},0)$. As well known, multiple wells may induce multistability [25,26]. In the following sections, we consider the case of three potential energy wells, and focus on how to make use of the external excitation and initial conditions to induce complex dynamics such as jump phenomena among periodic attractors, inter-well oscillation or even chaos, thus harvesting vibration energy effectively. All the values of system parameters are dimensionless for analysis. Based on the physical properties of the energy harvester in reference [10], some invariable parameters can be set as:

$$k_1=8,k_2=12,\mu =0.1,\eta =0.1,\gamma =0.1,\beta =0.1.$$

(10)

The parameters $f_0$ and $\omega$ will be changed to study the influence mechanism of dynamical response characteristics. It can be calculated that the horizontal coordinates for the two nontrivial centers $C_{\pm }(\pm x_c,0)$ and saddles $S_{\pm }(\pm x_s,0)$ are $x_c=\pm \frac{\sqrt{2}}{2}$ and $x_s=\pm \frac{\sqrt{6}}{6}$, respectively.

3. Periodic Responses

3.1. Periodic Solutions near the Origin (0, 0)

In this subsection, we consider the periodic solutions near the equilibrium (0, 0). First, by introducing a small parameter ε satisfying 0 < ε≪1 and rescaling the parameters $\mu ,\gamma ,\eta$ and $f_0$ in Equation (5) as

$$\mu =\epsilon \tilde{\mu },\gamma =\epsilon \tilde{\gamma },\eta =\epsilon \tilde{\eta },f_0={\epsilon }^2\tilde{f},$$

(11)

one can rewrite the system (5) as

$$\ddot{x}+x=-\epsilon \tilde{\mu }\dot{x}+{\epsilon }^2\tilde{f}\cos \omega T+\epsilon \tilde{\eta }v+k_1{x}^3-k_2{x}^5,\dot{v}+\epsilon \tilde{\gamma }v=-\beta \dot{x}.$$

(12)

The Method of Multiple Scales (MMS) [15] is employed to obtain the approximate-analytical solution of Equation (12). To this end, the displacement x and voltage response v are expanded as

$$x=\epsilon x_1+{\epsilon }^2{x}_2+{\epsilon }^3{x}_3+\cdots , v=\epsilon v_1+{\epsilon }^2{v}_2+{\epsilon }^3{v}_3+\cdots .$$

(13)

The time derivatives can be rewritten as

$$T_i={\epsilon }^{i}T,D_i=\frac{\partial }{\partial T_i},\frac{d}{dT}={\displaystyle \sum _{i=0}^n{\epsilon }^i{D}^i}(i=0,1,\cdots ).$$

(14)

Considering harmonic resonance, it is assumed that

$$\omega =1+\epsilon \sigma$$

(15)

where the detuning parameter $\sigma =O(1)$. Introducing Equations (13)–(15) into Equation (12) and separating the coefficients of ε¹, ε² and ε³ leads to

$${\epsilon }^1:D_0{}^2{x}_1+{\omega }^2{x}_1=0,D_0{v}_1=-\beta D_0{x}_1,$$

(16)

$${\epsilon }^2:D_0{}^2{x}_2+{\omega }^2{x}_2=-2D_0{D}_1{x}_1-\tilde{\mu }{D}_0{x}_1+\tilde{f}\cos \omega T_0+\tilde{\eta }{v}_1+2\sigma \omega x_1,D_0{v}_2=-\tilde{\gamma }{v}_1-\beta (D_1{x}_1+D_0{x}_2)-D_1{v}_1,$$

(17)

and

$${\epsilon }^3:D_0{}^2{x}_3+{\omega }^2{x}_3=-D_1{}^2{x}_1-2D_2{D}_0{x}_1-2D_1{D}_0{x}_2-\tilde{\mu }(D_1{x}_1+D_0{x}_2)+\tilde{\eta }{v}_2+2\sigma \omega x_2-x_1{\sigma }^2+k_1{x}_1{}^3.$$

(18)

The solution of Equation (16) can be assumed as

$$x_1=A(T_1,T_2)e^{i\omega T_0}+cc,\hspace{1em}{v}_1=-\beta A{e}^{i\omega T_0}+cc$$

(19)

where cc represents the complex conjugate of the preceding terms, and the complex slow-times dependent amplitude $A(T_1,T_2)$ is expressed as $A=\frac{a(T_1,T_2)}{2}{e}^{i\theta }$. Here, the $a(T_1,T_2)$ and $\theta$ are real, representing the real amplitude and the phase difference of the periodic solution $x_1$, respectively. By substituting the solution (19) into Equation (17) and separating the secular terms, one has

$$D_{1}A=-\frac{\tilde{\mu }A}{2}-\frac{i\tilde{f}}{4\omega }+\frac{\tilde{\eta }\beta Ai}{2\omega }-i\sigma A.$$

(20)

Since $A$ is the function of $T_1$ and $T_2$, it can be solved from Equation (17) that

$$x_2=0, v_2=-\frac{\beta \tilde{\gamma }Ai}{\omega }+cc.$$

(21)

Submitting Equations (19) and (21) into Equation (13), one can approximately express the periodic solution x as $x=\epsilon a\cos (\omega T_0+\theta )$. Substituting Equations (19)–(21) into Equation (18) and separating its secular terms yields

$$D_{2}A=\frac{\tilde{f}\beta \tilde{\eta }i}{16{\omega }^3}-\frac{A{\beta }^2{\tilde{\eta }}^{2}i}{8{\omega }^3}-\frac{A\beta \tilde{\eta }\tilde{\gamma }}{2{\omega }^2}+\frac{\tilde{f}\tilde{\mu }}{16{\omega }^2}-\frac{\tilde{f}\sigma i}{8{\omega }^2}+\frac{A\beta \tilde{\eta }\sigma i}{2{\omega }^2}-\frac{3k_1{A}^2\overline{A}}{2\omega }-\frac{A{\tilde{\mu }}^{2}i}{8\omega }.$$

(22)

Letting the amplitude of x be $\widehat{a}=\epsilon a$, considering $\dot{A}\approx D_{0}A+\epsilon D_{1}A+{\epsilon }^2{D}_{2}A$, substituting Equations (20) and (22) into it and returning its parameters to the original dimensionless parameters of Equation (5), one obtains

$$\begin{array}{l}\dot{\widehat{a}}=-\frac{\mu }{2}\widehat{a}-\frac{\eta \beta \gamma \widehat{a}}{2{\omega }^2}+\frac{f_0\mu \cos \theta }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega )f_0\sin \theta }{8{\omega }^3},\\ \widehat{a}\dot{\theta }=-(\omega -1)\widehat{a}-\frac{{\eta }^2{\beta }^2\widehat{a}}{8{\omega }^3}+\frac{\eta \beta (\omega -1)\widehat{a}}{2{\omega }^2}-\frac{3k_1{\widehat{a}}^3}{8\omega }+\frac{\eta \beta \widehat{a}}{2\omega }-\frac{{\mu }^2\widehat{a}}{8\omega }-\frac{f_0\mu \sin \theta }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega )f_0\cos \theta }{8{\omega }^3}\end{array}$$

(23)

Letting the right side of Equation (23) be zero leads to

$$\begin{array}{l}\frac{\mu }{2}\widehat{a}+\frac{\eta \beta \gamma \widehat{a}}{2{\omega }^2}=\frac{f_0\mu \cos \theta }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega )f_0\sin \theta }{8{\omega }^3},\\ (\omega -1)\widehat{a}+\frac{{\eta }^2{\beta }^2\widehat{a}}{8{\omega }^3}-\frac{\eta \beta (\omega -1)\widehat{a}}{2{\omega }^2}+\frac{3k_1{\widehat{a}}^3}{8\omega }-\frac{\eta \beta \widehat{a}}{2\omega }+\frac{{\mu }^2\widehat{a}}{8\omega }=-\frac{f_0\mu \sin \theta }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega )f_0\cos \theta }{8{\omega }^3}.\end{array}$$

(24)

Eliminating the triangulation function of Equation (24), one can get

$$\frac{f_0{}^2{\mu }^2}{64{\omega }^4}+\frac{{(\beta \eta -6{\omega }^2+2\omega )}^2{f}_0{}^2}{64{\omega }^6}={(\frac{\mu }{2}+\frac{\eta \beta \gamma }{2{\omega }^2})}^2{\widehat{a}}^2+{(\omega -1+\frac{{\eta }^2{\beta }^2}{8{\omega }^3}-\frac{\eta \beta (\omega -1)}{2{\omega }^2}+\frac{3k_1{\widehat{a}}^2}{8\omega }-\frac{\eta \beta }{2\omega }+\frac{{\mu }^2}{8\omega })}^2{\widehat{a}}^2.$$

(25)

To determine the stability of the periodic solution, one can get its characteristic equation. Based on Equations (23)–(25), the characteristic equation can be written as

$${\lambda }^2+(\mu +\frac{\eta \beta \gamma }{{\omega }^2})\lambda +\frac{f_0{}^2{\mu }^2}{64{\omega }^4}+\frac{{(\beta \eta -6{\omega }^2+2\omega )}^2{f}_0{}^2}{64{\omega }^6}+((\omega -1)\widehat{a}+\frac{{\eta }^2{\beta }^2\widehat{a}}{8{\omega }^3}-\frac{\eta \beta (\omega -1)\widehat{a}}{2{\omega }^2}+\frac{3k_1{\widehat{a}}^3}{8\omega }-\frac{\eta \beta \widehat{a}}{2\omega }+\frac{{\mu }^2\widehat{a}}{8\omega })\frac{3k_1\widehat{a}}{4\omega }=0$$

(26)

According to the equation above, the stability of the periodic solution in the neighborhood of the origin will be changed when λ = 0, namely

$$\frac{f_0{}^2{\mu }^2}{64{\omega }^4}+\frac{{(\beta \eta -6{\omega }^2+2\omega )}^2{f}_0{}^2}{64{\omega }^6}+((\omega -1)\widehat{a}+\frac{{\eta }^2{\beta }^2\widehat{a}}{8{\omega }^3}-\frac{\eta \beta (\omega -1)\widehat{a}}{2{\omega }^2}+\frac{3k_1{\widehat{a}}^3}{8\omega }-\frac{\eta \beta \widehat{a}}{2\omega }+\frac{{\mu }^2\widehat{a}}{8\omega })\frac{3k_1\widehat{a}}{4\omega }=0.$$

(27)

The frequency response curves in the neighborhood of the origin are plotted in Figure 3. In Figure 3, the branches of analytical amplitude $\widehat{a}$ for f₀ = 0.002, 0.01 and 0.05 are obtained from Equation (25). The solid curves and the dashing ones represent the stable and unstable periodic solutions, respectively. Their stability is determined by the positive-negative sign of the real parts of eigenvalues solved from Equation (26). Based on Equation (27), one gets the saddle-node point separating the solid curve and dashing curve. The numerical results for the amplitudes of periodic responses are presented by applying the 4th Runge-Kutta approach to simulate the numerical solutions of the dimensionless system (5) via MATLAB. Definitely, the numerical approach will only produce stable solutions. As illustrated in Figure 3, there is a good agreement between the theoretical and numerical response results, implying the approximate-analytical solutions are valid. Furthermore, it follows from Figure 3 that for f₀ = 0.05 and ω in the range (0.50, 0.76), the frequency response curves bend to the left and yield multivalued solutions. It means that the saddle-node bifurcation of the periodic solution around the origin O(0, 0), hence the coexistence of two intra-well attractors, can be triggered by the increase of f₀ and the variation of ω in a low-frequency range.

When two intra-well attractors around the well center O(0, 0) or $C_{\pm }(\pm x_c,0)$ is induced in the vicinity of the saddle-node bifurcation point, the initial conditions will determine the branch to which the response attracts. For the purpose of energy harvesting, the initial conditions should be chosen to follow the high-amplitude one in order to output high voltage consequently. That is why we present the sequences of the attractors and their basins of attraction for the variation of the dimensionless external excitation. Basin of attraction means the union of initial conditions leading to the same attractor [20,21]. If the boundary of the basin of attraction of one attractor is fractal and intermingled with another, a jump among multiple attractors may easily occur. In this study, the 4th Runge–Kutta approach and the point-mapping method [26] are employed to describe the basin of attraction. The basins of attraction are drawn in the initial-condition plane −1.0 ≤ x(0) ≤ 1.0, −1.0 ≤ y(0) ≤ 1.0 by generating a 600 × 600 array of starting conditions, for each of the initial points. The time step is taken as 0.01. For each figure of sequences of attractors and their basins of attraction, there are pairs of pictures on each line of which the left one is the phase map of the attractors, while the right one shows their basins of attraction. For each attractor, its basin of attraction is marked in the same color as its phase map.

Given $f_0=0.05$, the sequences of attractors and their basins of attraction with the variation of the dimensionless excitation frequency ω are depicted in Figure 4. For ω = 0.3 (see Figure 4(a1)), there coexist three intra-well attractors around three different well centers O(0, 0) and $C_{\pm }(\pm x_c,0)$, respectively. Even though their basins of attraction in some areas entangle each other (see Figure 4(a2)), the vicinity of three potential-well centers is single-colored with clear basin boundaries, meaning that near the potential-well centers, the phenomenon jump will not occur. Only if the initial condition of the structure changes dramatically from the neighborhood of one center to another, there will be a jump between two intra-well attractors. It is worth mentioning that the attractors are around different potential-well centers, thus jump is an inter-well jump. The amplitudes of the three attractors are very low, contributing little to harvest energy. Comparatively, the inter-well jump caused by the change of initial conditions contributes more to energy harvesting. As ω increases to 0.55, there are four intra-well attractors coexisting (see Figure 4(b1)). Apart from the three intra-well attractors, a new intra-well attractor around O(0, 0) with a much higher amplitude appears. However, as can be observed in Figure 4(b2), its basin of attraction is fractal and hard to be detected. Hence, it is a so-called rare attractor [27] whose occurring probability is very low. Compared with the basin map in Figure 4(a2), a better point is that outside of the vicinity of the three well centers, the fractality extent of the basins of attraction becomes more severe, meaning that jump can be triggered more easily in this region. When ω grows to 0.75, the occurrence probability of the higher-amplitude intra-well attractor around O(0, 0) is much higher, because its basin of attraction grows larger, especially in the neighborhood of O(0, 0) (see Figure 4(c2)). Note that the phenomenon of jump between the two intra-well attractors around O(0, 0), namely the intra-well jump, more easily occurs than the inter-well jump, which will be beneficial for energy harvesting. As ω increases to 1.0, the lower-amplitude intra-well attractor around O(0, 0) disappears; still there coexist three attractors, as shown in Figure 4(d1). Based on their basins of attraction in Figure 4(d2), it is obvious that a better energy harvesting performance can be given by the higher-amplitude intra-well oscillation around O(0, 0) or the inter-well jump led by the dramatic change of initial conditions.

3.2. Periodic Solutions near the Nontrivial Equilibria ($\pm x_c$, 0)

Around the nontrivial equilibria ($\pm x_c$, 0), the periodic vibration of the cantilever beam structure can be induced by the perturbation of the two nontrivial centers of the unperturbed system (6) (see C₋ and C₊ in Figure 2). To begin with, supposing $x=\pm x_c+\widehat{x}$, and rescaling the system parameters by Equation (11) yields

$$\ddot{\widehat{x}}+{\widehat{\omega }}^2\widehat{x}=-\epsilon \tilde{\mu }\dot{\widehat{x}}+{\epsilon }^2\tilde{f}\cos \omega T+\epsilon \tilde{\eta }v\mp Q_1{\widehat{x}}^2-Q_2{\widehat{x}}^3-k_2(\pm 5x_c{\widehat{x}}^4+{\widehat{x}}^5),\hspace{1em}\dot{v}+\epsilon \tilde{\gamma }v=-\beta \dot{\widehat{x}}.$$

(28)

where

$${\widehat{\omega }}^2=1-3k_1{x}_c{}^2+5k_2{x}_c{}^4, Q_1=10k_2{x}_c{}^3-3k_1{x}_c, Q_2=10k_2{x}_c{}^2-k_1.$$

(29)

It can be calculated that $\widehat{\omega }=2$, $Q_1=25.46$ and $Q_2=52$, all positive. To apply the Method of Multiple Scales, assuming harmonic resonance in the system (28) that

$$\omega =\widehat{\omega }+\epsilon \widehat{\sigma },$$

(30)

where $\widehat{\sigma }$ is the detuning parameter, rescaling $\widehat{x}$ and voltage response $v$ in Equation (28) as

$$\widehat{x}=\epsilon {\widehat{x}}_1+{\epsilon }^2{\widehat{x}}_2++{\epsilon }^3{\widehat{x}}_3\cdots ,\hspace{1em}v=\epsilon {\widehat{v}}_1+{\epsilon }^2{\widehat{v}}_2+{\epsilon }^3{\widehat{v}}_3\cdots ,$$

(31)

and comparing the coefficients of ε¹, ε² and ε³ in the system (28), respectively, one has

$${\epsilon }^1:D_0{}^2{\widehat{x}}_1+{\omega }^2{\widehat{x}}_1=0,D_0{\widehat{v}}_1=-\beta D_0{\widehat{x}}_1,$$

(32)

$$\begin{array}{c}{\epsilon }^2:D_0{}^2{\widehat{x}}_2+{\omega }^2{\widehat{x}}_2=-2D_0{D}_1{\widehat{x}}_1-\tilde{\mu }{D}_0{\widehat{x}}_1+\tilde{f}\cos \omega T_0+\tilde{\eta }{\widehat{v}}_1+2\sigma \omega {\widehat{x}}_1\mp Q_1{\widehat{x}}_1{}^2,\\ D_0{\widehat{v}}_2=-\tilde{\gamma }{\widehat{v}}_1-\beta (D_1{\widehat{x}}_1+D_0{\widehat{x}}_2)-D_1{\widehat{v}}_1,\end{array}$$

(33)

and

$${\epsilon }^3:D_0{}^2{\widehat{x}}_3+{\omega }^2{\widehat{x}}_3=-D_1{}^2{\widehat{x}}_1-2D_2{D}_0{\widehat{x}}_1-2D_1{D}_0{\widehat{x}}_2-\tilde{\mu }(D_1{\widehat{x}}_1+D_0{\widehat{x}}_2)+\tilde{\eta }{\widehat{v}}_2+2\sigma \omega {\widehat{x}}_2-{\widehat{x}}_1{\sigma }^2\mp 2Q_1{\widehat{x}}_1{\widehat{x}}_2-Q_2{\widehat{x}}_1{}^3.$$

(34)

The solution of Equation (32) can be written as

$${\widehat{x}}_1=B(T_1,T_2)e^{i\omega T_0}+cc,\hspace{1em}{\widehat{v}}_1=-\beta B{e}^{i\omega T_0}+cc,$$

(35)

where the complex slow-times dependent amplitude $B(T_1,T_2)$ is given by

$$B=\frac{b(T_1,T_2)}{2}{e}^{i\phi (T_1,T_2)}.$$

(36)

In Equation (36), $b(T_1,T_2)$ and $\phi (T_1,T_2)$ are real, representing the real amplitude and the phase difference of the periodic solution ${\widehat{x}}_1$, respectively. Substituting Equation (35) into Equation (33) and separating its secular terms leads to

$$D_{1}B=\frac{-\tilde{\mu }B}{2}-\frac{i\tilde{f}}{4\omega }+\frac{\tilde{\eta }\beta Bi}{2\omega }-i\sigma B,$$

(37)

and

$${\widehat{x}}_2=\pm (-\frac{2Q_{1}B\overline{B}}{{\omega }^2}+\frac{Q_1{B}^2}{3{\omega }^2}{e}^{2i\omega T_0}+\frac{Q_1{\overline{B}}^2}{3{\omega }^2}{e}^{-2i\omega T_0}),{\widehat{v}}_2=\mp \frac{\beta Q_1{B}^2}{3{\omega }^2}{e}^{2i\omega T_0}-\frac{\tilde{\gamma }\beta Bi}{\omega }{e}^{i\omega T_0}+cc,$$

(38)

substituting Equations (35), (37) and (38) into Equation (34) and separating the secular terms there yields

$$D_{2}B=\frac{\tilde{f}\beta \tilde{\eta }i}{16{\omega }^3}-\frac{B{\beta }^2{\tilde{\eta }}^{2}i}{8{\omega }^3}-\frac{B\beta \tilde{\eta }\tilde{\gamma }}{2{\omega }^2}+\frac{\tilde{f}\tilde{\mu }}{16{\omega }^2}-\frac{\tilde{f}\sigma i}{8{\omega }^2}+\frac{B\beta \tilde{\eta }\sigma i}{2{\omega }^2}+\frac{3Q_2{B}^2\overline{B}i}{2\omega }-\frac{5Q_1{}^2{B}^2\overline{B}i}{3{\omega }^3}-\frac{B{\tilde{\mu }}^{2}i}{8\omega }.$$

(39)

Considering $\dot{B}\approx D_{0}B+\epsilon D_{1}B+{\epsilon }^2{D}_{2}B$ and $\widehat{b}=\epsilon b$, substituting Equations (36), (37) and (39) into it and expressing it by the original dimensionless parameters of Equation (5), one has:

$$\begin{array}{l}\dot{\widehat{b}}=-\frac{\mu }{2}\widehat{b}-\frac{\eta \beta \gamma \widehat{b}}{2{\omega }^2}+\frac{f_0\mu \cos \phi }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })f_0\sin \phi }{8{\omega }^3},\\ \widehat{b}\dot{\phi }=-(\omega -\widehat{\omega })\widehat{b}-\frac{{\eta }^2{\beta }^2\widehat{b}}{8{\omega }^3}+\frac{\eta \beta (\omega -\widehat{\omega })\widehat{b}}{2{\omega }^2}+\frac{3Q_2{\widehat{b}}^3}{8\omega }-\frac{5Q_1{}^2{\widehat{b}}^3}{12{\omega }^3}+\frac{\eta \beta \widehat{b}}{2\omega }-\frac{{\mu }^2\widehat{b}}{8\omega }-\frac{f_0\mu \sin \phi }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })f_0\cos \phi }{8{\omega }^3}\end{array}$$

(40)

For the periodic solution around C₋ or C₊, the right side of Equation (40) is zero, i.e.,

$$\begin{array}{l}\frac{\mu }{2}\widehat{b}+\frac{\eta \beta \gamma \widehat{b}}{2{\omega }^2}=\frac{f_0\mu \cos \theta }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })f_0\sin \theta }{8{\omega }^3},\\ (\omega -\widehat{\omega })\widehat{b}+\frac{{\eta }^2{\beta }^2\widehat{b}}{8{\omega }^3}-\frac{\eta \beta (\omega -\widehat{\omega })\widehat{b}}{2{\omega }^2}-\frac{3Q_2{\widehat{b}}^3}{8\omega }+\frac{5Q_1{}^2{\widehat{b}}^3}{12{\omega }^3}-\frac{\eta \beta \widehat{b}}{2\omega }+\frac{{\mu }^2\widehat{b}}{8\omega }=-\frac{f_0\mu \sin \theta }{8{\omega }^2}+\frac{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })f_0\cos \theta }{8{\omega }^3}.\end{array}$$

(41)

Eliminating the triangulation function of Equation (41), one obtains

$$\frac{f_0{}^2{\mu }^2}{64{\omega }^4}+\frac{{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })}^2{f}_0{}^2}{64{\omega }^6}={(\frac{\mu }{2}+\frac{\eta \beta \gamma }{2{\omega }^2})}^2{\widehat{b}}^2+{(\omega -\widehat{\omega }+\frac{{\eta }^2{\beta }^2}{8{\omega }^3}-\frac{\eta \beta (\omega -\widehat{\omega })}{2{\omega }^2}-\frac{3Q_2{\widehat{b}}^2}{8\omega }+\frac{5Q_1{}^2{\widehat{b}}^2}{12{\omega }^3}-\frac{\eta \beta }{2\omega }+\frac{{\mu }^2}{8\omega })}^2{\widehat{b}}^2.$$

(42)

The periodic solution near the nontrivial equilibria ($\pm x_c$, 0) can be expressed as

$$x=\pm x_c\mp \frac{2Q_1{\widehat{b}}^2}{3{\omega }^2}+\widehat{b}\cos (\omega T+\phi )\pm \frac{Q_1{\widehat{b}}^2}{3{\omega }^2}{\cos }^2(\omega T+\phi )).$$

(43)

As in Section 3.1, the stability of the periodic solutions can be determined by the following characteristic equation

$$\begin{array}{l}{\lambda }^2+(\mu +\frac{\eta \beta \gamma }{{\omega }^2})\lambda +\frac{f_0{}^2{\mu }^2}{64{\omega }^4}+\frac{{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })}^2{f}_0{}^2}{64{\omega }^6}\\ +((\omega -\widehat{\omega })\widehat{b}+\frac{{\eta }^2{\beta }^2\widehat{b}}{8{\omega }^3}-\frac{\eta \beta (\omega -\widehat{\omega })\widehat{b}}{2{\omega }^2}-\frac{3Q_2{\widehat{b}}^3}{8\omega }+\frac{5Q_1{}^2{\widehat{b}}^3}{12{\omega }^3}-\frac{\eta \beta \widehat{b}}{2\omega }+\frac{{\mu }^2\widehat{b}}{8\omega })(\frac{5Q_1{}^2\widehat{b}}{6{\omega }^3}-\frac{3Q_2\widehat{b}}{4\omega })=0.\end{array}$$

(44)

Hence, saddle-node bifurcation will occur when

$$\frac{f_0{}^2{\mu }^2}{64{\omega }^4}+\frac{{(\beta \eta -6{\omega }^2+2\omega \widehat{\omega })}^2{f}_0{}^2}{64{\omega }^6}+((\omega -\widehat{\omega })\widehat{b}+\frac{{\eta }^2{\beta }^2\widehat{b}}{8{\omega }^3}-\frac{\eta \beta (\omega -\widehat{\omega })\widehat{b}}{2{\omega }^2}-\frac{3Q_2{\widehat{b}}^3}{8\omega }+\frac{5Q_1{}^2{\widehat{b}}^3}{12{\omega }^3}-\frac{\eta \beta \widehat{b}}{2\omega }+\frac{{\mu }^2\widehat{b}}{8\omega })(\frac{5Q_1{}^2\widehat{b}}{6{\omega }^3}-\frac{3Q_2\widehat{b}}{4\omega })=0$$

(45)

According to Equations (42)–(45), the frequency response curves in the neighborhood of the nontrivial equilibria are presented in Figure 5 where the branches of analytical amplitude $\widehat{b}$ for f₀ = 0.002, 0.01 and 0.05 are obtained from Equation (42). The solid curves and the dashed ones show the stable and unstable periodic solutions, respectively. Their stability is determined by the positive-negative sign of the real parts of eigenvalues solved from Equation (44). The saddle-node point separating the solid curve and dashing curve is determined by Equation (45). By applying 4th Runge–Kutta approach to simulate the numerical solutions of the dimensionless system (5), the numerical results are presented. As shown in Figure 5, theoretical results totally agree with the numerical ones. Similar to Figure 3, it also can be observed in Figure 5 that for f₀ = 0.05, the frequency response curves bend to the left and yield bistability, indicating that due to saddle-node bifurcation of the periodic solution, in the vicinity of the saddle-node bifurcation point, there are two intra-well periodic attractors around each nontrivial equilibrium.

Based on the analysis of the periodic responses near the origin and the nontrivial equilibria, the evolution of the periodic solutions of system (5) with the amplitude of f₀ for ω = 1.5 is illustrated in Figure 6. Obviously, when f₀ increases from 0, the nontrivial equilibria and the origin lose their stability; instead, three periodic attractors coexist, attributed to the disturbance of the tristable equilibria. When f₀ varies from 0 to 0.02, the numerical simulation is in great agreement with the theoretical solution. However, when f₀ continues to increase, the theoretical prediction of the periodic responses does not match the numerical results well, which may be due to the limitation of the Method of Multiple Scale, namely the value of the excitation frequency ω should be closed to the inherent frequency of the perturbed systems (12) and (28) near the origin and the nontrivial equilibria ($\pm x_c$, 0). Actually, the given value of the excitation frequency is ω = 1.5, while the inherent frequencies of the systems (12) and (28) are 1 and 2, respectively. It means that the derivations of ω from the two inherent frequencies are not that small. Although the prediction is not that quantitatively accurate, the theoretical results and numerical simulation both show that the increase of the excitation amplitude triggers bistability around each nontrivial well center ($\pm x_c$, 0). According to the theoretical analysis, it can be ascribed to the saddle-node bifurcation of the periodic solutions. For instance, for f₀ = 0.05, there coexist five intra-well periodic attractors, i.e., two around each nontrivial potential-well center and one around O(0, 0).

Fixing ω = 1.5, the evolution of periodic responses and their basins of attraction with the increase of the dimensionless excitation amplitude f₀ is shown in Figure 7. For f₀ = 0, the coexisting triple attractors are the three well centers whose basin boundaries are smooth, as shown in Figure 7(a1,a2). For f₀ = 0.045, the three point attractors are perturbed and replaced by three intra-well periodic attractors around them (see Figure 7(b1)). Although the single-colored neighborhood of each well center in Figure 7(b2) implies that the initial conditions chosen there must attract to the intra-well attractor around it, the inter-well jump may be triggered when the initial conditions are disturbed dramatically. As f₀ increases a bit to reach 0.046, two new intra-well periodic attractors with higher amplitudes suddenly appear in the neighborhood of the nontrivial well centers. As predicted in Figure 6, this can be attributed to the saddle-node bifurcation. Now, there are five intra-well attractors in Figure 7(c1). However, as shown in Figure 4(c2), the basins of attraction of the new attractors are fractal and scattered. Thus, the occurring probability for these two high-amplitude intra-well responses is very low. As f₀ continues to increase, the basins of attraction of the two new periodic attractors are enlarged (see Figure 7(d2)). It means that the occurrence probability of the higher-amplitude intra-well attractors around the well centers C₋ and C₊ is much higher. Moreover, the jump between the two intra-well attractors around each nontrivial well center, namely the intra-well jump, is of higher probability to occur. Therefore, not only the intra-well jump due to saddle-node bifurcation of the periodic solution around the nontrivial well centers and a small perturbation of the initial condition but also the inter-well jump among attractors around different wells that is led by dramatic disturbance of the initial condition can make the structure harvest vibration energy more conveniently.

4. Global Bifurcation

In this section, we discuss the necessary conditions for heteroclinic and homoclinic bifurcation and the complex dynamics that these global bifurcations induce. To begin with, the heteroclinic orbits and homoclinic orbits are written as [20]:

$$x_{hetero}=\pm \frac{\sqrt{2}{x}_s\sinh (\frac{\chi T}{2})}{{(\kappa +\cosh (\chi T))}^{\frac{1}{2}}},y_{hetero}=\pm \frac{\sqrt{2}{x}_s\chi (1+\kappa )\cosh (\frac{\chi T}{2})}{2{(\kappa +\cosh (\chi T))}^{\frac{3}{2}}}.$$

(46)

and

$$x_{\mathit{homo}}=\pm \frac{\sqrt{2}{x}_s\cosh (\frac{\chi T}{2})}{{(-\kappa +\cosh (\chi T))}^{\frac{1}{2}}},y_{\mathit{homo}}=\mp \frac{\sqrt{2}{x}_s\chi (\kappa +1)\sinh (\frac{\chi T}{2})}{2{(-\kappa +\cosh (\chi T))}^{\frac{3}{2}}}.$$

(47)

where

$$\chi =x_s{}^2\sqrt{2k_2(\frac{x_c{}^2}{x_s{}^2}-1)},\hspace{1em}\kappa =\frac{3x_c{}^2-5x_s{}^2}{3x_c{}^2-x_s{}^2}>0.$$

(48)

The expression of the component $v_{hetero}$ and $v_{\mathit{homo}}$ of the heteroclinic orbits and homoclinic ones obtained by solving Equation (5) is presented as follows

$$v_{hetero}(T)=-\beta e^{-\gamma T}{\displaystyle {\int }_0^T{e}^{\gamma \tilde{T}}{y}_{hetero}(\tilde{T})}d\tilde{T},\hspace{1em}{v}_{\mathit{homo}}(T)=-\beta e^{-\gamma T}{\displaystyle {\int }_0^T{e}^{\gamma \tilde{T}}{y}_{\mathit{homo}}(\tilde{T})}d\tilde{T}.$$

(49)

The dimensionless Equation (5) is approximated to

$$\begin{array}{l}\dot{x}=y,\\ \dot{y}=-x+k_1{x}^3-k_2{x}^5+\epsilon (-\tilde{\mu }y+\widehat{f}\cos \omega T+\tilde{\eta }v),\\ \dot{v}=-\gamma v-\beta y.\end{array}$$

(50)

Since the system (50) is a time-periodic perturbation of a Hamiltonian system, the Melnikov method can be used to describe how the heteroclinic/homoclinic orbits break up in the presence of the perturbation [28]. According to the Melnikov method, the Melnikov function is a signed measure of the distance between the stable and unstable manifolds for the perturbed system [29,30]. If it has simple zeros, there will be intersection of heteroclinic orbits or homoclinic orbits, corresponding to heteroclinic bifurcation or homoclinic bifurcation, respectively. Submitting Equations (46)–(49) to its Melnikov functions under the two types of orbits yield

$$\begin{array}{l}{M}_{hetero}^{\pm }(T_0)=\epsilon {\displaystyle {\int }_{-\infty }^{+\infty }{y}_{hetero}(T)}(-\tilde{\mu }{y}_{hetero}(T)+\widehat{f}\cos (\omega (T-T_0))+\tilde{\eta }{v}_{hetero}(T))dT\\ =-\frac{\mu x_s{}^2\chi }{4(1-\kappa )\sqrt{1-{\kappa }^2}}{I}_{1hetero}-\frac{\eta \beta x_s{}^2{\chi }^2{(1+\kappa )}^2}{2}{I}_{2hetero}+\frac{2\pi \omega }{\chi }{f}_0{x}_s\mathrm{csch}\frac{\pi \omega }{\chi }\cos (\omega T_0)\end{array}$$

(51)

and

$$\begin{array}{l}{M}_{\mathit{homo}}^{\pm }(T_0)=\epsilon {\displaystyle {\int }_{-\infty }^{+\infty }{y}_{\mathit{homo}}(T)}(-\tilde{\mu }{y}_{\mathit{homo}}(T)+\widehat{f}\cos (\omega (T-T_0))+\tilde{\eta }{v}_{\mathit{homo}}(T))dT\\ =-\frac{\mu x_s{}^2\chi }{4(1-\kappa )\sqrt{1-{\kappa }^2}}{I}_{1\mathit{homo}}-\frac{\eta \beta x_s{}^2{\chi }^2{(1+\kappa )}^2}{2}{I}_{2\mathit{homo}}-2f_0{x}_s\sin (\frac{2\omega }{\chi })\sin (\omega T_0)\end{array}$$

(52)

where

$$\begin{array}{l}{I}_{1hetero}=\left(2-\kappa \right)\sqrt{1-{\kappa }^2}+\left(2-4\kappa \right)\arctan \frac{1-\kappa }{\sqrt{1-{\kappa }^2}}, I_{2hetero}={{\displaystyle \int }}_{-\infty }^{-\infty }\frac{\cos \mathrm{h}(\frac{\chi T}{2})}{e^{\gamma T}\left(\kappa +\cos \mathrm{h}(\chi T\right){)}^{\frac{3}{2}}}({{\displaystyle \int }}_0^T\frac{e^{\gamma \tilde{T}}\cos \mathrm{h}(\frac{\chi \tilde{T}}{2})}{\left(\kappa +\cos \mathrm{h}(\chi \tilde{T}\right){)}^{\frac{3}{2}}}d\tilde{T})dT,\\ I_{1\mathit{homo}}=\left(2-\kappa \right)\sqrt{1-{\kappa }^2}-\left(2-4\kappa \right)\arctan \frac{\kappa +1}{\sqrt{1-{\kappa }^2}}, I_{2\mathit{homo}}={{\displaystyle \int }}_{-\infty }^{-\infty }\frac{\sin \mathrm{h}(\frac{\chi T}{2})}{e^{\gamma T}\left(-\kappa +\cos \mathrm{h}(\chi T\right){)}^{\frac{3}{2}}}({{\displaystyle \int }}_0^T\frac{e^{\gamma \tilde{T}}\sin \mathrm{h}(\frac{\chi \tilde{T}}{2})}{\left(-\kappa +\cos \mathrm{h}(\chi \tilde{T}\right){)}^{\frac{3}{2}}}d\tilde{T})dT.\end{array}$$

(53)

In the above equation, the integrals $I_{2hetero}$ and $I_{2\mathit{homo}}$ can be evaluated numerically [20]; $I_{1hetero}$, $I_{1\mathit{homo}}$, $I_{2hetero}$ and $I_{2\mathit{homo}}$ are all positive.

For

$$f_0>f_0^{hetero}=\frac{\mu x_s{\chi }^2\sinh (\frac{\pi \omega }{\chi })}{8\pi \omega (1-\kappa )\sqrt{1-{\kappa }^2}}{I}_{1hetero}+\frac{\eta \beta x_s{\chi }^3{(1+\kappa )}^2\sinh (\frac{\pi \omega }{\chi })}{4\pi \omega }{I}_{2hetero},$$

(54)

there will be a real number of $T_0$ satisfying $M_{hetero}^{\pm }(T_0)=0$ and $M_{hetero}^{\pm }{}^{\prime }(T_0)\ne 0$, meaning that the equilibria of Equation $M_{hetero}^{\pm }(T_0)=0$ are simple, enabling the existence of the transverse heteroclinic orbits. Accordingly, $f_0^{hetero}$ is the amplitude threshold of heteroclinic bifurcation. Similarly, it follows from the Melnikov function (52) that if

$$f_0>f_0^{\mathit{homo}}=\frac{\mu x_s\chi }{8(1-\kappa )\sqrt{1-{\kappa }^2}}\csc (\frac{2\omega }{\chi })I_{1\mathit{homo}}+\frac{\eta \beta x_s{\chi }^2{(1+\kappa )}^2}{4}\csc (\frac{2\omega }{\chi })I_{2\mathit{homo}},$$

(55)

the equation $M_{homo}^{\pm }(T_0)=0$ will have simple equilibria, corresponding to an intersection of the stable and unstable manifolds and the occurrence of homoclinic bifurcation. Given ω = 1.5, critical values can be calculated that $f_0^{\mathit{homo}}\approx 0.089$ and $f_0^{hetero}\approx 0.27$.

A question is then aroused: what dynamical behaviors will be induced by the global bifurcations? To get the answer, we continue to depict the evolution of attractors and their basins of attraction with the excitation amplitude $f_0$ for ω = 1.5 (see Figure 8). For $f_0$ = 0.089 (see Figure 8(a1)), the same as the case of $f_0$ = 0.067 (see Figure 7(d1)), there are still five periodic attractors. To be different, the two higher-amplitude intra-well attractors, i.e., the periodic attractors in yellow and blue, are expanded so that their phases nearly touch the homoclinic obits (see the light grey dashing curves in Figure 8(a1)); meanwhile, it is hard to observe their basins of attraction in Figure 8(a2), as they are eroded seriously. It indicates that for $f_0$ = 0.089, i.e., the theoretical threshold for homoclinic bifurcation, the two higher-amplitude intra-well attractors around two nontrivial well centers become rare attractors. When $f_0$ increase a bit to 0.090, it can be observed in Figure 8(b1) that these two attractors disappear; three lower-amplitude intra-well attractors are left; the intra-well attractor around the origin becomes dominant, as the most area is its basin of attraction (see the green region in Figure 8(b2)). It follows from the comparison of Figure 8(a1,b1) that the homoclinic bifurcation leads to the disappearance of two higher-amplitude intra-well attractors around nontrivial well centers.

When $f_0$ continues to increase, the intra-well attractor around O(0, 0) will be enlarged, and its basin of attraction will erode the basins of attraction of the other two attractors seriously (see Figure 8(c2)). For $f_0$= 0.11, two intra-well attractors around the nontrivial equilibria disappear (see Figure 8(d1)), and the whole initial-condition plane of Figure 8(d2) is green, showing the global stability of the intra-well attractor around O(0, 0). There is no jump now, unfavorable for energy harvesting. As $f_0$ continues to increase, the situation still maintains. For $f_0$ = 0.27, i.e., the theoretical threshold for heteroclinic bifurcation, the phase orbit of the only attractor is enlarged and nearly touches the heteroclinic obits (see the light grey dashing curves in Figure 8(e1)); still the intra-well attractor is globally attractive, as shown in the green map in Figure 8(e2). Then, for $f_0$ increasing a bit to 0.30, a new periodic attractor suddenly appears (see Figure 8(f1)). It is an inter-well periodic attractor whose amplitude is much higher than the intra-well one. It can be concluded that the heteroclinic bifurcation triggers the occurrence of inter-well response. Due to the large amplitude of the inter-well attractor in Figure 8(f1), we have to expand the phase map plane to display it. Note that we still keep the initial plane −1.0 ≤ x(0) ≤ 1.0, −1.0 ≤ y(0) ≤ 1.0 which is totally green, meaning that initial conditions in this region cannot lead to this attractor (see Figure 8(f2)). As $f_0$ grows to 0.31, one can observe its basin of attraction in the initial plane (see Figure 8(g2)). Of course, the occurrence of inter-well oscillation is advantageous for energy harvesting.

Based on Equations (54) and (55), i.e., the theoretical threshold of the dimensionless amplitude $f_0$ for homoclinic bifurcation and heteroclinic bifurcation, we present Figure 9 to show the variation of $f_0$ with the increase of excitation frequency ω. The blue curve and the red one represent the thresholds for homoclinic bifurcation and heteroclinic one, respectively. The numerical values of $f_0^{\mathit{homo}}$ and $f_0^{hetero}$ are obtained at the points for the disappearance of two higher-amplitude intra-well attractors around nontrivial well centers and the occurrence of the inter-well attractor, respectively. For the purpose of energy harvesting, the former is unwanted for reducing jump, while the latter is desirable for inducing high-amplitude inter-well oscillation. Each numerical critical value of $f_0^{\mathit{homo}}$ and $f_0^{hetero}$ is kept to two decimal places. We make sure that if $f_0$ is less than the numerical results $f_0^{\mathit{homo}}$ or $f_0^{hetero}$, the mentioned intra-well attractors will not disappear, or the inter-well periodic attractor will not appear, respectively. In Figure 9, the numerical results for the critical values of $f_0$ are in agreement with analytical ones, verifying the validity of the analysis. Figure 9 also demonstrates that the thresholds of the amplitude $f_0$ for homoclinic bifurcation and heteroclinic bifurcation will increase monotonically with ω. Comparatively, $f_0^{hetero}$ increases more rapidly than $f_0^{\mathit{homo}}$. In addition, Figure 9 shows that the lower the excitation frequency is, the lower threshold of the excitation amplitude to trigger inter-well oscillation will be, which will be helpful to harvest vibration energy.

As $f_0$ continues to increase, the variation of the dynamical behaviors and their basins of attraction can be observed in Figure 10. For $f_0$ = 0.37, the basins of attraction of the inter-well attractor are enlarged and intermingled with the intra-well one, as can be observed from the grey region and grey dots scattered in the green region of Figure 10(a2). That means that the occurrence probability of inter-well oscillation and inter-well jump between the two attractors increases, which is surely good for energy harvesting. When $f_0$ increases to 0.45, the fractality of basins of attraction becomes more obvious, and the basins of attraction of the intra-well attractor are seriously eroded by the basin of the inter-well attractor (see the intermingled green and brown regions in Figure 10(b2)). For $f_0$ = 0.52, there is very little basin of attraction of the intra-well attractor left, showing that the intra-well attractor becomes a rare attractor. As $f_0$ reaches 0.60, the whole initial-condition plane is brown, illustrating that the inter-well response becomes globally attractive, and favorable for the structure to harvest energy.

It is known that the Poincaré map is one of the most useful methods of investigating continuous-time nonlinear systems that involves a discretization technique [30]. For f₀ continuing to increase, we present the bifurcation of the system (5) in the Poincaré map in Figure 11. Here, the points on the Poincaré map are collected from a cross-section at y = 0 and x > 0 in a sufficiently large time interval $5000\le T\le 8000$ to ensure that the response is stable. As can be observed, when $f_0$ is more than 3.0, the system (5) evolves large-range inter-well chaotic motion through the routine of the period-doubling bifurcation. The chaotic vibration for $f_0$ = 4.2 is displayed via its phase map, Poincaré map and spectrum diagram in Figure 12. Note that points of the Poincaré map in Figure 12b are obtained from the cross-section at $T=2n\pi$ satisfying $5000\le T\le 8000$; here n represents integers. It follows from Figure 11 and Figure 12 that global bifurcations finally lead to this large-range inter-well chaotic motion, which is unquestionably useful for energy harvesting.

5. Conclusions

It is well known that energy harvesters favor the appearance of complex dynamics. In this paper, bifurcations and complex dynamics of the electromechanical system of a piezoelectric energy harvester comprising a cantilever piezoelectric beam with an attached tip mass exposed to a harmonic tip force are investigated analytically and numerically.

First, the dimensionless ordinary differential equations governing the electromechanical system are obtained. The equilibrium bifurcation and stability have been investigated for the unperturbed system. The bifurcation sets of the equilibrium in parameter space have been constructed to demonstrate that the number and shapes of potential well depend on the coefficients of the polynomial form of the magnetic force: the smaller the penta power stiffness coefficient of nonlinear magnetic force is, the higher possibility for the occurrence of a triple potential well will be.

Then, fixing the physical parameters of the structure and varying the excitation amplitude and frequency, the case of a triple potential well is discussed. The Method of Multiple Scales is applied to provide analytical intra-well solutions. The oscillator is found to exhibit saddle-node bifurcation leading to bistable intra-well attractors around the same well center. Thus, the system may undergo four intra-well attractors or five intra-well attractors. The numerical integration method was then utilized to verify the solutions.

Consequently, we perform a detailed investigation of basins of attraction of multiple attractors that had still relatively little consideration in the literature. The basins of attraction are obtained by the point-mapping method. The numerical results further confirm the coexistence of these attractors, in good agreement with the theoretical ones, showing the validity of the analysis. It is found that in the regions with evidence of multistability basins of attraction with fractal structures occur quite frequently. We have also found extremely intermingled basins and rare attractors. It also reveals that for the excitation amplitude and frequency varying in the vicinity of saddle-node bifurcation points, the intra-well jump between the two intra-well attractors around the same well center is much easier to trigger than the inter-well jump among intra-well attractors around different well centers. The latter can be induced by a dramatic change of initial position or velocity. Both the two jumps are beneficial for energy harvesting.

Furthermore, analytical expressions for homoclinic orbits and heteroclinic orbits of the unperturbed system are derived. The Melnikov method is successfully employed to detect the analytical criteria for homoclinic bifurcation and heteroclinic bifurcation successively. These results obtained are verified by the numerical simulations in the form of phase maps, basins of attraction, bifurcation diagrams and Poincaré maps where fine agreement is achieved. It is found that the increase of the excitation amplitude can induce homoclinic bifurcation and heteroclinic bifurcation successively; the thresholds of the excitation amplitude for homoclinic bifurcation and heteroclinic bifurcation both increase monopoly with the increase of the excitation frequency. Homoclinic bifurcation induces the disappearance of the intra-well attractors around the nontrivial equilibria. An increase in the excitation amplitude can break the basins into discrete pieces or points. In contrast, heteroclinic bifurcation leads to the occurrence of a large-range jump between inter-well attractor and intra-well one, the globally attractive large-amplitude inter-well attractor, and eventually inter-well chaos through the routine of the period-doubling bifurcation. It implies that for the purpose of energy harvesting, homoclinic bifurcation is undesirable, while heteroclinic bifurcation is favorable.

In this paper, the effect of the excitation amplitude, the excitation frequency as well as initial conditions on the nonlinear behavior of this triple-well piezoelectric energy harvester has been studied. Based on these investigations, the next steps will consist of designing an actual energy harvesting system featuring the calculated position of the equilibria, assessing the performance and optimizing the device. Hence, nonlinear magnetic interactions through the use of two magnets on a base interacting with a magnet at the free end of a cantilever beam can be considered. From an applicative point of view, a comparative analysis with moveable magnets on the base can be considered (see the structure in Figure 13).

It should be pointed out that the results in this paper are limited to the fixed physical properties of the structure. We have not considered the effect of physical properties on inducing the large-scale complex responses which can be beneficial for further optimization of triple-well energy harvesters. Jump phenomena, inter-well oscillation and large-range chaos that we have discussed in detail can make the structure under frequent strain. Considering the material characteristics of the piezoelectric generators, there may be brittles in piezo layers, leading to a challenge for the working lifespan of the devices. For better energy harvesting performance and reliability, these should be taken into account in theoretical study as well as experiments in practical application, which will be included in our future work.