Theoretical Study on Widening Bandwidth of Piezoelectric Vibration Energy Harvester with Nonlinear Characteristics

In order to make a piezoelectric vibration energy harvester collect more energy on a broader frequency range, nonlinearity is introduced into the system, allowing the harvester to produce multiple steady states and deflecting the frequency response curve. However, the harvester can easily maintain intra-well motion rather than inter-well motion, which seriously affects its efficiency. The aim of this paper is to analyze how to take full advantage of the nonlinear characteristics to widen the bandwidth of the piezoelectric vibration energy harvester and obtain more energy. The influence of the inter-permanent magnet torque on the bending of the piezoelectric cantilever beam is considered in the theoretical modeling. The approximate analytical solutions of the primary and 1/3 subharmonic resonance of the harvester are obtained by using the complex dynamic frequency (CDF) method so as to compare the energy acquisition effect of the primary resonance and subharmonic resonance, determine the generation conditions of subharmonic resonance, and analyze the effect of primary resonance and subharmonic resonance on broadening the bandwidth of the harvester under different external excitations. The results show that the torque significantly affects the equilibrium point and piezoelectric output of the harvester. The effective frequency band of the bistable nonlinear energy harvester is 270% wider than that of the linear harvester, and the 1/3 subharmonic resonance broadens the band another 92% so that the energy harvester can obtain more than 0.1 mW in the frequency range of 18 Hz. Therefore, it is necessary to consider the influence of torque when modeling. The introduction of nonlinearity can broaden the frequency band of the harvester when it is in primary resonance, and the subharmonic resonance can make the harvester obtain more energy in the global frequency range.


Introduction
The battery is a traditional energy storage device. Many devices (such as intelligent home equipment, implantable medical equipment, environmental monitors) rely on batteries for energy supply. However, the battery cannot achieve long-term energy supply due to limited life, energy storage, and power efficiency and often needs frequent replacement or maintenance. An increasing number of researchers have paid attention to topics such as collecting energy from the environment to eliminate the demand for batteries or extend batteries' life. Vibration energy is typical mechanical energy in the environment. There are many ways to produce it, such as the surrounding environment, wind flow, fluid movement, mechanical work, and human behavior. Piezoelectric vibration energy harvesters are widely used in research for low power electronic devices such as embedded electronic devices, implantable biomedical devices, wireless sensor nodes, and portable electronic devices due to their simple structure, long life, high power density, and no need for an initial voltage source.
The piezoelectric cantilever beam is the most commonly studied and discussed structure in a piezoelectric vibration energy harvester because of its simple structure and convenient modeling. In the linear piezoelectric cantilever beam structure, the mass is usually added to the free end of the cantilever beam [1], which can reduce and adjust the system's resonance frequency and improve the mechanical response and output power under low-frequency excitation. However, the main limitation of linear vibration energy harvester is that the natural frequency is single, and the resonance peak is very narrow. When the excitation frequency deviates from the resonant frequency of the harvester, very little power can be generated [2]. The vibration frequency of the environment usually fluctuates in a specific range. In order to obtain better output performance and maintain higher acquisition efficiency, a broadband harvester is required. Introducing nonlinearity into the harvester can broaden the frequency bandwidth of the system, increase the mechanical response amplitude and power output. Nonlinear spread spectrum technology makes use of nonlinear stiffness [3,4], the piecewise linear effect generated by the collision process [5,6], and multi-steady-state motion [7,8] to deflect the amplitude-frequency curve of the energy harvester, ensure a large amplitude in a wide frequency range, achieve the effect of broadening the frequency band, and improve the acquisition efficiency of the energy harvester.
The realization of the multi-steady state of the piezoelectric vibration energy harvester mainly uses the magnetic field force provided by permanent magnets. A different number of steady states can be obtained by using different arrangement modes of permanent magnets. Erturk et al. [9] propose a bistable piezoelectric vibration energy harvester with two stable positions. The device consists of a piezoelectric cantilever beam and three permanent magnets. A permanent magnet is installed at the tip of the cantilever beam, and two permanent magnets are fixed near the tip of the cantilever beam. The experimental results show that compared with the linear harvester, the system has a wider effective frequency band and higher output power. Since then, bistable piezoelectric vibration energy harvesters have been studied by many researchers. Zhou et al. [10] propose a triple steady-state piezoelectric vibration energy harvester, which consists of a piezoelectric cantilever beam with a permanent magnet at the tip and two external permanent magnets. Compared with the bistable energy harvester with deeper potential wells, the triple steady-state arrangement has a shallow potential well, which allows the harvester to easily produce large amplitude vibration between wells and higher energy output on a broader frequency range.
For the nonlinear piezoelectric energy harvester, the forward sweep can obtain a wider effective bandwidth than the reverse sweep when the sweep is conducted. Therefore, the nonlinear harvester cannot always obtain a higher amplitude as the external excitation frequency varies. Qing et al. [11] propose a lambda-shaped piezoelectric energy harvester, which can produce much higher power than the cantilever piezoelectric energy harvester but still has the shortcomings of the linear energy harvester. Li et al. propose [12] a generalized multi-mode piezoelectric energy harvester. The model consists of the main cantilever beam and multiple branches. The effective bandwidth of the harvester can be adjusted by changing the branch parameters and tip mass. It can overcome the drawbacks of the linear piezoelectric energy harvester and avoid the drawbacks mentioned above of nonlinear energy harvesters. Sun et al. [13] propose a horizontal asymmetric U-shaped piezoelectric energy harvester model by combining multi-mode and nonlinear force. The asymmetric U-shaped structure gives the model multi-mode properties, while the nonlinearity introduced by magnetic force improves the energy output of the harvester, reduces the resonant frequency, and broadens the effective bandwidth. The above studies improved the efficiency of energy acquisition by changing the structure of the piezoelectric energy harvester. However, the complex structure may have problems in practical application and model manufacturing.
The frequency sweep of a nonlinear piezoelectric energy harvester cannot reveal all the steady-state behaviors it can achieve, and frequency sweep does not reveal the subharmonic motion. Syta et al. [14] studied the subharmonic vibration of a bistable harvester at a specific frequency, and the system showed different motion behaviors under different initial conditions. Huguet et al. [15][16][17] propose a complete subharmonic orbit analysis to predict the contribution to the global bandwidth of bistable energy harvesters and study the effects of different parameters on the response through experiments. Syta and Huguet mainly studied subharmonic motion by numerical simulation and experimental methods but did not study the influence of subharmonics on broadening the frequency band under different external excitations. Huguet also used the harmonic balance method to study the robustness of subharmonic motion, but for strongly nonlinear systems, sufficient harmonic terms are needed to obtain more accurate asymptotic solutions. Therefore, this paper theoretically investigates how to make full use of nonlinear characteristics to broaden the bandwidth of the energy harvester of the piezoelectric cantilever beam structure.
The theoretical method for solving strongly nonlinear vibration has always been a difficulty. In order to solve this problem, Wang et al. [18] proposed a single-degree-of-freedom complex dynamic frequency (CDF) method, applied it to the study of strongly nonlinear vibration energy harvester, and experimentally verified the method's effectiveness. This paper extends this method based on this theory and studies the above strongly nonlinear vibration problem.
In modeling the piezoelectric cantilever beam harvester, the magnetic repulsion between permanent magnets is generally considered [19,20]. The torque between permanent magnets also affects the bending of piezoelectric cantilever beams, thereby affecting the piezoelectric effect. Kim et al. [21] established a model considering the torque between permanent magnets but did not analyze the influence of torque. This paper discusses the influence of torque on the piezoelectric cantilever vibration energy harvester.
In this paper, the energy harvester of the piezoelectric cantilever beam structure is studied. The permanent magnets provided nonlinearity for the harvester, and the magnetic repulsion and torque between them are considered. In order to achieve a bistable state, the distance between the free-end permanent magnet of the cantilever beam and the fixed permanent magnet is determined according to the static analysis. The approximate analytical solution and amplitude-frequency relationship of the strongly nonlinear harvester's steady-state motion are obtained using the CDF method. The broadening effect of primary resonance and subharmonic resonance on the bandwidth of the harvester is compared and analyzed.
The rest of this paper is organized as follows. Section 2 gives the piezoelectric cantilever beam dynamics model considering the torque between permanent magnets. Section 3 outlines the static analysis of the system. In Section 4, the approximate analytical solutions and amplitude-frequency relations of the system's primary and subharmonic resonance are obtained using the CDF method. Section 5 analyzes the system's dynamic characteristics and compares the differences in the energy obtained by the primary resonance and the 1/3 subharmonic resonance in the global operating frequency range and the variation of the subharmonic response with different external excitations. Section 6 summarizes the full-text results.

