Modeling and analysis of magnetically coupled piezoelectric dual beam with an annular potential energy function for broadband vibration energy harvesting

Conventional piezoelectric cantilever-based vibration energy harvesters have narrow bandwidth. In this article, we develop a dual-beam piezoelectric energy harvester featuring an annular potential energy function that can harvest vibration energy over a wide spectrum under small amplitude excitations. The proposed harvester contains two conventional piezoelectric cantilevers placed orthogonal to each other which are coupled by repulsive magnetic force. We demonstrate analytically and numerically that a new annular potential energy function can be built with proper configuration. In the new annular stable state, the harvester can detour around the potential barrier rather than jump over it, yielding large amplitude voltage outputs throughout a wide spectrum. Case studies were carried out, and it is proved that the proposed annular stable harvester has a bandwidth of 3.9 Hz and a voltage output performance 3.01 times better than that of a conventional bistable one under excitations of 3 m/s2. The nonlinear dynamics of the proposed harvester are analyzed in detail.


Introduction
Wireless sensors and portable devices have received significant attention due to the rapid growth of the Internet of things (IoT). They rely on chemical batteries as the source of power. However, the replacement of batteries in some situations is tedious and expensive. Vibration energy harvesting is one of the potential solutions for powering IoT sensors [1]. A typical vibrational energy harvester is a second-order resonator consisting of a cantilever beam, a piezoelectric plate, and a proof mass [2]. On the other hand, when the resonant frequency of the device matches that of the excitation, the system will perform optimally. A slight resonance frequency shift would yield drastically performance decrease [3].
Previously, researchers have proposed various approaches to improve the energy harvesting performance in terms of boarding the operating bandwidth [4], including resonant frequency tuning [3][4][5][6], introducing nonlinearity [7][8][9][10], frequency-up conversion [11,12], multimodal configuration [13], and internal resonance [14][15][16]. Among these, nonlinear dynamics received the most significant attention. As a typical nonlinear piezoelectric energy harvester (PEH), the Duffing-type PEHs have attracted great attention and have been investigated extensively [17]. This is because the voltage output of a piezoelectric energy harvester is dependent on the displacement of the cantilever. Therefore, the large amplitude motion of a nonlinear system can yield significantly enhanced output over a wide spectrum. Theoretical investigations, simulations, and experimental studies were carried out. Stanton et al. [18] found that both softening and hardening responses of a monostable PEH could contribute to a remarkable increase in the bandwidth. In addition, it was revealed that the inter-well oscillations of a bistable PEH could achieve a large amplitude response over a wider frequency range than that of a monostable PEH [19][20][21]. And Litak et al. [22] investigated the features of a bistable PEH under random excitation. In particular, the bistable PEH shows promising advantages for harvesting low-frequency vibration energy [19][20][21][22]. Masana et al. [23] demonstrated the performance of a harvester in mono-and bistable states and concluded that the features of the system are dependent on the shape of the potential functions.
On the other hand, the conventional bistable systems work optimally under large amplitude excitations only. This is because the large amplitude motion in a bistable system is usually induced by the inter-well motion, where large excitations are necessary for a jump motion from one of the stable points to the other. To overcome this limitation, research was carried out on modifying the shape of the potential barriers. For example, Leadenham et al. [24] proposed an M-shaped piezoelectric cantilever for ultra-wideband energy harvesting under small excitations by reducing the barrier of the potential well. For another example, a potential well with multiple stable states would also enable large amplitude motion [25][26][27]. Zhou et al. [27] presented a broadband piezoelectric energy harvester with a triple-well potential. Compared to conventional bistable energy harvester with a deeper potential well, the tristable one passes easily through potential barrier for realizing the inter-well oscillation. Inter-well modulation is used to reduce the barrier height and make inter-well oscillation easier to achieve [28][29][30][31][32][33][34][35]. Wang et al. [32] designed a springbased bistable energy harvester to achieve high working efficiency through inter-well modulation. Besides, magnetic force was utilized to dynamically alter the potential energy threshold in the nonlinear system, which contributes to the system maintaining inter-well motion and enhancing the working performance under lower intensity excitation. Fang et al. [35] proposed a tuned bistable energy harvester to outperform the conventional bistable energy harvester. The tuning effect can dynamically reduce the elastic potential barrier. In addition, it was demonstrated that an external impact force can trigger a jump phenomenon between two potential wells under small amplitude excitations [36].
Typically, the potential barriers in conventional nonlinear energy harvesters are essentially one-dimensional. The conventional bistable PEH has to jump over the potential barriers between the stable points. In this research, we extend the conventional one-dimensional potential well into two-dimensional space. We hypothesize that by extending the potential well into two-dimensional space, the harvester can circumambulate around the potential barrier for a large amplitude motion. In such a case, the PEH does not need to gather a large amount of mechanical energy for the large amplitude motion. And the piezoelectric energy harvester would have significantly enhanced performance under excitations with small amplitudes.
The paper is organized as follows: the concept and system configurations are presented in Sect. 2. The detailed derivation of the analytical solution of the model by the harmonic balance method is presented in Sect. 3. Case studies and discussions of the performances of the system in monostable and annular stable states are presented in Sect. 4. Finally, the conclusions are provided in Sect. 5.

