Parametric Analysis of Nonlinear Bi-Stable Piezoelectric Energy Harvester Based on Multi-Scale Method

Piezoelectric energy harvesting technology has advanced quickly in recent years due to the widespread application of wireless sensor networks in industries like architectural structures and health surveillance. Because of its limited operating frequency range, the early linear piezoelectric energy harvesting system [1], [2] was unable to efficiently capture the vibration energy present in the environment.

A bi-stable piezoelectric energy harvester was developed by Erturk et al. [3], [4] using a ferromagnetic cantilever beam and two permanent magnets. The nonlinear mechanical properties of the system were investigated through numerical simulation and experiments. Stanton et al. [5], [6] adopted harmonic balancing and numerical simulation to analyze how different system factors affect the performance of a bi-stable piezoelectric energy harvester. Lan and Qin [7] improved the bi-stable piezoelectric cantilever beam to reduce the potential barrier and shallow the potential well. They also the performance of the energy harvester under random excitation.

Liu et al. [8], [9] attempted to optimize the parameters of the bi-stable piezoelectric energy harvesting system. Panyam et al. [10] measured the influence of three key design parameters on the effective bandwidth of the bi-stable piezoelectric energy harvesting system. Yang and Towfighian [11] presented a new bi-stable piezoelectric energy harvester capable of generating internal resonance, established the mechanical model of the system, and derived the analytical solution. Kim et al. [12] designed a magneto-piezo-elastic vibration energy harvester with reversible hysteresis, and relied on different methods to analyze its resonant response, according to the strength of external excitation.

Firoozy et al. [13] established a complete model for the distributed parameters of the bi-stable piezoelectric cantilever beam, and analyzed the dynamics of the beam numerically. Zhang et al. [14] studied how the hysteresis properties of materials and external excitation affect bi-stable piezoelectric energy harvesters. Zhao et al. [15] designed a new bi-stable piezoelectric energy harvester, and demonstrated that the output performance of the system can be improved by increasing the weight mass block and decreasing the cantilever beam stiffness. Through simulations, Song et al. [16] found that the energy harvesting rate in the bi-stable energy harvester system can be greatly improved by adjusting the position of the beam under external excitation.

Singh et al. [17] accurately modeled a bi-stable sensor based on a nonlinear vibration energy harvester, and proved the validity of the model. Through finite-element analysis, Kumar et al. [18] found that the bi-stable variable-width piezoelectric energy harvester is better than the uniform-width piezoelectric energy harvester, and carried out harvester optimization accordingly. Yao et al. [19] improved the L-shaped piezoelectric cantilever beam, and experimentally verified the bi-stability of the improved beam. Yao et al. [20] created a novel bi-stable piezoelectric-electromagnetic energy harvester, designed mono-stable and bi-stable piezoelectric cantilever beams, and analyzed the beam performance through experiments.

Wang et al. [21] proposed a bi-stable piezoelectric energy harvester with an elastic magnifier. Experimental results show that their harvester is more efficient in energy harvesting than the conventional bi-stable piezoelectric energy harvester. Zhou et al. [22] proposed a new quad-stable energy harvester with higher energy harvesting efficiency than the bi-stable energy harvester. Zhou et al. [23] designed a new quadratic energy harvester that can change the position of the magnets and the distance between the magnets. This new harvester overcomes the low energy conversion rate, a common problem of bi-stable energy harvesters under random excitation. Udani et al. [24] investigated the effect of laminates with different design parameters on bi-stable energy harvesters.

This paper develops a distributed parameter model of the nonlinear bi-stable piezoelectric cantilever energy harvesting system based on the Euler Bernoulli beam hypothesis and geometric nonlinearity of a piezoelectric beam, using the nonlinear magnetic force model and Hamilton variational equation. Galerkin modal decomposition and the multi-scale technique were used to produce analytical representations of the system dynamic response. The relevant data were then obtained by numerically simulating the system performance. This study focuses on how the nonlinear bi-stable cantilever piezoelectric energy harvester system performs as a function of the distance between magnets, damping ratio, external excited amplitude, and external load resistance.