Theoretical Model
The piezoelectric material as an energy transfer device directly affects the collection efficiency of the harvester. Currently, commonly used piezoelectric materials [22][23][24] include piezoelectric ceramics, piezoelectric polymer, and piezoelectric composites. Piezoelectric ceramics have a high electromechanical coupling coefficient and low cost but are brittle and easily fractured, piezoelectric polymers are flexible but have low electromechanical coupling coefficients, and piezoelectric composites combine the advantages of piezoelectric ceramics and piezoelectric polymers with excellent piezoelectric properties and strong flexibility. Therefore, in this paper, piezoelectric composites are used as piezoelectric materials of the piezoelectric energy harvester.
The mechanical model of the piezoelectric cantilever beam is shown in Figure 1. L is the length of the composite beam, which consists of a substrate beam and two piezoelectric plates. The materials of the substrate and piezoelectric plate are beryllium bronze QBe1.9 and FMC M2807 P2. A layer of the piezoelectric plate is pasted on the substrate's top and bottom surfaces, and the substrate is rigidly connected to the piezoelectric plates without considering the influence of the multilayer beam. One end of the composite beam is fixed on base C, and the other is free. Permanent magnet A is fixed to the free end of the composite beam, and permanent magnet B is fixed to another side of base C. The two permanent magnets are mutually exclusive and horizontally aligned. When the system is in equilibrium, the gravity of the composite beam and permanent magnet A has no effect on the deformation of the piezoelectric cantilever beam. P(t) is the basic motion; x is the axial spatial coordinate; t is the time coordinate; z(x,t) is the transverse deformation of the position on the cantilever beam at any time; and F x , F z , and M c are the magnetic force in the x-direction, the magnetic force in the z-direction, and the torque in the y-direction of the permanent magnet A under the action of the permanent magnet B, respectively. The physical parameters and material properties of the piezoelectric cantilever beam vibration system are shown in Table 1. and strong flexibility. Therefore, in this paper, piezoelectric composites are used as piezoelectric materials of the piezoelectric energy harvester. The mechanical model of the piezoelectric cantilever beam is shown in Figure 1. L is the length of the composite beam, which consists of a substrate beam and two piezoelectric plates. The materials of the substrate and piezoelectric plate are beryllium bronze QBe1.9 and FMC M2807 P2. A layer of the piezoelectric plate is pasted on the substrate's top and bottom surfaces, and the substrate is rigidly connected to the piezoelectric plates without considering the influence of the multilayer beam. One end of the composite beam is fixed on base C, and the other is free. Permanent magnet A is fixed to the free end of the composite beam, and permanent magnet B is fixed to another side of base C. The two permanent magnets are mutually exclusive and horizontally aligned. When the system is in equilibrium, the gravity of the composite beam and permanent magnet A has no effect on the deformation of the piezoelectric cantilever beam. P(t) is the basic motion; x is the axial spatial coordinate; t is the time coordinate; z(x,t) is the transverse deformation of the position on the cantilever beam at any time; and Fx, Fz, and Mc are the magnetic force in the x-direction, the magnetic force in the z-direction, and the torque in the y-direction of the permanent magnet A under the action of the permanent magnet B, respectively. The physical parameters and material properties of the piezoelectric cantilever beam vibration system are shown in Table 1.    Due to the sizeable length/thickness ratio of the beam, this paper only considers the bending deformation of the piezoelectric cantilever beam; regardless of the shear deformation and the influence of the rotational inertia of the section around the neutral axis, the piezoelectric cantilever beam is regarded as an Euler-Bernoulli beam. The nonlinear magnetic force at the free-end boundary is regarded as a concentrated load, as shown in Figure 2. 0 Due to the sizeable length/thickness ratio of the beam, this paper only consid bending deformation of the piezoelectric cantilever beam; regardless of the shear mation and the influence of the rotational inertia of the section around the neutr the piezoelectric cantilever beam is regarded as an Euler-Bernoulli beam. The no magnetic force at the free-end boundary is regarded as a concentrated load, as sh Figure 2. The influence of temperature and residual stress on the piezoelectric cantileve is not considered in the modeling. According to the Hamilton principle, the gov equation of the harvester can be written as where T is the system's kinetic energy, U is the system's potential energy, W is the of the magnetic repulsion and torque, and δ is the variational symbol. According to 1, the kinetic energy of the substrate is where s A is the section area of the substrate, and ρ s is the substrate's material d The kinetic energy of the piezoelectric plate is where p A is the section area of the piezoelectric plate, and ρ p is the piezoelectric material density. The kinetic energy of free-end permanent magnet A is where a m is the mass of permanent magnet A, and a J is the rotational inertia manent magnet A. The bending strain energy of the substrate is where s E is the substrate's Young' s modulus, and s I is the section inertia mom the substrate, 3 / 12 The bending strain energy of the piezoelectric plate is The influence of temperature and residual stress on the piezoelectric cantilever beam is not considered in the modeling. According to the Hamilton principle, the governing equation of the harvester can be written as where T is the system's kinetic energy, U is the system's potential energy, W is the energy of the magnetic repulsion and torque, and δ is the variational symbol. According to Figure 1, the kinetic energy of the substrate is where A s is the section area of the substrate, and ρ s is the substrate's material density. The kinetic energy of the piezoelectric plate is where A p is the section area of the piezoelectric plate, and ρ p is the piezoelectric plate's material density. The kinetic energy of free-end permanent magnet A is where m a is the mass of permanent magnet A, and J a is the rotational inertia of permanent magnet A. The bending strain energy of the substrate is where E s is the substrate's Young' s modulus, and I s is the section inertia moment of the substrate, I s = bh 3 s /12. The bending strain energy of the piezoelectric plate is where E p is the piezoelectric plate's Young' s modulus, I p = (4h 2 p + 6h p h s + 3h 2 s )bh p /12 is the inertia moment of the piezoelectric plate, e 31 = E p d 31 is the effective piezoelectric stress constant, b is the width of the substrate, h s and h p are the substrate and piezoelectric thickness, λ(t) is the flux linkage, C p = bε S 33 L p /h p is the capacitance, ε S 33 = ε 31 ε 0 is the dielectric constant of medium, and ε 31 is the relative dielectric constant of the piezoelectric plate. By considering the first-order vibration mode of the piezoelectric cantilever beam, we can set where φ(x) is the vibration mode function of the beam, and η(t) indicates the vibration mode of the beam as a function of the generalized time coordinate. Since the research object is a variable-section beam, the vibration mode function of the beam needs to be expressed in segments, where c i1 ∼ c i4 (i = 1, 2) are determined by the piezoelectric cantilever beam boundary conditions and continuity conditions, and β i represents eigenvalues. Let the virtual work of moment M y0 , M yL and shear Q z0 , Q zL acting on the ends of the cantilever beam at the corresponding virtual displacement on the boundary be where δz( (1) and organized to obtain the system control equation ..
where ξ is the damping ratio. The damping source in the system consists of mechanical damping and air damping, and the specific value of the damping ratio is determined by the experiment; ξ = 0.0178 is taken in this paper. ω is the natural frequency of the system, is the base excitation coefficient, and θ = e 31 b(h p + h s )φ 1 (L p )/2 is the electromechanical coupling terms. The normalization condition is In order to obtain the natural frequency and mode of the system, the Euler-Bernoulli beam transverse vibration equation without damping free vibration [25] is used, where EI(x) is flexural stiffness of composite section, and m(x) is mass per unit length of composite beam. Substituting Equation (7) into Equation (15) yields where ω is the natural frequency of the system. Equation (16) can be rewritten as where β 4 = m(x)ω 2 /EI(x). For the segmented beam whose vibration mode function is shown in Equation (8), we have fixed-end boundary conditions: free-end boundary conditions: and continuity conditions: By combining the eight Equations in (19)~(21), we can obtain the linear homogeneous equations for c 11 ∼ c 14 and c 21 ∼ c 24 . In order to obtain the non-zero solution, the determinant of the coefficient matrix is set to zero, and the first-order natural frequency of the system can be obtained as ω = 154 rad/s. By using the above normalization conditions, all the unknown variables c 11 ∼ c 14 and c 21 ∼ c 24 can be obtained, and the modal function of the piezoelectric cantilever beam is