Schematic of design
The schematic and the lumped parameter model of the proposed PEH are presented in Fig. 1a and b. The PEH consists of two identified piezoelectric cantilever beams placed orthogonal to each other: in particular, beam A vibrations in the x-direction and beam B vibrations in the y-direction. Besides, the two cantilevers are coupled by magnetic force. Stainless steel sheets were chosen as the substrate for the two beams. Permanent magnets and two pieces of lead zirconate titanate (PZT) sheets are attached to the tip and root of the two beams, respectively. Here, repulsive force between the two magnets is chosen for generating a two-dimensional interaction between the two cantilevers. Indeed, similar configurations were analyzed previously [28][29][30][31][32][33][34][35]37]. On the other hand, these researches focused on one-dimensional features of adjustable potential wells or multiple directional harvesting. In other words, the harvesters have to jump over the potential barriers for the large amplitude motion.
Notably, the configuration Fig. 1a potentially yields an annular two-dimensional potential well in Fig. 1c with proper parameters combination. We hypothesize that with the annular two-dimensional potential well, the harvester can circumambulate around the potential barrier for a large amplitude motion. Hence, the PEH does not need to gather a large amount of mechanical energy for the large amplitude motion. To this end, this research aims to derive the analytical solution of the magnetically coupled dualbeam PEH system. We demonstrate that, with the proper arrangement of the permanent magnets placed on the tip of the cantilevers, the harvester may have a complicated two-dimensional potential well by adjusting the distance between the magnets. Hence, the voltage output performance of both beams is expected to be enhanced over a wide frequency range under small amplitude excitations. In the following analysis, the dynamic characteristics of the proposed system are studied in monostable, and annular stable states with a comprehensive parametric analysis.