As shown in Figure 1, the nonlinear bi-stable piezoelectric cantilever energy harvester with magnetic coupling consists of a piezoelectric cantilever beam and two magnets (denoted as A and B). The piezoelectric cantilever beam, which is fixed at the base, is composed of a metal base and a pair of piezoelectric layers (PZTs). Two identical PZTs are tightly bonded on the upper and lower surfaces of the base. The two PZTs have opposite polarities in the thickness direction and are electrically connected in series with an equivalent load resistance (R), which represents a low-power electronic device. Magnet A is attached at the tip of the cantilever beam, while magnet B is fixed at the frame. The horizontal gap between the two magnets is denoted by d. In addition, l and b denote the length and width of the beam, respectively; h_s and t_p denote the thickness of the metal base and the PZTs, respectively; e denotes the eccentricity of the tip magnet.

Let $v_b(t)$ be the vibration displacement of the base; $v(s, t)$ be the displacement of the beam at s relative to its fixed end. Then, the nolinear constitutive equations of the piezoelectric beam can be established as:

where, T is stress; S is strain; Y is Young's modulus; $D_3$ is electric displacement; E is electric field strength; $d_{31}$ is piezoelectric constant; $\varepsilon_{33}^T$ is the dielectric constant of piezoelectric materials when constant stress is unchanged; the superscripts s and p represent the parameters related to the base layer and the piezoelectric layer, respectively; the subscripts 1 and 3 represent the $x$ and $y$ directions, respectively; $E_3=-V(t) / 2 t_p$, with $V(t)$ being the voltage. Considering the geometric nonlinearity of the beam [25], the nonlinearity between stress and strain can be expressed as:

where, T is stress; S is strain; Y is Young's modulus; D₃ is electric displacement; E is electric field strength; d₃₁ is piezoelectric constant; $\varepsilon_{33}^T$ is the dielectric constant of piezoelectric materials when constant stress is unchanged; the superscripts s and p represent the parameters related to the base layer and the piezoelectric layer, respectively; the subscripts 1 and 3 represent the x and y directions, respectively; $E_3=-V(t) / 2 t_p$, with V(t) being the voltage. Considering the geometric nonlinearity of the beam [25], the nonlinearity between stress and strain can be expressed as:

The Lagrange function of the system can be expressed as:

where, T_k is the kinetic energy; U_e is the potential energy; W_e is the electrical energy; U_m is the magnetic potential energy. The kinetic energy can be solved by:

where, $m=2 \rho_p t_p b+\rho_s h_s b$, with ρ_p and ρ_s being the density of the piezoelectric layer and the base layer, respectively; M_t and J are mass and the rotary inertia of the tip mass, respectively. The potential energy can be solved by:

where, $h=\frac{h_s}{2} ; \mathrm{YI}=\frac{2}{3}\left[Y_s b h^3+Y_p b\left(3 h^2 t_p+3 h t_p^2+t_p^3\right)\right]$.

Under a constant stress, the strain dielectric constant satisfies $\mathcal{E}_{33}^{{S}_{33}^{T_{31}^2{_p}}}$. The electrical energy can be solved by:

Using the Galerkin method, the displacement $v(s, t)$ can be expressed as a linear combination of the first n-order modal vibration functions:

where, $\varphi_r(s)$ and $\eta_r(t)$ are the r-order mode shape function and the generalized mode coordinates of the cantilever beam, respectively. Then, we have:

where, $\delta_{r s}$ is the Kronecker delta; $\omega_r=\lambda_r^2 \sqrt{Y I /\left(m l^4\right)}$ is the resonance frequency of the r-th mode; $\lambda_r$ is the eigenvalue. The values $\varphi_r(s)$ and $\lambda_r$ can be calculated by the methods of Stanton et al. [6].

The magnetic potential energy can be solved by:

where, $k_0=\frac{2 \kappa}{d}$; $k_1=\frac{\kappa\left(12 \varphi_1(l)^2+2 d^2 \varphi_1^{\prime}(l)^2+6 d \varphi_1(l) \varphi_1^{\prime}(l)\right)}{d^5}$; $k_2=\frac{12 \kappa\left(3.75 \varphi_1(l)^4+0.25 d^4 \varphi_1^{\prime}(l)^4+d^2 \varphi_1(l)^2 \varphi_1^{\prime}(l)^2+0.5 d^3 \varphi_1(l) \varphi_1^{\prime}(l)^3+2.5 d \varphi_1(l)^3 \varphi_1^{\prime}(l)\right)}{d^7}$; $\kappa=\frac{\mu_0 M_A V_A M_B V_B}{4 \pi}$.

Only considering the first order mode, Eq. (3) can be substituted with the following Lagrange variational equation under the normalization condition of mode function:

where, $F_1(t)=-2 \xi_1 \omega_1 \dot{\eta}_1(t)$ is the generalized dissipation force; $\xi_1$ is the damping ratio; $\dot{Q}(t)=-V(t) / R$ is the generalized output charge. Based on Eq. (11), the electro-mechanical coupling equation of the piezoelectric energy harvesting system can be derived as:

where, $\quad p_1=2 Y I \int_0^l \varphi_1^{\prime 2}(s) \varphi_1^{\prime \prime 2}(s) d s \quad ; \quad \theta_1=Y_p b d_{31}\left(h+\frac{t_p}{2}\right) \int_0^l \varphi_1^{\prime \prime}(s) d s \quad ; \quad \theta_2=\frac{3}{2} Y_p b d_{31}(h+$ $\left.\frac{t_p}{2}\right) \int_0^l \varphi_1^{\prime 2}(s) \varphi_1^{\prime \prime}(s) d s$ ; $\Gamma_1=m \int_0^l \varphi_1(s) d s+M_t\left(\varphi_1(l)+e \varphi_1^{\prime}(l)\right)$ ; $C_p=\frac{b l \varepsilon_{33}^S}{2 t_p}$.

Next, the base incentive is configured as $\ddot{v}_b(t)=\bar{v}_b \cos \left(\omega_e t\right)$, where $\bar{v}_b$ is the amplitude of acceleration excitation, and $\omega_e$ is the circular frequency of excitation. Through a dimensionless transform: $x=\eta_1 / l, \bar{V}=V / e$, $e=l \theta_1 / C_p$, and $\tau=\omega_1 t$, we have:

The potential energy of the target system totals:

Figure 2 shows the potential energy curves of the system with different distances between the two magnets. When $1-K_1>0$ and $K_2=0$, the potential energy curve has only one potential well, making the system linear and monostable. When $1-K_1>0$ and $K_2 \neq 0$, the system belongs to the nonlinear monostable state. When $1-K_1<0$ and $K_2>0$, the potential energy curve has a potential barrier (unstable equilibrium point, coordinate $x_s$) and two potential wells (stable equilibrium point, coordinate $x \sqrt{\left(K_1-1\right) / K_2}$), making the system bi-state. The following analysis focuses on the intra-well and inter-well motions of bi-stable systems.

By introducing a new variable $x s_t$ into Eqs. (12) and (13), the motion equations at the point of stable equilibrium can be expanded as:

where, $\widetilde{\omega}_1=\sqrt{2\left(K_1-1\right)}$ is the natural frequency of linearization of the single potential well; $\delta=$ $3 \sqrt{\left(K_1-1\right) K_2}$ is the quadratic nonlinear term coefficient; $x_t$ is the motion trajectory around the nontrivial equilibrium point. The remaining coefficients can be described in turn as $q_1=\theta_2\left[\left(K_1-1\right) / K_2\right], q_2=$ $2 \Theta_2 \sqrt{\left(K_1-1\right) / K_2}, q_3=\beta\left[\left(K_1-1\right) / K_2\right]$, and $q_4=2 \beta \sqrt{\left(K_1-1\right) / K_2}$.

By introducing a small perturbation parameter $\mathcal{E}$, the new independent time variable T_n can be expressed as:

Taking the derivative with respect to $\tau$:

The displacement and output voltage response of the system can be respectively expressed as:

Next, the constant parameters in Eq. (15), namely, nonlinear term coefficients, electro-mechanical coupling coefficients, and the excitation force, are scaled to unify the order of the viscous damping effect in the perturbation problem. That is, we let

If the excitation frequency of the intra-well motion is very close to the natural frequency of a single well, then

where, $\sigma$ is the tuning parameter of excitation frequency.

Substituting Eqns. (19), (20), (21) and (22) into Eqns. (17) and (18), omitting the terms above $\varepsilon^2$, and setting the coefficients of $\varepsilon^0, \varepsilon^1$ and $\varepsilon^2$ as zero, we have:

The solutions of Eqns. (22) and (23) are as follows:

where, $A\left(T_1, T_2\right)$ is a complex-valued function; cc is the complex conjugate of the preceding term. Substituting Eqns. (28) and 29) into Eqns. (24) and (25) and eliminating the secular term, we have:

On this basis, the solutions of Eqns. (33) and (34) can be expressed as:

Substituting Eqns. (28), (29) and (30) into Eq. (26) and eliminating the secular term, we have:

Substituting $A=\frac{1}{2} a e^{i \theta}$ and $\gamma=\sigma T_2-\theta$ into Eq. (33) and separating the real part and imaginary part, we have:

where, $c_1=\widetilde{\xi}_1+\frac{\left(\widetilde{\Theta}_1+\tilde{q}_1\right)\left(1+q_3\right) \alpha}{2\left(\widetilde{\omega}_1^2+\alpha^2\right)} ; c_2=\frac{\beta\left(\widetilde{\Theta}_1+\tilde{q}_1\right) \alpha}{8\left(\widetilde{\omega}_1^2+\alpha^2\right)}+\frac{\widetilde{\Theta}_2\left(1+q_3\right) \alpha}{8\left(\widetilde{\omega}_1^2+\alpha^2\right)}+\frac{\tilde{q}_2 q_4 \alpha}{8\left(4 \widetilde{\omega}_1^2+\alpha^2\right)} ; c_3=\frac{\widetilde{\Theta}_2 \beta \alpha}{32\left(\widetilde{\omega}_1^2+\alpha^2\right)}+\frac{\widetilde{\Theta}_2 \beta \alpha}{32\left(9 \widetilde{\omega}_1^2+\alpha^2\right)} ; c_4=\frac{\tilde{f}}{2 \widetilde{\omega}_1} ;$$c_5=\sigma-\frac{\left(\widetilde{\Theta}_1+\tilde{q}_1\right)\left(1+q_3\right) \widetilde{\omega}_1}{2\left(\widetilde{\omega}_1^2+\alpha^2\right)} ; c_6=\frac{\left(\widetilde{\Theta}_1+\tilde{q}_1\right) \beta \widetilde{\omega}_1}{8\left(\widetilde{\omega}_1^2+\alpha^2\right)}-\frac{5}{12 \widetilde{\omega}_1}\left(\frac{\tilde{\delta}}{\widetilde{\omega}_1}\right)^2+\frac{3}{8 \widetilde{\omega}_1} \widetilde{K}_2+\frac{\widetilde{\Theta}_2\left(1+q_3\right) \widetilde{\omega}_1}{8\left(\widetilde{\omega}_1^2+\alpha^2\right)}+\frac{\tilde{q}_2 q_4 \widetilde{\omega}_1}{4\left(4 \widetilde{\omega}_1^2+\alpha^2\right)} ; c_7=\frac{\widetilde{\Theta}_2 \beta \widetilde{\omega}_1}{32\left(\widetilde{\omega}_1^2+\alpha^2\right)}+\frac{3 \widetilde{\Theta}_2 \beta \widetilde{\omega}_1}{32\left(9 \widetilde{\omega}_1^2+\alpha^2\right)}$.

For the steady-state response, setting $D_2 a$ and $a D_2 \gamma$ to zero, squaring Eqs. (36) and (37) and adding up the resulting equations, we have:

where, a is the amplitude of steady-state displacement response, which can be obtained by Eq. (36). The steady-state solutions for the displacement and the voltage can then be obtained by:

where, $\phi_0=\arctan \left(\frac{c_1+c_2 a^2+c_3 a^4}{c_5-c_6 a^2-c_7 a^4}\right) ; \varphi=\arctan \left(\frac{\alpha}{\widetilde{\omega}_1}\right)$.

The steady-state output power amplitude of the intra-well motion can be written as:

For the inter-well motion, the local stiffness near the saddle point (the unstable equilibrium point) is negative $\left(1-K_1<0\right)$. In this case, the conventional multi-scale method is not applicable. Therefore, the authors set the damping and external excitation terms in Eq. (20) and coefficients $\beta_2, \beta_3$ and $\beta_4$ as $\mathcal{E}-orders$ [26]. Then, we have:

It is assumed that $x=\bar{A} \cos (\widetilde{\omega} \tau)$ is the approximate solution of Eq. (40). Substituting this term into Eq. (40), we have:

Adding and subtracting the $\omega^2 x$ term to the left side of Eq. (14), the following can be derived from Eq. (43) [26]:

where, $\varepsilon \bar{\xi}_1=\xi_1 ; \varepsilon \widetilde{\omega}^2=\omega^2 ; \varepsilon \vartheta=1-K_1 ; \varepsilon \bar{K}_2=K_2 ; \varepsilon \bar{\Theta}_1=\Theta_1 ; \varepsilon \bar{\Theta}_2=\Theta_2 ; \varepsilon \bar{f}=f$.

According to the derivation process in Section 3.1, the steady-state solutions for the displacement and the voltage of the inter-well motion can then be expressed as