Static Analysis
In the study of nonlinear vibration energy harvesters, many researchers [26][27][28] provide nonlinear force for the vibration energy harvester by introducing magnetic force. The magnetic dipole method and the magnetization current method are common methods to calculate magnetic force. Tan et al. [29] prove that both methods produce errors when the permanent magnet spacing is sufficiently small. However, the magnetic dipole method results in closer calculations to experimental measurements than the magnetization current method. This section firstly uses the magnetic dipole method to derive the magnetic force and torque expressions [30]. Secondly, the equilibrium point of the harvester is determined. Finally, the effect of torque on the harvester is analyzed.

Magnetic Force and Torque
The structure of the piezoelectric cantilever beam is shown in Figure 1, where the permanent magnets A and B are simplified as magnetic dipoles A and B, respectively. The magnetic induction intensity generated by dipole B at the location of dipole A is given by where µ 0 is the vacuum permeability. m B is the magnetic moment of magnetic dipole B, and the size of the magnetic moment is related to the volume of the permanent magnet.
where m A is the magnetic moment of the magnetic dipole A. According to the vector differential principle, the following can be obtained where m A = m AmA , m B = m BmB and r = rr.m A ,m B , andr are unit vectors. Figure 3 shows the geometric relationship diagram when the free-end displacement of the piezoelectric cantilever beam is z. The magnetic force in the z-direction and the torque in the y-direction can be obtained as follows