Governing equation
In this section, a mathematical model is presented to analyze the magnetically coupled piezoelectric cantilevers system. The assumed-mode approach is adopted for analytical modeling, and its vibratory motion relative to the base can be represented by a uniformly and convergent series of eigenfunctions [2] where / i;r and g i;r are the mode shape functions and modal displacement expressions of the r-th vibration mode, respectively. We adopt the first beam bending mode / i;1 to characterize the beam moves along the bending direction, that is, w i ¼ / i;1 g i;1 . Besides, the linear constitutive relation of the piezoelectric transducer is [38].
where d 1 , r 1 , D 3 ; and E 3 represent, the strain, stress, electrical displacement (charge/area), and electrical field (voltage/length) of the piezoelectric transducer, respectively, and Y 11 , d 31 , and e 33 are Young's modulus, piezoelectric constant, and dielectric constant of the material. Without loss of generality, here the two cantilevers have identical parameters, i.e., we have the same model shape and parameters of the two cantilevers. Therefore, the lumped masses, lumped stiffnesses, and electromechanical coupling coefficients of the two cantilevers have the same value. In the following analysis, we use a single lumped mass, stiffness, and coupling coefficient for the system dynamics analysis. Specifically, the lumped coefficients related to material constants and structural geometry for beam A can be derived using the assumed modal method and F p represent, respectively, the length of the piezoelectric transducer, the length of the beam, the width of the composite cantilevers, the cross-sectional area of the piezoelectric transducers, the cross-sectional area of the beams, the mass density of the beams, and the mass density of the transducers, the proof mass, Young's modulus of beams, the moments of inertia of the beams, and the first moment of area of the transducers of the two beams. And beam B has the same parameters. The bending stiffness term of the composite cantilever E a I a can be obtained as where h p and h b represent, respectively, the thickness of the transducer and the beam. a ¼ between the connecting layer between the transducer and the beam and the neutral axis with zero strain. The piezoelectric transducer is simplified as a current source and a shunt capacitor [2], as shown in Fig. 1c.
Note that beam B has the same configuration as beam A. Hence, the lumped parameters of beam B have the same values as that of beam A. Besides, x ¼ w 1 l b ; t ð Þ and y ¼ w 2 l b ; t ð Þstand for the displacement of the tips of the two beams. The magnetic coupling dual-beam PEH governing equation can be written as In Eq. (5a) and (5c), c 1 and c 2 are the mechanical damping terms of two beams, respectively, and a x and a y are the amplitude of excitations in the x-and ydirections, respectively. U m x; y ð Þ is the two-dimensional potential energy of the magnetic field. The correction factor of the forcing function l 1 for beam A can be obtained as [2] Beam B has the same parameters as that beam A due to the identity configuration. The main difference between the different magnet arrangements results in the different potential energy of the magnetic field, which can be described by the dipole-dipole model [39] U m x; y ð Þ ¼ where M 1 and M 2 are the moments of the magnetic dipoles of the magnets attached to the two cantilevers, respectively; s is the vacuum permeability; and D is the distance between the magnetic dipoles. For attractive magnets, we have M 1 ¼ ÀM 2 , and for repulsive magnets, we have M 1 ¼ M 2 . In this article, the repulsive magnetic force is chosen for generating complicated potential well. And M 1 ¼ M 2 ¼ 1:146A Á m 2 , s ¼ 4p Â 10 À7 H=m are chosen. Thus, the total potential energy and the restoring forces in different directions are given It can be obtained from Eq. (8) that restoring forces in two-dimensional space is depending on the magnetic field and the stiffness of the beams. Adjusting the distribution of the magnetic field can effectively yield a complicated potential well in the two-dimensional space.
2.3 Influence of the distance on the potential energy shapes and equilibria In this section, the features of the two-dimensional potential barrier are analyzed. The static characteristics of magnetic coupled PEH, the shape of the potential well, the equilibria, and natural frequencies are analyzed analytically. The calculated lumped parameters of the two beams are listed in Table 1.
Notably, resistance loads of 50 MX are chosen for evaluating the open circuit responses of the system. The potential energy of the system is evaluated. Figures 2 and 3 depict the potential energy of the system with different distances D between the two magnets. The figures of the shape of the potential wells are obtained from calculations using Eq. (8a). The potential energy shapes and contour lines show that the PEH is almost a linear one and the magnetic force is negligible when D = 30 mm, as shown in Figs. 2a and 3a. When D decreases to 18 mm, there is a twodimensional valley in the potential well only. On the other hand, the bottom of the potential well becomes flatter than that of the linear one, as shown in Figs. 2b and 3b. In this case, the system has one stable point at the valley of the potential well with nonlinearity introduced. This case may have a similar feature to that of a monostable system in the conventional nonlinear PEH with a one-dimensional potential well. In addition, the height of the potential well at the central point arises when D decreases to 14 mm and 15 mm, as shown in Figs. 2c, d and 3c, d. In these cases, the single-point potential valley changes into an annular one. Besides, it can be obtained that reducing the value of D can effectively boost the height of the peak at the center of the potential well, as shown in Fig. 2c and d. The trend is similar to that of a conventional nonlinear PEH, where proper parameters adjustment of the magnetic field would switch the monostable potential well into a bistable one.
It is also worth mentioning that the bistable PEHs with dynamic potential wells [28][29][30][31][32][33][34][35] may have similar potential functions in two-dimensional space. On the other hand, the cantilevers in the abovementioned research have large or different stiffnesses. In particular, the systems with large stiffness cantilevers are essentially belonging to the case in a monostable state, as shown in Fig. 2b. Besides, the systems with cantilevers of different stiffnesses may have an elliptical potential energy function in the twodimensional space. The system would still have large potential energy variation during the vibrations, i.e., requiring a large amount of external energy for the large amplitude motion. In this research, the proposed system has cantilevers with the same parameters. And a proper parameter combination is chosen for creating an annular stable state. In the new stable state, the PEH does not need to absorb a lot of external energy for circumambulating around the potential barrier. Therefore, this configuration has the potential to operate optimally under small amplitude excitations.
Note that beams A and B have one-dimensional freedom only, the potential curve for each beam is a dynamic one. In other words, one beam has a dynamic potential energy curve depending on the displacement of the other one. The real-time potential well for a single beam would be a one-dimensional cross-section curve of the two-dimensional potential wells. We take an example of the system in an annular stable state when D = 14.5 mm. It can be obtained that changing the displacement of beam B results in different shapes of the potential well, as shown in Fig. 4. It can be obtained that the potential well has a large peak when the displacement Y is 0 mm and 1 mm. Meanwhile, the shape of the potential well becomes a flat one when Y = 5 mm. Consider that the displacement of beam B, the Y, is a dynamic value. Hence, the potential well for beam A is indeed a dynamic one. It is worth emphasizing that beams A and B have the same resonant frequency. The flexible potential well is fundamentally hinged upon the resonant coupling of the two beams and thereby maximizes the efficiency. This phenomenon demonstrates the novel nonlinearity of the system that the beam may vibrate under a potential well with a large and a small barrier. Conventionally, the shape of the potential well is constant in the monostable and bistable PEH system. It would be difficult for a conventional bistable PEH system to jump over the barrier in the potential well. In contrast, beam A may have a large amplitude motion easily under a flatter potential well when beam B has a large displacement. The state of the stability of the system, i.e., the monostable and the annular stable state, is dependent on the distance between the two cantilevers. The types of stable states can be analyzed by calculating the equilibria of the two-dimensional potential well. And the equilibria of the dual-beam PEH are calculated from Eq. (8b) and (8c) where X 0 and Y 0 are the equilibria positions of beam A and beam B, respectively. Considering the symmetric of the system, here we define stiffness Subtracting K into Eq. (9), the relations between these two beams' equilibrium positions are given The differences between the cases in monostable and annular stable states stem from the location of the lowest point of the potential well. In addition, the monostable and annular stable states are dependent on the distance between the two magnets. Firstly, we consider a case in the monostable state. Here, the system has a single-point valley in the potential well. In such a case, the solutions of Eq. (10) should be Therefore, the distance D for the system in monostable can be obtained as D ! ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The critical value for the case in this research is obtained as D 0 = 15.75 mm. It is worth mentioning that under the critical value condition of D 0 , the system remains in the monostable state. Notably, the system has system dynamics approximate a linear one when D has a value large enough. However, there is no critical value for the system switching from the monostable state into the semi-linear one.  . In the annular stable state, the system has features of the circularshaped potential well. And the equilibrium positions can be obtained by solving.
where the solution of Eq. (12a) denotes the stable equilibrium positions and the solution of Eq. (12b) stands for an unstable point of the system, as shown in Fig. 5. Particularly, the positions of the equilibrium point solved by Eq. (12a) stand for the lowest points of potential energy in Figs. 2 and 3, standing for the circular-shaped potential well of the system in the annular stable state. Figure 5 shows the stability point of the system varies with the distance of the magnet. When the magnet distance is large, the system is in a monostable state and the stable point is (X, Y) = (0, 0). As the distance decreases, the system becomes an annular steady state when it is less than the critical value D 0 , and the stability points are distributed on an arc and the unstable point is (X, Y) = (0, 0). Moreover, as the distance decreases further, the radius of stability increases, which means that the amplitude of the inter-well oscillation is also larger. When the distance becomes zero, the radius reaches its maximum value, which is the same as the critical value.