where, $\bar{\phi}_0=\arctan \left(\frac{\bar{c}_1+\bar{c}_2 \bar{a}^2+\bar{c}_3 \bar{a}^4}{\bar{c}_5-\bar{c}_6 \bar{a}^2-\bar{c}_7 \bar{a}^4}\right) ; \bar{\varphi}=\arctan \left(\frac{\alpha}{\omega}\right) ; \bar{a}$ is the amplitude of the steady-state amplitude response for the inter-well motion:

where, $\bar{c}_1=\bar{\xi}_1+\frac{\bar{\theta}_1 \alpha}{2\left(\omega^2+\alpha^2\right)} ; \quad \bar{c}_2=\frac{\bar{\theta}_1 \beta \alpha}{8\left(\omega^2+\alpha^2\right)}+\frac{\bar{\theta}_2 \alpha}{8\left(\omega^2+\alpha^2\right)} ; \quad \bar{c}_3=\frac{\bar{\theta}_2 \beta \alpha}{32\left(\omega^2+\alpha^2\right)}+\frac{\bar{\theta}_2 \beta \alpha}{32\left(9 \omega^2+\alpha^2\right)} ; \bar{c}_4=\frac{\bar{f}^2}{2 \omega} ; \bar{c}_5=\frac{\widetilde{\omega}^2}{2 \omega}-\frac{\vartheta}{2 \omega}-$ $\frac{\bar{\theta}_1 \omega}{2\left(\omega^2+\alpha^2\right)} ; \bar{c}_6=\frac{\bar{\theta}_1 \beta \omega}{8\left(\omega^2+\alpha^2\right)}+\frac{3 \bar{K}_2}{8 \omega}+\frac{\bar{\theta}_1 \omega}{8\left(\omega^2+\alpha^2\right)} ; \bar{c}_7=\frac{\bar{\theta}_2 \beta \omega}{32\left(\omega^2+\alpha^2\right)}+\frac{3 \bar{\theta}_2 \beta \omega}{32\left(9 \omega^2+\alpha^2\right)}$.

The steady-state output power amplitude of the inter-well motion can be expressed as:

This section discusses the analytical expressions of the steady-state response of the intra-well and inter-well motion of the bi-stable system developed in the preceding part. The analysis and discussion center on how several factors, including load impedance, excitation amplitude, damping ratio, and distance between magnets, affect the system performance. The physical parameters of the system were configured as follows: $l=50 \mathrm{~mm}$, $b=20 \mathrm{~mm}$, $h \mathrm{~mm}_s, t \mathrm{~mm}_p$, $l_a=5 \mathrm{~mm}$, $\rho \mathrm{kgm}_s^3$, $Y 9 \mathrm{Nm}_s^2$, $\rho \mathrm{kgm}_p^3$, $Y 9 N^2 p$, $d_{31}=-3.2 \times 10^{-10} \mathrm{C} / \mathrm{N}$, $\varepsilon_{33}^{s^{-8} F m}$, $M-3 \mathrm{~kg}_t$, $V_A=V_B=1 \times 10^{-6} \mathrm{~m}^3$, $M_A=M_B=1.16 \times 10^6 \mathrm{~A} / \mathrm{m}$, and $\xi_1=0.01$.

The multi-scale method may yield many solutions to the intra-well and inter-well steady-state responses, but not all of them are stable solutions. Therefore, it is necessary to analyze the stability of the solutions. Taking intra-well solution as an example, Eqns. (36) and (37) are linearized at $(a, \gamma)$ to form an autonomous differential equation with respect to the perturbations $\Delta a$ and $\Delta \gamma$ [27]:

Using Eqns. (36) and (37) to eliminate $\gamma$ in the above formula, the feature equation can be obtained as:

Expanding the determinant, we have:

Following the Routh-Hurwitz criterion, the root of the above algebraic equation has a negative real part under the following condition:

According to the Lyapunov theory, the solutions of Eqns. (36) and (37) are stable. The stability of the inter-well motion solution can be determined by the same method.

In the following figures on calculation results, the stable and unstable solutions are represented by solid and dotted lines, respectively.

Figure 3 shows the frequency response curves of the system displacement and output power at the resistance R=200kΩ, f=0.0021 and different distances between magnets (d=[12mm, 14mm, 16mm]). Figure 3 reveals that the inter-well motion's frequency response curve is inclined to the high frequency band, exhibiting hardening features, while the intra-well motion's frequency response curve is inclined to the low frequency band, exhibiting softening features. While the frequency band width of the inter-well motion diminishes as the distance between the magnets increases, the displacement and peak output power of the motion rise. As the distance between the magnets grows, the frequency band of the frequency response curve of the intra-well motion shifts to the lower frequency band. In addition, Figure 3 contrasts the calculation results of the two models, which considers and ignores geometric nonlinearity, respectively. The results demonstrate that: when geometric nonlinearity is taken into account, the inter-well motion's displacement, output power peak, and frequency band width all significantly increase, and the frequency response curve is skewed toward the high frequency band, exhibiting more hardening features. The displacement and peak output power of the intra-well motion do not significantly change, but the frequency response curve is pushed to the high frequency band.