Equilibrium Points of the System
By analyzing the equilibrium point and stability of the autonomous system corresponding to Equations (11) and (12), the static bifurcation characteristics of the piezoelectric cantilever beam system are obtained. Fourier expansion omits higher-order terms for magnetic force and torque applied to the free end of a piezoelectric cantilever beam ( ) where a 1 , a 2 , b 1 , and b 2 are the corresponding Fourier expansion coefficients, respectively. (11) and (12), we obtain

Equilibrium Points of the System
By analyzing the equilibrium point and stability of the autonomous system corresponding to Equations (11) and (12), the static bifurcation characteristics of the piezoelectric cantilever beam system are obtained. Fourier expansion omits higher-order terms for magnetic force and torque applied to the free end of a piezoelectric cantilever beam where a 1 , a 2 , b 1 , and b 2 are the corresponding Fourier expansion coefficients, respectively. Letting η, and z 3 = v in Equations (11) and (12), we obtain where The system has three fixed points: (0,0,0), ± (α 1 − ω 2 )/α 2 , 0, 0 . According to the Routh-Hurwitz criterion, there is a stable zero fixed point when α 1 < ω 2 , and there are two stable nonzero fixed points and an unstable zero fixed point when α 1 > ω 2 , so α 1 = ω 2 is a bifurcation point. By substituting the parameters, it can be determined that when the distance between permanent magnets d < 25 mm, the system produces bistable motion. If the term containing damping does not appear in Equation (32), the energy function of the system is The contour diagram of the energy function L(z 1 , z 2 ) shown in Figure 4 is obtained by assigning different initial values to z 1 and z 2 . There are four different types of curves (or points) in the figure: The contour diagram of the energy function L z z   Figure 5 shows the static bifurcation diagram of the equilibrium from the figure that when the distance between the permanent magn is almost no difference in the equilibrium position in both cases. W effect of torque on the equilibrium position is greater as d increases. In influence of torque on the energy output of the harvester, the effect of voltage and power is analyzed below when d > 20 mm.  Figure 5 shows the static bifurcation diagram of the equilibrium point. It can be seen from the figure that when the distance between the permanent magnets d is small, there is almost no difference in the equilibrium position in both cases. When d > 20 mm, the effect of torque on the equilibrium position is greater as d increases. In order to study the influence of torque on the energy output of the harvester, the effect of torque on the RMS voltage and power is analyzed below when d > 20 mm. from the figure that when the distance between the permanent magn is almost no difference in the equilibrium position in both cases. Wh effect of torque on the equilibrium position is greater as d increases. In influence of torque on the energy output of the harvester, the effect of voltage and power is analyzed below when d > 20 mm.  Figure 6a gives the RMS voltage versus load resistance for the h voltage increases as the load resistance increases and eventually incr Figure 6b shows the corresponding power versus load resistance, wh increases and then decreases with the increasing load resistance. It can Figure 6 that torque affects voltage and power. Figure 7 shows the re the relative error of the output power and the load resistance witho torque. The results show that when the load resistance is small and th permanent magnets is large, not considering that the torque causes a maximum error can reach 18.96% when d = 24 mm.  Figure 6a gives the RMS voltage versus load resistance for the harvester. The RMS voltage increases as the load resistance increases and eventually increases to saturation. Figure 6b shows the corresponding power versus load resistance, where the power first increases and then decreases with the increasing load resistance. It can be observed from Figure 6 that torque affects voltage and power. Figure 7 shows the relationship between the relative error of the output power and the load resistance without considering the torque. The results show that when the load resistance is small and the distance between permanent magnets is large, not considering that the torque causes a large error, and the maximum error can reach 18.96% when d = 24 mm.

Approximate Analytical Solution
Since the system is a strongly nonlinear vibration system, the approximate analytical solution of the system is obtained using the complex dynamic frequency (CDF) method and compared with the multiscale method and the numerical solution.

Approximate Analytical Solution
Since the system is a strongly nonlinear vibration system, the approximate analytical solution of the system is obtained using the complex dynamic frequency (CDF) method and compared with the multiscale method and the numerical solution.

Approximate Analytical Solution
Since the system is a strongly nonlinear vibration system, the approximate analytical solution of the system is obtained using the complex dynamic frequency (CDF) method and compared with the multiscale method and the numerical solution.