Influence of the distance on the natural frequencies
The strength of the magnetic coupling influences the shape of the two-dimensional potential well. In addition, the natural frequencies of the system are dependent on the magnetic force. Here, we linearize the equation of motion to analyze the effect of magnetic forces on the natural frequencies. The nonlinear restoring forces are expanded using the Taylor series at (X, Y) = (0, 0). Omitting the higher orders, the magnetic forces can be simplified as oU m x; y ð Þ By omitting the cubic term of the magnetic force in Eq. (13), the natural frequencies of the system are given where k i , m i ; and f i represent, respectively, the beam A (i = 1) and beam B (i = 2) of mass, stiffness, and resonant frequency of the first mode. And the influences of magnet distance on the natural frequencies for the monostable case are shown in Fig. 6. As for comparison, the natural frequency of a linear PEH, an independent beam A, is presented. It can be obtained from Fig. 6 that the natural frequencies of the magnetically coupled dual-beam PEH cover a wide range. It indicates that the proposed system has the potential to operate under excitations with a board spectrum. Particularly, the natural frequency of the system is almost equal to that of the uncoupled piezoelectric cantilever when D = 30 mm due to weak magnetic coupling. As the distance between the two magnets decreases, the natural frequency of the system has a trend of decrease. For example, when D is decreased from 30 to 15.75 mm, the natural frequency will decrease to zero gradually. In addition, the system will be in an annular stable state with negative stiffness characteristics at the location of (X, Y) = (0, 0) when the distance D is smaller than the critical value D 0 .

Harmonic balance analysis
In this section, the harmonic balance method is adopted to solve the analytical responses of the magnetically coupled dual-beam PEH. Firstly, the nonlinear restoring force is expanded by the Taylor series to three orders. Substituting Eq. (13) into Eq. (5) and the governing equations of the dual-beam system are given , the base accelerations a x ¼ A x cos Xt ð Þ and a y ¼ A y cos Xt ð Þ. A x and A y are the amplitudes of the base accelerations in x-and y-directions. The governing Eq. (15) of motion can be normalized as where the system parameters in Eq. (16) are given which represent the damping coefficients of beams, the natural frequencies of the beams, and the electromechanical coupling coefficients, respectively. It can be obtained from Eq. (16) that the vibration motions of the two beams are coupled together through the coefficients k a1 , k a2 , k b1 ; and k b2 . Moreover, the coupling features two-dimensional characteristics with strong nonlinearity. In the following analysis, the Harmonic balance method is adopted to solve the nonlinear equations. The harmonic balance method has been widely utilized to predict the dynamic responses of nonlinear PEHs [40][41][42][43]. Here, we derive the appropriate analytical solutions to Eq. (15) using the harmonic balance method by considering cases in monostable and annular stable states with different shapes of potential wells. The solutions of Eq. (16) are assumed to have the form as Substituting Eqs. (18) and (19) into Eq. (16), neglecting high-order harmonics and balancing the terms of sinðxtÞ and cosðxtÞ, we obtain

Discussion of the simplified magnetic force
In the theoretical analysis, the nonlinear magnetic forces are simplified into a nonlinear term with a cubic stiffness, as presented in Eq. (13). It is necessary to evaluate this simplification. Here, Fig. 7 presents a comparison between the analytical magnetic forces in Eq. (5) and the simplified ones in Eq. (13). It can be obtained that the analytical magnetic force and the simplified one have minor differences in the vicinity at (X, Y) = (0, 0). The discrepancy between the analytical magnetic force and the simplified one becomes larger when the displacements of the two cantilevers X and Y increase. Therefore, all the numerical simulations in this article are performed based on Eq. (5), i.e., the one with analytical magnetic force, for better accuracy. In addition, the numerical simulation can be utilized to validate the results calculated using the harmonic balance method using the approximated magnetic force.
By comparing the frequency responses of the system with analytical magnetic force with that of the system with the simplified force (the one with cubic nonlinearity), it is found that the results have minor differences from each other under small displacement, as shown in Fig. 7. On the other hand, the results of Eq. (15) may not be capable to represent that of Eq. (5) and may result in errors in the analytical solutions with large displacements of the two beams. Notably, the displacement range of ± 10 mm shown in Fig. 7 fits the analysis of the system dynamics in the monostable state or the annular stable state with small amplitude excitations. The responses of the system with strong nonlinearity under large amplitude excitations are beyond the scope of this article.

Case studies and discussions
The nonlinear dynamics and energy harvesting performance of the dual-beam PEH are numerically simulated and analytically analyzed in this section. Here, we consider three cases, namely monostable with weak nonlinearity, monostable with strong nonlinearity, and annular stable states of the system by changing the distance between the two magnets D. Forward and backward frequency sweeping responses are simulated using the Runge-Kutta method in MATLAB. And the harmonic balance analysis is adopted to demonstrate the dynamics of the system. The parameters adopted in the following analysis are presented in Table 1. Without loss of generality, responses of the system under x-direction excitation are focused on considering the symmetry of the PEH. For comparison, the frequency responses of a conventional nonlinear single-degree-of-freedom (SDOF) PEH are calculated. The conventional nonlinear PEH has the same configuration as that of the proposed dual-beam one, while the freedom of beam B is fixed. In this case, the potential energy of the SDOF PEH is the cross section of the two-dimensional potential well when Y = 0, as shown in Fig. 8.