Figure 4 shows the frequency response curves of the system displacement and output power at the resistance R=200kΩ, d=16mm and different excitation force amplitudes (f=0.0017,0.0021,0.0027]). The results show that the displacement, peak output power and frequency band width of the inter-well motion increased greatly with the growing excitation amplitude. The displacement and peak output power of the intra-well motion increased with the excitation amplitude, but the increase was very small.

Figure 5 shows the variation of the peak displacement and peak output power for inter-well motion with excitation force amplitudes, at resistanceR=200kΩ and with different distances between magnets (d=[12mm, 16mm]). It can be seen that There is a critical excitation force amplitude for the inter-well motion; when this critical value is reached, the displacement peak and output power peak of the inter-well motion are both small, and little changes occur as the excitation force amplitude increases; The displacement peak and output power peak of the inter-well motion rapidly increase when the excitation force amplitude exceeds the critical value, and they continue to rise as the excitation force amplitude rises. The critical value of the inter-well motion's excitation force amplitude can be decreased by shortening the distance between magnets. But this also results in decreased peak displacement and peak output power. It can also be seen that, when geometric nonlinearity is taken into account, the critical value of the inter-well motion excitation force amplitude decreases. This is demonstrated by comparing the calculation results of the two models considering geometric nonlinearity and ignoring geometric nonlinearity.

Figure 6 shows the frequency response curves of the system displacement and output power at the resistance R=200kΩ, d=16mm, f=0.0033, and with different damping ratios (ξ₁=[0.01,0.02,0.05]). It can be seen that reducing the damping ratio of the system can effectively broaden the inter-well motion’s frequency band of the system. By contrast, the intra-well motion is not sensitive to the change of damping ratio.

Figure 7 shows the variation of the peak displacement, peak output power, and frequency band width of the inter-well motion at different distances between the magnets (d=[12mm,14mm,16mm]), along with changes in the load impedance when f=0.0021. As shown in Figure 7(a) and Figure 7(c), the peak displacement and frequency band width of the inter-well motion first decreased and then increased with the rising load impedance. The growing rate of peak displacement picks up speed, and that of the band width slows down, with the increase of the magnet spacing. Comparing the calculation results of the two models considering geometric nonlinearity and ignoring geometric nonlinearity, it is clear that the deviation of peak displacement and frequency band width calculated by the two models decreases with the load impedance. As suggested in Figure 7(b), as the load impedance increases, a load impedance appears in the short-circuit and open-circuit states, respectively. As a result, the output power of the inter-well motion reaches the maximum. Hence, increasing the distance between the magnets can significantly increase the optimal load impedance in the open circuit state, but has little effect on the optimal load impedance in the short circuit state.

The geometric nonlinearity of the beam is taken into account to produce the electromechanical coupling equations of the nonlinear bi-stable cantilevered piezoelectric energy captive system based on the Euler-Bernoulli beam assumption. The Galerkin modal decomposition approach and the multi-scale method are also used to reach the analytical solution of the equation. The following conclusions are reached after comparing and analyzing how the piezoelectric material affect system performance in the two situations of considering and ignoring geometric nonlinearity:

(1) As the distance between two magnets increases, the intra-well and inter-well motion's displacement and peak output power both rise, while the inter-well motion's frequency bandwidth narrows. The inter-well motion frequency response curve's hardening effect is evident, and the peak values of displacement and output power of the intra-well and inter-well motions taking geometric nonlinearity into account are obviously larger than those of the model disregarding geometric nonlinearity. The operational frequency band can be expanded and the output power of the inter-well motion can be increased by varying the distance between the two magnets.

(2) The inter-well motion has a critical exciting force amplitude. The displacement peak and output power peak of the inter-well motion grow obviously with the increase in the excitation force amplitude when it exceeds the critical value. As the distance between two magnets grows, so does the amplitude of the critical excitation force.

(3) To maximize the output power of the inter-well motion, there is an ideal load impedance that may be changed by adjusting the distance between two magnets.