When 1/3 subharmonic resonance occurs in systems Equations
Combined with Equation (7), the actual response of the harvester when 1/3 subharmonic resonance occurs can be obtained. Figure 8 shows the amplitude-frequency response diagram of the system, where Figure 8a is the primary resonance. The purple line is the amplitude-frequency relationship obtained with the multi-scale method. The multi-scale method [31,32] has been widely used in the research of energy harvesters, so the iterative process of the multi-scale method is not listed in this paper. The multi-scale method effectively solves weakly nonlinear vibration problems, but accurate results cannot be obtained for strongly nonlinear vibration problems. It can be seen from Figure 8a that for the object studied in this paper, the results obtained using the multi-scale method differ significantly from the numerical solutions. The results obtained using the CDF method are in better agreement with the numerical solution. Therefore, it is feasible to use the CDF method for the theoretical analysis in this paper. As shown from Figure 8a when the excitation frequency is between 8.5 and 17.4 Hz, one frequency corresponds to multiple amplitudes in the system, one unstable solution, and two stable solutions. The larger values represent the inter-well motion in the stable solution, and the smaller values represent the intra-well motion. The trajectories of the response amplitudes during forward and reverse sweep are indicated by arrows 1 and 2 , respectively. The steady-state response at the former excitation frequency during the sweep is the initial condition for the vibration at the latter excitation frequency. It is easily obtained that the steady-state motion amplitude of the system depends on the initial condition in the frequency range of single-frequency multi-amplitude. As shown in Figure 4, the system has two centers and a saddle point. When the initial conditions are different, some orbits in the system surround the center, and the other orbits are around the three fixed points, which leads to the results in Figure 8a. The amplitude-frequency response of the corresponding linear system is shown in Figure 8b. It can be seen from the figure that the linear system does not have a single-frequency multi-amplitude phenomenon, and the forward and reverse sweep results are the same. In Figure 8a, the frequency range of amplitude greater than 4 mm is 0~17.4 and 0~8.5 Hz for the forward sweep and reverse sweep, respectively. In Figure 8b, the frequency range of amplitude greater than 4 mm is 23.5~25.8 Hz. It can be seen that the introduction of strong nonlinearity in the harvester not only reduces the resonant frequency but also expands the effective frequency range by at least 270%. Figure 8c gives the amplitude-frequency response relationship of the nonlinear system when 1/3 subharmonic resonance occurs. It can be seen from the figure that only in a specific frequency range, i.e., 10.1~17.9 Hz, the system is likely to undergo 1/3 subharmonic resonance. Comparing Figure 8a,c, it is found that the frequency range of the singlefrequency multi-amplitude occurring in the primary resonance is roughly the same as that of the 1/3 subharmonic resonance. In this frequency range, due to different initial conditions of the system, three different steady-state responses can be generated: large-amplitude vibration between wells, small-amplitude vibration in wells, and 1/3 subharmonic vibration.