Monostable state with weak nonlinearity
In this section, we first analyze the characteristics of the system in a monostable state with weak nonlinearity. Here, we calculated the frequency response of the system with distances between the two magnets from 30 to 18 mm. The frequency responses of the displacement and voltage output of beam A are shown in Fig. 9.  Figure 9 shows the analytical and simulation results of the dual-beam PEH under x-direction excitation with A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 . Forward and backward frequency sweep is adopted to analyze the performance of the system. Here, the displacements between the two magnets are chosen as D = 30, 20, and 18 mm for illustration. Besides, the frequency responses of displacements and voltage output performance of a linear system, i.e., without magnetic coupling between the two cantilevers, are presented for comparison. It can be obtained from Fig. 9 that the frequency responses of the proposed PEH have almost the same performance as that of the linear one when D = 30 mm. In this case, the two cantilevers have weak magnetic coupling and the proposed PEH is almost a linear one. Moreover, the natural frequency of the proposed system with magnetic coupling (25.5 Hz) is slightly lower than that of the linear one (26 Hz). This phenomenon follows the trend of the resonant frequency reduction effect of the magnetic coupling, as illustrated in Fig. 6. Moreover, the resonant frequency of the system decreases to 21.9 Hz when the distance D decreases to 20 mm. The trend agrees well with the theoretical predictions presented in Fig. 6. In addition, the dynamics of the PEH have a slight hardening effect when D is reduced to 18 mm. This is because of the increased nonlinearity of the system because of the reduced distancer D. Moreover, the forward and backward frequency sweep responses have a minor difference from each other. It can also be obtained that the voltage outputs of the PEH increase by decreasing the distance D. Moreover, the simulation results agree with the analytical ones well.
A key influence factor is the amplitude of excitation. Here, we calculate the analytical frequency response of the system with D = 20 mm and 18 mm under excitations with amplitudes of 1, 2, and 3 m=s 2 , as shown in Fig. 10. Figure 10a shows the frequency response of voltage output when D = 20 mm. It can be obtained that an increase in the amplitude of excitation can effectively enlarge the output voltage. Besides, the increased excitations can also trigger the hardening effect of the system. A slight hardening effect can be observed when A x ¼ 3 m=s 2 . As for comparison, a significant hardening effect can be observed when D = 18 mm and A x ¼ 3 m=s 2 , as shown in Fig. 10b. In addition, the system with D = 18 mm has a peak voltage output of 73.6 V, which is larger than the case with D = 20 mm (66.5 V). In general, the PEH exhibits strong distance D-dependent characteristics. Changing the distance between the permanent magnets can effectively tune its resonant frequency. In addition, the voltage output of the PEH increases as the frequency decreases. This phenomenon demonstrates a promising feature of the system for applications in different vibration scenarios.

Monostable state with strong nonlinearity
In this section, we proceed to investigate the feature of the system in monostable states with large nonlinearity. In this state, the potential well of the system has only one valley at the central point. Here, cases are considered where the distances between the two magnets are 16 and 16.5 mm, i.e., the system has a larger nonlinearity. The frequency responses of the displacements and voltage outputs of beam A and beam B are shown in Fig. 11. Forward and backward frequency sweep responses are calculated analytically and numerically.
It can be obtained from Fig. 11 that the system with D = 16.5 mm exhibits a strong hardening effect in the forward and backward frequency sweeping. The frequency-response curves were obtained from harmonic balance analysis and numerical simulation under A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 . Analytically, the responses of the beam have multiple solutions (with different response amplitudes) in the forward and backward frequency sweeping. This phenomenon follows that of a conventional monostable PEH. Besides, the forward frequency sweeping response has a larger frequency range and voltage output than that of the backward frequency sweeping. The effect stems from the increased nonlinearity of the magnetic coupling. Moreover, the bandwidth of the system is broadened and the amplitudes of displacement/voltage outputs are also increased compared with the responses of the system in a monostable state with weak nonlinearity. For example, the bandwidth of the proposed PEH is 4.1 Hz under excitation of A x ¼ 2 m=s 2 . In addition, the motion of beam A induces a larger amplitude oscillation of beam B in the vicinity of resonances. And the system dynamics of beam B also have multiple solutions in the forward and backward frequency sweeping. On the other hand, the peak voltage outputs of beam B have the same amplitudes in the forward and backward frequency sweep.
We further reduce the distance between the two magnets to D = 16 mm. Similarly, forward and backward sweep response curves were obtained from harmonic balance analysis and numerical simulation under A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 , as shown in Fig. 12. It can be obtained the responses of the beam have multiple solutions in the forward and backward frequency sweep. In addition, the response has a sharp peak in the vicinity of the jumping point. Besides, the frequency responses of the system with D = 16 mm have peak voltage output (48.1 V) smaller than that of the case with D = 16.5 mm (65.8 V) due to the increased nonlinearity and reduced amplitude of motion. Besides, the reduced distance between the two magnets has decreased the frequency of peak voltage output from 14.7 to 10.5 Hz. Moreover, the motion of beam A has also induced the oscillation of beam B. And both the forward and backward frequency sweep responses of beam B have a jumping phenomenon.
For comparison, the frequency responses of a conventional monostable PEH are calculated by the harmonic balance method and numerical simulation, as shown in Fig. 13. The conventional nonlinear PEH has the same configuration as that of the proposed dual-beam one while the freedom of beam B is fixed. Figure 13 shows the responses of the conventional SDOF PEH in a monostable state under excitations of A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 with D = 16.5 mm and D = 16 mm, respectively. It can be obtained that the strong hardening effect can be observed in the responses of the conventional PEH for the two cases with D = 16.5 mm and D = 16 mm. The analytical predictions agree well with that of the numerical simulations. In addition, the amplitude of the conventional PEH has a slightly larger voltage output. For instance, the conventional monostable PEH has a peak voltage output of 66.8 V, while the value is 65.7 V (Beam A) and 34 V (Beam B) for the proposed dualbeam PEH when D = 16.5 mm. Besides, the frequency responses of the conventional one with a hardening effect cover a larger frequency range than that of the proposed dual-beam one. The magnetic coupling between the two beams has significantly reduced the frequency of the jumping point. In particular, the jumping point of the proposed PEH is 10.5 Hz, while it is 13 Hz for the conventional one when D = 16 mm. That is, the proposed system has a narrower bandwidth than that of the conventional ones in the monostable state with strong nonlinearity. For further demonstration of the proposed dual-beam PEH, the frequency responses under excitations with different amplitudes are performed, as shown in Fig. 14. Figure 14 illustrates the analytical responses of the dual-beam PEH in the monostable state and the conventional monostable PEH with D = 16.5 mm.
Here, forward frequency sweep responses were calculated under excitations with different amplitudes. It can be obtained that the operational bandwidth of beam A becomes wider with a higher voltage output as the amplitude of excitation increases. Moreover, the proposed PEH and the conventional one have similar responses when A x = 1 m=s 2 . This is because beam A has a small amplitude motion and weak magnetic coupling between the two beams. Hence, beam B has vibration with minor amplitude in such a case, as shown in Fig. 14b. Moreover, both the responses of beams A and beam B have a significant hardening effect when the amplitudes of excitation increase.
It can be obtained that the proposed dual-beam PEH in a monostable state has voltage output with a smaller amplitude and weaker hardening effect compared to that of the conventional PEH. The reduced amplitude can be attributed to the coupling between the two beams. This is because the hardening effect is induced by a dramatic magnetic force and stiffness shifting due to the large amplitude resonance motion. On the other hand, the proposed dual-beam PEH has motion in twodimensional space. Such a two-dimensional motion between the two magnets than that of the conventional monos-