Dynamic Analysis
(a) (b) (c) Figure 8. Displacement versus frequency with p = 9 m/s 2 , d = 24 mm: (a) nonlinear system (primary resonance); (b) linear system (primary resonance); (c) nonlinear system (1/3 subharmonic resonance). Figure 9 is the relationship between the average power and the excitation frequency of the nonlinear system, where Figure 9a is the primary resonance and can be obtained from Equation (45), Figure 9b is 1/3 subharmonic resonance, similar to the primary resonance, which can be obtained by deformation of Equation (51). In the primary resonance, the maximum power is 1.8 mW when the excitation frequency is 17.3 Hz, and the minimum power is less than 0.01 mW when the excitation frequency is 15~17.3 Hz. However, according to Huguet et al. [16], the inter-well motion is less robust in the range of 8.5 to  Figure 8c gives the amplitude-frequency response relationship of the nonlinear system when 1/3 subharmonic resonance occurs. It can be seen from the figure that only in a specific frequency range, i.e., 10.1~17.9 Hz, the system is likely to undergo 1/3 subharmonic resonance. Comparing Figure 8a,c, it is found that the frequency range of the single-frequency multi-amplitude occurring in the primary resonance is roughly the same as that of the 1/3 subharmonic resonance. In this frequency range, due to different initial conditions of the system, three different steady-state responses can be generated: large-amplitude vibration between wells, small-amplitude vibration in wells, and 1/3 subharmonic vibration. Figure 9 is the relationship between the average power and the excitation frequency of the nonlinear system, where Figure 9a is the primary resonance and can be obtained from Equation (45), Figure 9b is 1/3 subharmonic resonance, similar to the primary resonance, which can be obtained by deformation of Equation (51). In the primary resonance, the maximum power is 1.8 mW when the excitation frequency is 17.3 Hz, and the minimum power is less than 0.01 mW when the excitation frequency is 15~17.3 Hz. However, according to Huguet et al. [16], the inter-well motion is less robust in the range of 8.5 to 17.3 Hz, so the maximum power obtainable is 0.7 mW. When 1/3 subharmonic resonance occurs, the maximum power obtained is 0.25 mW, which is 36% of the maximum power of the primary resonance and exceeds 2500% of the minimum power of the primary resonance. Although the power of the subharmonic resonance is smaller than the power of the interwell motion of the primary resonance, it is much larger than the intra-well motion of the primary resonance. Therefore, the 1/3 subharmonic resonance expands the band of the nonlinear energy harvester by 92% so that the energy harvester can obtain more than 0.1 mW power in the frequency range of 18 Hz. Meanwhile, Figures 8 and 9 prove that the CDF method can predict the primary resonance and 1/3 subharmonic resonance behavior. 17.3 Hz, so the maximum power obtainable is 0.7 mW. When 1/3 subharmonic resonance occurs, the maximum power obtained is 0.25 mW, which is 36% of the maximum power of the primary resonance and exceeds 2500% of the minimum power of the primary resonance. Although the power of the subharmonic resonance is smaller than the power of the inter-well motion of the primary resonance, it is much larger than the intra-well motion of the primary resonance. Therefore, the 1/3 subharmonic resonance expands the band of the nonlinear energy harvester by 92% so that the energy harvester can obtain more than 0.1 mW power in the frequency range of 18 Hz. Meanwhile, Figures 8 and 9 prove that the CDF method can predict the primary resonance and 1/3 subharmonic resonance behavior.
(a) (b) Figure 9. Average power versus frequency with p = 9 m/s 2 , d = 24 mm, RL = 1000 kΩ: (a) nonlinear system (primary resonance); (b) nonlinear system (1/3 subharmonic resonance). Figure 10 shows the spectrum of the system at different excitation frequencies when the 1/3 subharmonic resonance occurs. It can be found that the response contains two frequency components, indicating that when 1/3 subharmonic resonance occurs, the primary harmonic and the 1/3 subharmonic coexist. When the excitation frequency increases, the primary harmonic component of the response becomes smaller and smaller, and the 1/3 subharmonic component becomes larger and larger as a whole. Combined with Figure  8c, the same conclusion can be drawn by observing that the primary harmonic coefficient a increases and the 1/3 subharmonic coefficient b decreases as the excitation frequency increases in Equation (55).   Figure 10 shows the spectrum of the system at different excitation frequencies when the 1/3 subharmonic resonance occurs. It can be found that the response contains two frequency components, indicating that when 1/3 subharmonic resonance occurs, the primary harmonic and the 1/3 subharmonic coexist. When the excitation frequency increases, the primary harmonic component of the response becomes smaller and smaller, and the 1/3 subharmonic component becomes larger and larger as a whole. Combined with Figure 8c, the same conclusion can be drawn by observing that the primary harmonic coefficient a increases and the 1/3 subharmonic coefficient b decreases as the excitation frequency increases in Equation (55). occurs, the maximum power obtained is 0.25 mW, which is 36% of the maximum power of the primary resonance and exceeds 2500% of the minimum power of the primary resonance. Although the power of the subharmonic resonance is smaller than the power of the inter-well motion of the primary resonance, it is much larger than the intra-well motion of the primary resonance. Therefore, the 1/3 subharmonic resonance expands the band of the nonlinear energy harvester by 92% so that the energy harvester can obtain more than 0.1 mW power in the frequency range of 18 Hz. Meanwhile, Figures 8 and 9 prove that the CDF method can predict the primary resonance and 1/3 subharmonic resonance behavior.
(a) (b) Figure 9. Average power versus frequency with p = 9 m/s 2 , d = 24 mm, RL = 1000 kΩ: (a) nonlinear system (primary resonance); (b) nonlinear system (1/3 subharmonic resonance). Figure 10 shows the spectrum of the system at different excitation frequencies when the 1/3 subharmonic resonance occurs. It can be found that the response contains two frequency components, indicating that when 1/3 subharmonic resonance occurs, the primary harmonic and the 1/3 subharmonic coexist. When the excitation frequency increases, the primary harmonic component of the response becomes smaller and smaller, and the 1/3 subharmonic component becomes larger and larger as a whole. Combined with Figure  8c, the same conclusion can be drawn by observing that the primary harmonic coefficient a increases and the 1/3 subharmonic coefficient b decreases as the excitation frequency increases in Equation (55).    Figure 11 illustrates the components of the time-domain diagram when the system undergoes 1/3 subharmonic resonance. Figure 11c,f shows the time-domain diagrams of the system's displacement when 1/3 subharmonic resonance occurs at the excitation frequency f = 14 Hz and f = 17.5 Hz, respectively. When the system undergoes 1/3 subharmonic resonance, the response contains the primary harmonic and 1/3 harmonic. Figure 11c is superimposed from the results of Figure 11a,b, where Figure 11a is the primary harmonic and Figure 11b is the 1/3 harmonic. Since the period of the 1/3 harmonic is three times that of the primary harmonic, the superposition results in the same period as the 1/3 harmonic period. Similarly, Figure 11f is a superposition of the results of Figure 11d,e. Comparing Figure 11a,b,d,e, it can be seen that the higher the excitation frequency, the smaller the primary harmonic component and the larger the 1/3 harmonic component. The primary harmonic can be regarded as a small perturbation applied to the 1/3 harmonic. The smaller the perturbation, the more stable the superposition result. the system's displacement when 1/3 subharmonic resonance occurs at the excitation frequency f = 14 Hz and f = 17.5 Hz, respectively. When the system undergoes 1/3 subharmonic resonance, the response contains the primary harmonic and 1/3 harmonic. Figure  11c is superimposed from the results of Figure 11a,b, where Figure 11a is the primary harmonic and Figure 11b is the 1/3 harmonic. Since the period of the 1/3 harmonic is three times that of the primary harmonic, the superposition results in the same period as the 1/3 harmonic period. Similarly, Figure 11f is a superposition of the results of Figure 11d,e. Comparing Figure 11a,b,d,e, it can be seen that the higher the excitation frequency, the smaller the primary harmonic component and the larger the 1/3 harmonic component.
The primary harmonic can be regarded as a small perturbation applied to the 1/3 harmonic. The smaller the perturbation, the more stable the superposition result. The conditions for 1/3 subharmonic resonance of the system are discussed below. Equation (51) is rewritten to It can be seen from Equation (57) that the conditions for the 1/3 subharmonic resonance of the nonlinear system are related to excitation frequency, excitation amplitude, The conditions for 1/3 subharmonic resonance of the system are discussed below. Equation (51) is rewritten to where a = a 2 . According to Weda' s theorem, the condition for the existence of positive real roots is It can be seen from Equation (57) that the conditions for the 1/3 subharmonic resonance of the nonlinear system are related to excitation frequency, excitation amplitude, linear stiffness coefficient, and nonlinear stiffness coefficient. Figure 12 is derived from Equation (57) and shows the range of the excitation amplitudes and frequencies that can occur for the 1/3 subharmonic resonance with other parameters held constant. When p < 1.5 m/s 2 , the 1/3 subharmonic resonance does not occur at any excitation frequency; when f < 5.5 Hz, the 1/3 subharmonic resonance does not occur at any excitation amplitude. In the region where the 1/3 subharmonic resonance can occur in the system, the corresponding excitation frequency range increases with excitation amplitude. Similarly, as the excitation frequency increases, the corresponding excitation amplitude range increases. Therefore, for a piezoelectric energy harvester whose parameters have been determined, the method can be used to grasp the external excitation conditions for its generation of the 1/3 subharmonic resonance; if a piezoelectric energy harvester is to be designed to broaden the bandwidth using the 1/3 harmonic resonance under specific environmental conditions, the method can be used to invert the parameters of the harvester and provide theoretical guidance for the design of the harvester.
linear stiffness coefficient, and nonlinear stiffness coefficient. Figure 12 is Equation (57) and shows the range of the excitation amplitudes and freque occur for the 1/3 subharmonic resonance with other parameters held const 1.5 m/s 2 , the 1/3 subharmonic resonance does not occur at any excitation fre f < 5.5 Hz, the 1/3 subharmonic resonance does not occur at any excitation the region where the 1/3 subharmonic resonance can occur in the system, th ing excitation frequency range increases with excitation amplitude. Similar tation frequency increases, the corresponding excitation amplitude ran Therefore, for a piezoelectric energy harvester whose parameters have bee the method can be used to grasp the external excitation conditions for its gen 1/3 subharmonic resonance; if a piezoelectric energy harvester is to be broaden the bandwidth using the 1/3 harmonic resonance under specific e conditions, the method can be used to invert the parameters of the harveste theoretical guidance for the design of the harvester.  Figure 13 is the diagram of the relationship between the response amplitude and the excitation amplitude of the system (if not specified, i case of the primary resonance). When the excitation amplitude is less tha system has only a small vibration; when the excitation amplitude is 5.9~36. tem has a single excitation amplitude corresponding to multiple response am unstable solution, and two stable solutions; and when the excitation amplit than 36.7m/s 2 , the system can maintain large vibration. The trajectories of re tudes during forward and reverse sweep are shown by arrows ③ and ④ When sweeping, the steady-state response under the previous excitation am initial condition for the vibration at the next excitation amplitude. It can be e the steady-state motion of the system depends on the initial condition in th excitation amplitude region. This is similar to Figure 8a, where there is a sweeping. For the range of excitation amplitudes that can produce multi-v not guarantee that the system always vibrates with large amplitudes. If th forms small amplitude vibration, the efficiency of the energy harvester is gr  Figure 13 is the diagram of the relationship between the response displacement amplitude and the excitation amplitude of the system (if not specified, it refers to the case of the primary resonance). When the excitation amplitude is less than 5.9 m/s 2 , the system has only a small vibration; when the excitation amplitude is 5.9~36.7 m/s 2 , the system has a single excitation amplitude corresponding to multiple response amplitudes, one unstable solution, and two stable solutions; and when the excitation amplitude is greater than 36.7 m/s 2 , the system can maintain large vibration. The trajectories of response amplitudes during forward and reverse sweep are shown by arrows 3 and 4 , respectively. When sweeping, the steady-state response under the previous excitation amplitude is the initial condition for the vibration at the next excitation amplitude. It can be easily seen that the steady-state motion of the system depends on the initial condition in the multi-value excitation amplitude region. This is similar to Figure 8a, where there is a jump during sweeping. For the range of excitation amplitudes that can produce multi-values, we cannot guarantee that the system always vibrates with large amplitudes. If the system performs small amplitude vibration, the efficiency of the energy harvester is greatly reduced. We are interested in considering the effect of 1/3 subharmonic motion.
icromachines 2021, 12, x FOR PEER REVIEW Figure 13. Displacement amplitude versus frequency with f=14 Hz, d=24 mm. Figure 14 shows the relationship between response amplitude and e plitude under different excitation frequencies obtained using the CDF m the excitation amplitude is in the range of 0 to 10 m/s 2 and the excitation 4, 10, and 16 Hz, the excitation amplitude range of the system with larg 0.7~10, 3.2~10, and 7.7~10 m/s 2 , respectively; the range of excitation am can ensure the system vibrates greatly is 1.5~10 m/s 2 , 0, and 0 respectiv quency excitation can obtain a large amplitude response, thus obtaining put power. However, the excitation amplitude range of multi-values is w small vibration in the multi-values may be obtained, so as to obtain the s At this time, the low-frequency excitation can obtain a larger response a this limits the bandwidth of the system. If the 1/3 subharmonic motion when the excitation frequencies are 4, 10, and 16 Hz, the excitation amp of 1/3 subharmonic motion in the system are 0, 3.2~8.7, and 7.2~32.7 m/ tion amplitude range of 1/3 subharmonic motion is basically in the ra valued excitation amplitude generated by primary resonance, and it enla increase in excitation frequency. Although the amplitude of 1/3 subha nance is smaller than that of the large amplitude motion of the primary is larger than that of the small amplitude of the primary resonance. If 1/3 motion is considered instead of the small amplitude motion of the prima the harvester can obtain more energy. For f = 4 Hz, the system has no resonance, consistent with Figure 12. Combined with the previous analy 1.5 m/s 2 , f > 5.5 Hz, 1/3 subharmonic resonance occurs in the system, a larger, the frequency range of 1/3 subharmonic resonance is wider. At t 1/3 subharmonic resonance can replace intra-well motion, so that the obtain higher energy output in the broader frequency range. In practica  Figure 14 shows the relationship between response amplitude and excitation amplitude under different excitation frequencies obtained using the CDF method. When the excitation amplitude is in the range of 0 to 10 m/s 2 and the excitation frequency is 4, 10, and 16 Hz, the excitation amplitude range of the system with large vibration is 0.7~10, 3.2~10, and 7.7~10 m/s 2 , respectively; the range of excitation amplitude that can ensure the system vibrates greatly is 1.5~10 m/s 2 , 0, and 0 respectively. High-frequency excitation can obtain a large amplitude response, thus obtaining a higher output power. However, the excitation amplitude range of multi-values is wide, and only small vibration in the multi-values may be obtained, so as to obtain the smaller output. At this time, the low-frequency excitation can obtain a larger response amplitude, but this limits the bandwidth of the system. If the 1/3 subharmonic motion is considered, when the excitation frequencies are 4, 10, and 16 Hz, the excitation amplitude ranges of 1/3 subharmonic motion in the system are 0, 3.2~8.7, and 7.2~32.7 m/s 2 . The excitation amplitude range of 1/3 subharmonic motion is basically in the range of multi-valued excitation amplitude generated by primary resonance, and it enlarges with the increase in excitation frequency. Although the amplitude of 1/3 subharmonic resonance is smaller than that of the large amplitude motion of the primary resonance, it is larger than that of the small amplitude of the primary resonance. If 1/3 subharmonic motion is considered instead of the small amplitude motion of the primary resonance, the harvester can obtain more energy. For f = 4 Hz, the system has no 1/3 harmonic resonance, consistent with Figure 12. Combined with the previous analysis, when p > 1.5 m/s 2 , f > 5.5 Hz, 1/3 subharmonic resonance occurs in the system, and when p is larger, the frequency range of 1/3 subharmonic resonance is wider. At this point, the 1/3 subharmonic resonance can replace intra-well motion, so that the harvester can obtain higher energy output in the broader frequency range. In practical application, the CDF method can provide theoretical guidance for the parameter design of the harvester to obtain higher power output in the environment of low-frequency excitation and small excitation amplitude.

Conclusions
In this paper, the influence of nonlinear characteristics on the ba piezoelectric vibration energy harvester is studied using the CDF metho modified the model, considered the influence of torque on the bending o ver beam when modeling with the Hamiltonian principle, and compare situation without considering the torque. It is found that torque greatly i equilibrium point and piezoelectric output of the harvester. In the ran parameters, with the increase in permanent magnet spacing, the influe also increases. Therefore, the influence of torque should be considered i ical modeling. Secondly, the static analysis of the system is carried out static bifurcation behavior and determine the parameter conditions fo state of the harvester. Then, the approximate analytical solution and am quency relationship of the primary resonance and 1/3 subharmonic res system are obtained using the CDF method. Finally, the CDF method is lyze the dynamic behavior of the energy harvester. The introduction of makes the harvester produce a bistable state and increases the bandwid resonance. When the 1/3 subharmonic vibration occurs in the harvester, i tains the primary harmonic and the 1/3 subharmonic, and with the in excitation frequency, the larger the proportion of the subharmonic. At th it is found that the frequency range of 1/3 subharmonic vibration is w quency range of multi-values generated by the primary resonance, and monic motion's amplitude is greater than that of the intra-well motion. subharmonic resonance also increases the excitation amplitude's band harvester at a specific excitation frequency. Therefore, to a certain exten monic motion can compensate for the low energy collection efficiency c