table PEH
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi . In such a case, the range of magnetic force and stiffness shifting of the proposed system would be much smaller than that of the conventional one. Therefore, the proposed system has a hardening effect weaker than that of the conventional one. For further illustration, the real-time responses of the proposed PEH are presented at the response peak of the proposed PEH, as shown in Fig. 15. Figure 15 presents the time domain responses of beam A and beam B of the proposed dual-beam PEH in a monostable state. Here, excitation with A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 is selected for demonstration. The distances between the two magnets are D = 16.5 mm and D = 16 mm. It can be obtained that the excitation in the x-direction induced large amplitude vibration of beam B in the y-direction due to the magnetic coupling. Besides, the magnetic coupling between the two beams increases rapidly with the reduced magnet distance. Therefore, the displacement of beam B has a larger displacement amplitude when D = 16 mm than when D = 16.5 mm. It can also be obtained that the phase difference between X and Y is closer to p=2 when the distance D is reduced. In such a case, the distance between the two magnets would have little variations, and the distance ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi x 2 þ y 2 þ D 2 p between the magnets is almost a constant one, yielding minor variations of the equivalent stiffness of the system. Therefore, the proposed dual-beam PEH has a weaker hardening effect in the case with D = 16 mm than that with D = 16.5 mm. Besides, reducing the value of D can enhance the coupling between beams A and beam B. And a strong mechanical energy interchange can be found between the two beams, resulting in additional energy dissipation and reduced voltage output.