Conclusions
In this paper, the influence of nonlinear characteristics on the bandwidth of a piezoelectric vibration energy harvester is studied using the CDF method. Firstly, we modified the model, considered the influence of torque on the bending of the cantilever beam when modeling with the Hamiltonian principle, and compared it with the situation without considering the torque. It is found that torque greatly influences the equilibrium point and piezoelectric output of the harvester. In the range of bistable parameters, with the increase in permanent magnet spacing, the influence of torque also increases. Therefore, the influence of torque should be considered in the theoretical modeling. Secondly, the static analysis of the system is carried out to obtain the static bifurcation behavior and determine the parameter conditions for the bistable state of the harvester. Then, the approximate analytical solution and amplitude-frequency relationship of the primary resonance and 1/3 subharmonic resonance of the system are obtained using the CDF method. Finally, the CDF method is used to analyze the dynamic behavior of the energy harvester. The introduction of nonlinearity makes the harvester produce a bistable state and increases the bandwidth of primary resonance. When the 1/3 subharmonic vibration occurs in the harvester, it mainly contains the primary harmonic and the 1/3 subharmonic, and with the increase in the excitation frequency, the larger the proportion of the subharmonic. At the same time, it is found that the frequency range of 1/3 subharmonic vibration is within the frequency range of multi-values generated by the primary resonance, and the subharmonic motion's amplitude is greater than that of the intra-well motion. Similarly, 1/3 subharmonic resonance also increases the excitation amplitude's bandwidth of the harvester at a specific excitation frequency. Therefore, to a certain extent, the subharmonic motion can compensate for the low energy collection efficiency caused by the low robustness of inter-well motion. Making full use of the nonlinear characteristics of a system can effectively broaden the bandwidth of the piezoelectric vibration energy harvester and improve the efficiency of obtaining environmental energy. The theoretical method (CDF) can accurately predict the dynamic behavior of this kind of strongly nonlinear vibration energy harvester and provide guidance for its theoretical research.
In the next step, the parameters of the harvester model are optimized by considering the effects of residual stress, temperature, damping, and the dimensions of the multilayer beam. On this basis, the following studies will be carried out. 1 An optimal load resistance exists in the system. A more significant power can be output under the same conditions when the optimal load resistance is obtained. However, the variation of the load resistance may affect the dynamic behavior of the system. Determining the optimal load by integrating the mechanical system and the circuit system remains to be studied. 2 Subharmonic resonance occurs only within a specific range related to the excitation frequency, excitation amplitude, linear stiffness coefficient, and nonlinear stiffness coefficient. By adjusting these parameters, the subharmonic resonance can maximize the bandwidth of the harvester.
Author Contributions: Z.Q. and Y.Y. designed the model and research idea; W.W. and Y.Y. contributed theoretical analysis; Z.Q. and Y.Y. wrote the paper. All authors have read and agreed to the published version of the manuscript.