Annular stable state
In this section, we analyze the characteristics of the PEH in an annular stable state. The conventional bistable PEH has two stable points and a potential barrier between them. In the annular stable state, the potential well of the system has a circular-shaped potential valley in two-dimensional space, as illustrated in Fig. 2c and d. In the following analysis, the distance between the two magnets is chosen as 14.8 mm. The responses of the displacements and voltage outputs of beam A and beam B are shown in Fig. 16. Figure 16 presents the simulation results of the dual-beam PEH with D = 14.8 mm. Here, the amplitude of excitation is chosen as A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 . The displacement and voltage output of beam A are presented in Fig. 16a and b, and the displacement and voltage output of beam B are presented in Figs. 16c and 15d. Different from that in a monostable state, the frequency response of the dual-beam PEH has a softening effect. Most importantly, the displacement of beam A in the annular stable state has an amplitude significantly larger than that in the monostable state, as shown in Fig. 16a. For example, the dual-beam PEH has a peak voltage output of 103.5 V at 19.1 Hz in the annular stable state. And the bandwidth of the system is 3.1 Hz here. Meanwhile, the proposed PEH has a peak voltage output of 48.1 V at 10.5 Hz in the monostable stable state when D = 16 mm, i.e., the voltage output amplitude is increased by 215% in the annular stable states. Moreover, there is a strong coupling between beam A and beam B. Hence, beam B has a peak voltage output of 93.4 V at 19.9 Hz, at the same level as that of beam A. In addition, the forward and  Fig. 17. The conventional bistable PEH has the same parameters as that of the proposed one, while the displacement of beam B is fixed. Figure 17 shows the responses of the conventional SDOF PEH in the bistable state with D = 14.8 mm.
Here, the amplitude of excitation is chosen as A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 . It can be obtained from Fig. 17 that the conventional PEH has responses much smaller than that of the proposed PEH. In particular, the conventional PEH has a peak response of 37.6 V, while the proposed one has a peak response of 103.5 V. That is, the proposed dual-beam PEH has peak voltage output improved by 275% in this case study. Besides, the bandwidth of the conventional bistable PEH is 2.6 Hz. The conventional PEH has intra-well motion only due to the small amplitude excitation of A x ¼ 2 m=s 2 . In other words, the conventional PEH cannot jump over the large potential well as it received limited energy from the external vibration source.
It should be noted that beam B in the dual-beam PEH has large amplitude oscillation simultaneously. For further demonstration of the proposed dual-beam PEH, the frequency responses under excitations with different amplitudes are performed, as shown in Fig. 18. Figure 18a and b illustrates the numerical frequency responses of the dual-beam PEH in the annular stable state with D = 14.8 mm. Besides, the frequency response of a conventional bistable PEH is presented for comparison. It can be obtained that both beam A and beam B of the proposed dual-beam PEH have large amplitude voltage outputs in a wide frequency range. Moreover, the operational bandwidth of beam A and beam B is expanded with a higher voltage output as the amplitude of excitation increases. Meanwhile, the operational bandwidth and voltage outputs of the conventional bistable PEH increase slightly as the amplitude of excitation increases. It is worth noticing  that the proposed dual-beam PEH has a rectangularshaped large amplitude over a wide frequency range, i.e., the peak performance can be maintained throughout the frequency range. This is because the dual-beam PEH has a large amplitude circular-shaped motion due to the two-dimensional potential well. As for comparison, the conventional bistable PEH has a triangleshaped voltage output, i.e., with a peak voltage output at a single frequency point only.
In addition, the proposed dual-beam PEH in an annular stable state has significant advantages over the conventional one in the aspect of the amplitude of voltage output. For example, beam A outputs 103.5 V, 123.2 V, and 124.6 V under excitation of A x ¼ 2 m=s 2 , A x ¼ 3 m=s 2 , and A x ¼ 4 m=s 2 , respectively. Here, the bandwidths of the voltage output of beam A are 3.1 Hz, 3.9 Hz, and 5.6 Hz, respectively. Besides, the corresponding voltage outputs of beam B are 93.4, 105.1, and 111.2 V under the excitations, respectively. On the other hand, the voltage output of the conventional bistable PEH outputs 37.7 V, 40.9 V, and 43.1 V under excitation of A x ¼ 2 m=s 2 , A x ¼ 3 m=s 2 ; and A x ¼ 4 m=s 2 , respectively. That is, beam A has a performance 2.74, 3.01, and 2.89 times better than that of the conventional bistable one under excitation of A x ¼ 2 m=s 2 , A x ¼ 3 m=s 2 ; and A x ¼ 4 m=s 2 , respectively.
It can also be obtained that the proposed dual-beam PEH itself has a better performance in the annular stable state than that in the monostable state. For instance, the proposed dual-beam PEH has a triangleshaped voltage output with a peak of 51.8 V under excitation of A x ¼ 2 m=s 2 in the monostable with weak nonlinearity, as shown in Fig. 10. Besides, the proposed PEH has a triangle-shaped voltage output with a peak of 65.1 V under excitation of A x ¼ 2 m=s 2 in the monostable with large nonlinearity, as  Fig. 14. That is, the annular stable state can effectively enhance the voltage output as well as the working bandwidth due to the two-dimensional circular-shaped potential well.
To study the characteristics of the features of the oscillation of the proposed PEH, we presented the time domain responses and the corresponding phase of beam A and beam B in the vicinity of the resonant frequency. Figure 19 shows the time domain displacements, phase portrait, and potential energy in the twodimensional space of beam A and beam B, respectively. The distance between the two magnets is set as D = 14.8 mm. Figure 19d shows the top view of Fig. 19c, and the blue curve is the isotropic potential line of the system. And the amplitude of excitation is chosen as A x ¼ 2 m=s 2 and A y ¼ 0 m=s 2 . Besides, here we focus on analyzing the time domain response of the system at the peak of voltage output where f = 19.5 Hz, the frequency point with peak voltage output. It can be obtained from Fig. 19a that both beam A and beam B have large amplitude motions. In particular, the two beams have a combined high-frequency/small amplitude motion and low-frequency/large amplitude motion. Here, the high-frequency/small amplitude motion comes from the motion in the radial direction of the circular-shaped potential well, as shown in Fig. 19c and d. Moreover, the low-frequency/large amplitude motion stems from the motion in the tangential direction of the circularshaped potential well. Besides, the amplitudes of the motion of beam A and beam B oscillate from one region to the other periodically as shown in Fig. 19a and b. This phenomenon can be attributed to the combined motion in both the radial direction and the tangential direction of the circular-shaped potential well. In addition, the low-frequency component of the responses stems from the motion moving around the potential barrier.
It is worth mentioning that the conventional SDOF PEH may have a phase portrait with a similar shape as that presented in Fig. 19b. However, a jumping motion between two potential wells is necessary for the conventional bistable PEH which requires large amplitude excitations. On the other hand, the motion of the proposed PEH is indeed in a circular manner in c the potential energy of annular stable dual-beam PEH; and d top view of (c) the annular potential well. In other words, the proposed system bypassed the potential barrier rather than jumping through it. Therefore, the proposed system features advantages of operation under small amplitude excitations.
To further demonstrate the features of the system under small amplitude excitation, the responses of the system under A x ¼ 1 m=s 2 are calculated, while the other parameters are kept the same, as shown in Fig. 20. Figure 20 shows the responses of the dual-beam PEH in the annular stable state with D = 14.8 mm.
Here, the amplitude of the excitation is reduced to A x ¼ 1 m=s 2 and A y ¼ 0 m=s 2 . Figure 20a-d shows the time domain displacements, the phase portrait, and the potential energy in the two-dimensional space of beam A and beam B, respectively. Similarly, Fig. 20d shows the top view of Fig. 20c. It can be obtained from Fig. 20a and b that beam A exhibits large amplitude oscillation, while beam B vibrates in a small range of displacement. In particular, beam A oscillates in two states in the phase portrait, while the motion of beam B is limited to one state. This is because beam B does not have enough mechanical energy obtained under this scenario. Although the beams do not need to jump over the potential barrier, beam B still needs extra mechanical energy for the bending motion in the circular-shaped potential well. Nevertheless, beam A is capable to oscillate with large amplitude motions in the circular-shaped potential well, indicating significant voltage outputs. Figure 21 shows the simulated responses and phase portrait of the dual-beam PEH in the annular stable state with D = 14.8 mm. The operational frequency is chosen as f = 15 Hz, i.e., in off-resonance states. In the case studies, two excitations in the xdirection are considered that A x ¼ 1 m=s 2 and A x ¼ 2 m=s 2 . It can be obtained from Fig. 21a and b that beam A has limited vibration amplitudes in both cases, i.e., exhibiting intra-well oscillations. It is similar to the intra-well oscillations of the conventional SDOF bistable PEH. This is because the offresonance motion can provide limited energy for the circular-shaped motion of the beams. And the system moves at steady-state points of (X, Y) = (5.4 mm, 0). c the potential energy of annular stable dual-beam PEH; and d top view of (c) Meanwhile, beam B has a vibration with minor amplitudes, as shown in Fig. 21c-e.
It is worth emphasizing that the excitations with amplitudes of A x ¼ 2 m=s 2 and A x ¼ 1 m=s 2 are very small. In other words, a conventional bistable PEH cannot have inter-well motion under such an excitation. A typical bistable PEH may require excitation with amplitude up to 10 m=s 2 [23,36] and 0.8 g [20]. To further demonstrate this point as well as the advantages of the proposed system, case studies of a conventional bistable PEH are carried out. The conventional bistable PEH has the same configurations and parameters. The displacements of beam A are shown in Figs. 22 and 23. Here, the distance between the two magnets is kept as D = 14.8 mm. Besides, responses under excitations with amplitudes of A x ¼ 1 m=s 2 , 2 m=s 2 ; and 14.4 m=s 2 are analyzed. Figure 22 shows the simulated responses and phase portraits of the bistable PEH with D = 14.8 mm. The operational frequency is chosen as f = 19.5 Hz. In the case studies, two excitations in the x-direction are considered that A x ¼ 1 m=s 2 and A x ¼ 2 m=s 2 . It can  Fig. 22a and b that the conventional bistable PEH has limited vibration amplitudes in both cases, i.e., exhibiting intra-well oscillations. And the PEH oscillates at steady-state points of (X, Y) = (5.2 mm, 0) with a peak-peak amplitude of 3.3 mm. The inter-well motion of the conventional bistable PEH can be triggered when the amplitude of the excitation is increased to A x = 14.4 m=s 2 , as shown in Fig. 23a. In such a case, the response amplitude of the conventional bistable PEH has increased dramatically to 18.1 mm. It is worth noticing that large amplitude excitations are necessary for the conventional bistable PEH that features large amplitude inter-well motion. On the other hand, the proposed PEH is capable to yield large amplitude motion with low amplitude excitation. A detailed comparison is presented in Fig. 24. Figure 24 shows the simulated responses of the conventional SDOF PEH in a bistable state and the proposed dual-beam PEH in an annular stable state with D = 14.8 mm and A x ¼ 2 m=s 2 at f = 19.5 Hz. It can be obtained that the SDOF PEH in the bistable state only vibrates within one potential well and the vibration amplitude is limited. Meanwhile, the  dual-beam PEH can vibrate over a larger range as it can move along the circular-shaped potential well. Hence, it can be concluded that the proposed dualbeam PEH has significant advantages under small amplitude excitations. It is also worth mentioning that the dual-beam PEH has the potential to harvest vibration energy in two orthogonal directions. In general, the annular stable state enables large amplitude motion of the cantilever under excitations with small amplitudes. Besides, the annular stable state facilitates the oscillation of the beam along the bottom of the circular-shaped potential barrier, allowing the beam to operate in a wide frequency range.

Conclusion
In this manuscript, a new piezoelectric energy harvester with an annular stable potential well was presented for broadband vibration energy harvesting. To demonstrate the advantages of the proposed annular stable harvester, the response behavior of a dual-beam system was studied under excitations with small amplitudes. The two cantilevers are placed orthogonal to each other and be coupled by repulsive magnetic force with a two-dimensional potential well. A model of the magnetic field was utilized to demonstrate the distance-based characteristics of the potential well. And a mathematical model was derived, and the Harmonic balance method was adopted for the analytical solutions of the system dynamics. It is shown analytically and numerically that monostable and annular stable can be configurated with the capability of broadband energy harvesting with the proper assembling of the two cantilevers. In the monostable state, the dual-beam PEH exhibits a hardening effect that is feasible for broadband operation. Moreover, the operational frequency range is reduced with the reduction of the distance between two cantilevers. In the annular stable state, the dualbeam PEH is capable of detours around the potential barrier rather than jumping over it. This feature enables the system to maintain large amplitude voltage outputs throughout a wide spectrum under small amplitude excitations. Case studies were carried out, and it is proved that the proposed dual-beam PEH, with a bandwidth of 3.9 Hz, has a performance 3.01 times better than that of a conventional bistable one under 3 m/s 2 amplitude excitations. The annular potential energy function retains the advantage of large inter-well oscillation in the bistable system for generating high-voltage output over a wide frequency range. And it is not necessary to have a larger excitation to jump through the potential barrier. The proposed method may be useful in applications where large amplitude excitation are unavailable, or when large amplitude output should be maintained over a wide frequency range. Note that the two beams can absorb vibrational energy in different directions. This Data availability All necessary data needed to generate the results obtained in this study have been included in the article.

Declarations
Conflict of interest The author declares that he has no conflicts of interest to report regarding the present study. Appendix A Q i ði ¼ 1; 2; 3; 4; 5; 6Þ in Eqs. (20) and (21) can be expressed as: