Performance Analysis of an Electromagnetically Coupled Piezoelectric Energy Scavenger

: The deliberate introduction of nonlinearities is widely used as an effective technique for the bandwidth broadening of conventional linear energy harvesting devices. This approach not only results in a more uniform behavior of the output power within a wider frequency band through bending the resonance response, but also contributes to energy harvesting from low ‐ frequency excitations by activation of superharmonic resonances. This article investigates the nonlinear dynamics of a monostable piezoelectric harvester under a self ‐ powered electromagnetic actuation. To this end, the governing nonlinear partial differential equations of the proposed harvester are order ‐ reduced and solved by means of the perturbation method of multiple scales. The results indicate that, according to the excitation amplitude and load resistance, different responses can be distinguished at the primary resonance. The system behavior may involve the traditional bending of response curves, Hopf bifurcations, and instability regions. Furthermore, an order ‐ two superharmonic resonance is observed, which is activated at lower excitations in comparison to order ‐ three conventional resonances of the Duffing ‐ type resonator. This secondary resonance makes it possible to extract considerable amounts of power at fractions of natural frequency, which is very beneficial in micro ‐ electro ‐ mechanical systems (MEMS) ‐ based harvesters with generally high resonance frequencies. The extracted power in both primary and superharmonic resonances are analytically calculated, then verified by a numerical solution where a good agreement is observed between the results.


Introduction
Recent major progress has been made in Micro-Electro-Mechanical Systems (MEMS) technology, which has led to the development of various sensors and actuators on the one hand, and a vast reduction in the size and power consumption of Complementary Metal-Oxide-Semiconductor (CMOS) circuitry on the other hand, leading to the idea of developing small-scale harvesters to provide the required energy of intelligent systems and supersede former heavy batteries with a limited lifetime. Utilization of these harvesters in combination with conventional systems leads to energy autonomous devices with the capability of self-energy generation. These autonomous devices can be exploited in many applications, such as Wireless Sensor Networks (WSNs) [1], Implantable Medical Devices (IMDs) [2], and active Radio Frequency Identification systems (active-RFIDs) [3].
Among the various sources of environmental energy, such as thermal, electromagnetic and solar energies, the kinetic energy stored in mechanical vibrations has been shown to be a promising candidate for energy harvesting, considering the two important factors of availability and efficiency [4]. Vibrational energy harvesters consist of an either linear or nonlinear mechanical resonator composed of a mechanical spring as their main part, along with some lumped or distributed mass and an energy conversion mechanism which damps the mechanical vibration and transduces the energy to the electrical domain. Some energy harvesters exploit the electrostatic force between opposite charges to the convert kinetic type energy to an electrical type one, stored in an electrical field [5]. This occurs while the others employ Lorentz force to induce a voltage based on Faraday's law in a winding coil and convert the mechanical energy to an electrical one through the generation of time-varying electromagnetic fields [6]. However, the most prevalent method for energy conversion on a small scale is using electroelastic and smart materials. Piezoelectric [7,8] and triboelectric [9,10] materials are among these types of materials, whose application in energy harvesting is widely studied by the researchers in this area. The straightforward fabrication in a small scale, along with considerable amounts of harvested power, are the main advantages of using these materials in kinetic energy harvesters. Small-scale harvesters generally employ lumped masses connected to a mechanical beam to realize the mechanical resonator. Employed beams are used in cantilevered or doubly clamped boundary conditions. In electromagnetic and electrostatic harvesters, the beam plays the role of a spring, which suspends an electro or permanent magnet [11] or the electrodes of a capacitor [12], respectively. In other words, in a piezoelectric harvester the beam also performs as a base for the piezoelectric material, generating tension and compression in the piezoelectric material when bended [13].
While in a linear harvester the efficiency of the device can be enhanced by limiting the mechanical damping effects and increasing the quality factor of the system, these high efficiency harvesters usually have very narrow bandwidths, which restricts their deployment in real applications [14][15][16][17][18]. Several bandwidth broadening strategies have been adopted to circumvent this so-called "gain-bandwidth dilemma" [19], and employing nonlinear effects is one of them. Adding nonlinearities to a resonant harvester can lead to bandwidth improvement in response to singleharmonic excitations as a result of bending of the frequency response curves [20][21][22][23][24]. Additionally, when subjected to random excitations, a bi-stable nonlinear harvester can generate considerable amounts of power through the activation of noise-assisted jumps between its energy wells [25][26][27]. The nonlinear responses of a resonant harvester, due to the connection to nonlinear circuits of energy harvesting [28] and the nonlinear behavior of the piezoelectric material [29], have been investigated previously. The effect of geometric nonlinearities due to mid-plane stretching, nonlinear curvature, and nonlinear inertia in energy harvesting beams with doubly clamped [30] and cantilevered [31] configurations have also been studied.
The nonlinear behavior of piezoelectric harvesters may arise from geometrical nonlinearities, the nonlinear response of the piezoelectric ceramics, or coupling with nonlinear energy harvesting circuits. It has been shown that the behavior of a harvester under the effect of any type of these nonlinearities may be modeled using a Duffing-type governing equation [32][33][34]. It should be noted that a Duffing oscillator used as an energy scavenger suffers from a major drawback due to the coexistence of low-and high-energy stable equilibria, making its wideband operation dependent on the trajectory of the excitation [21]. As a result, the enhancement of the useful bandwidth only occurs in the case that the response converges to the high-energy output of the system [35]. Due to this fact, several previous research projects have been carried out on possible procedures forcing the harvester output to lie on the high-energy branch of the frequency response curve [36][37][38][39].
In addition to naturally occurring materials possessing piezoelectric properties, artificial materials such as metamaterials have been vastly used in energy harvesters in recent years. Metamaterials are artificial materials engineered to have properties that do not exist in nature. These properties arise from the structure and not the composition of the material itself, and have led to the evolution of a new generation of micro-scale devices [40,41]. The application of these types of materials in energy harvesting from electromagnetic Radio Frequency (RF) waves has been investigated by many researchers in this area [42,43].
As evidenced by the literature, all of the above-mentioned sources of nonlinearity lead to third order nonlinear effects. As a result, while they are able to extract energy from the third harmonic of excitation through the activation of a superharmonic resonance, they remain ineffective when used in environments in which the second harmonic of the excitation dominates. In this paper, the nonlinear behavior of a piezoelectric harvester composed of a cantilevered unimorph beam with an attached tip magnet under an electromagnetic force at both primary-and second-order superharmonic resonances is investigated. Nonlinear, electromechanically coupled, partial differential equations of the system are presented and order-reduced by means of Galerkin's method. For excitations in the neighborhood of the first-mode resonance frequency, the governing equations are solved by implementing the perturbation method of multiple scales. The system response and harvested power are analyzed in both primary and superharmonic resonances.

Harvester Model and Governing Equations
A micro harvester composed of a cantilevered silicon microbeam of length L , width W and thickness B h is shown in Figure 1  As shown in Figure 2, energy extraction is performed by the placement of two electrodes, which are very thin layers of conductive material on the both sides of the piezoelectric layer. The effect of the electrodes' stiffness is negligible on the mechanical properties of the system. The potential difference between these electrodes is shown by V in the figure, which is actually the output voltage of the harvester. P in this figure denotes the initial polarization vector of the ferroelectric-type piezoelectric ceramic, and shows that the piezoelectric layer is initially polarized along the thickness. This voltage is applied to two circuits. The first one is an amplifier circuit supplying the electromagnet, which may contain a simple step-up transformer or an op-amp. The second circuit is a switching-type AC-DC converter which supplies the target load, which may be a sensor node-for example, in wireless sensor network (WSN) applications. The total effective input impedance of these two circuits, which is observed by the harvester, will be substituted by a pure resistance R in the calculations. To increase the total area of the piezoelectric layer and consequently enhance the amount of harvested power, wide beams are normally used in the structure of these types of harvesters. Therefore, the polar moment of inertia of the beam cross-section is much higher than its area moment of inertia, and consequently the effect of torsion may be neglected in comparison to the bending effect within a deformation study. According to this fact, only the y -component of the magnetic force, which leads to lateral bending of the structure, will be considered in the present study. Due to the assumed symmetry, the y -component of the magnetic force has to be an odd function of the beam deflection at its free end and, therefore, one can write where  denotes the displacement of the beam free end in y direction, and the coefficients i b depend on the characteristics of the amplifier circuit, the coil used in the electromagnet, and the permanent magnet attached to the beam. An equivalent circuit model of the electromechanically coupled piezoelectric generator presented in Figures 1 and 2 is shown in Figure 3. As can be seen in the figure, the nonlinear impedance in the left-hand side of the circuit is composed of three terms, where the ones proportional to 1  and 4  represent the effect of nonlinear curvature and nonlinear inertia due to the large amplitude vibration of the beam [31], and 3  represents the effect of a nonlinear magnetic force applied to the permanent magnet. Since the beam is assumed to be a wide one, its polar moment of inertia is much higher than a bending one and, therefore, the torsion effect may be neglected. It should be noted that all the nonlinear terms are brought into account up to the third order. According to this figure, the coupled electromechanical equations governing the dynamics of the system can be written as Introducing the following parameters, Equation (2) can be simplified as where hat superscripts are omitted for the sake of simplicity.

Primary Resonance
In this section, the nonlinear system of governing equations of the harvester is solved by employing the perturbations method of multiple scale. Similar to any other perturbative technique, this analysis may be applied to nonlinear systems perturbed from a linear one by a small value. Therefore, the displacements in the structure must remain within a short interval around the equilibrium position. To this end, a small parameter  is introduced and the vibration amplitude is supposed to be comparable with this parameter. Furthermore, since the electromechanical coupling of piezoelectric structures is normally weak, the coupling parameters 2  and 5  are assumed to be of order  , too. In order that the nonlinear term, due to electromagnetic actuation, is comparable with nonlinear terms due to the large curvature of the beam, 3  is assumed to be of order 1 Therefore, one can let By substitution of these relations into Equation (4), one obtains: where the following assumptions are made to guarantee that damping and external forces appear in the same order as the nonlinearity term [44]: New independent variables, called time scales, may be introduced accordingly as Representing the operator denoting partial differentiation with respect to n T by n D (i.e.
The solution can be represented by expanding u and  as follows Since the excitation is introducing the detuning parameter  as Substituting Equations (11) and (12) into Equation (7), then equating terms of the same powers of  , the following sets of equations are obtained The solution of Equation (13) can be written as Substituting the above solution into Equation (14a), one obtains where cc denotes the complex conjugate terms. The secular terms in Equation (17) will be eliminated if Substituting Equations (16) and (19) into (14b), the solution is obtained as where NST stands for non-secular terms with frequency 3 E  . Therefore, eliminating the secular terms results in Representing A in polar form: Substituting into Equation (22), then separating real and imaginary parts, one obtains To transform the equations to an autonomous system, the following change of variable is used in the above equation: The steady state solution of the dynamical system described in Equation (24) According to the above equations, the linear stability of equilibrium points and location of bifurcation points can be found through eigenvalues of the following Jacobian matrix    (27) in which the equilibrium point of the system is shown by   Therefore, one can expand the steady state solution of Equation (24) as a Fourier series, as follows where  is the angular frequency of the periodic solution. When the solution is a fixed point and no limit cycle exists, the only nonzero coefficient in Equation (28) will be The coefficient i c may be derived by performing a Fast Fourier Transform (FFT) on the temporal solution of the dynamical system represented by Equation (24), found by a direct numerical integration. Substituting Equations (28) and (25) Averaging the simultaneous output power of the harvester over one period of vibration results in the following average power: It should be noted that the above solution is valid while   goes beyond  and Equation (34), which sums the power of all harmonics assuming they have different frequencies, will not be accurate anymore. Fortunately, since a good convergence will be obtained by values of N far below this limit, Equation (34) is applicable with a very good accuracy.
For the case in which no limit cycles are involved, Equations (31) and (34)

Superharmonic Resonance
In this section, we analyze the superharmonic resonance of the system generated by the secondorder nonlinearity. Because this resonance is activated at lower excitation levels in comparison to conventional order-three superharmonic resonance of the Duffing-type nonlinear beam equation, it can more effectively contribute to energy generation. Assuming K k   and    , then substituting together with (5) and (6) into (4), the governing equations will be of the following form: The equations of order The solution of Equation (38) may be written as Substituting this solution into Equation (39a), one obtains the condition of secular terms elimination is where the steady sate solution as t   is Substituting the above solution into Equation (40b), the output voltage is obtained as 2  2  5  3 5  1  2  2  2  2  2  2  2  2  2  2   2   2  2  2  2   cos  1  2  4 2  1  1 cos 2 Therefore, using Equation (33)

Results and Discussion
In this section, the nonlinearities effect on the dynamic response of the system and amount of power extracted is explored using derived relations in the case of an MEMS-based harvester. The harvester has a double-layered cantilever beam of width b and the mechanical quality factor of the system are assumed to be 0.25 and 300 , respectively. Table 1 presents the required material properties of silicon and PZT-5H (which is poled along the thickness) [17].
The frequency response of the system at an excitation level defined by a base acceleration Figure 4, where the external load is 2 k R  . As can be seen in this figure, the response curve is composed of two parts: the first is a simple, right-bended curve which traditionally appears in the response of Duffing-type oscillators with nonlinearities of the hardening type. Two stable high-and low-energy branches, connected through a curve with saddle stability, which is the locus of system saddle points, may be distinguished in this part of the response. Two points marked by SN1 and SN2 define saddle-node bifurcation points of the system, in which the saddle point collides with a stable node. The second part of the response involves one unstable branch and another branch with saddle stability. Points SN3 and SN4 define bifurcation points due to the incidence of a saddle point with an unstable node. Since the second part contains no stable nodes, it cannot contribute to power generation, and plays no important role in the harvester's behavior.
Increasing the load resistance leads to the two parts of the response approaching each other and, at a critical load, points SN2 and SN3 collide, resulting in a Hopf bifurcation. Figure 5 illustrates the system response at the same excitation level as the previous case, where the load resistance is increased to 3 k R  . Based on the calculations, the first Lyapunov coefficient for the dynamical system, represented by Equation (24), is found to be negative for the Hopf bifurcation point, indicating that it is of the supercritical type. Therefore, stable limit cycles are formed with an unstable focus in their interior. As the excitation frequency increases after the Hopf bifurcation occurrence, these limit cycles grow, and at point X1 they collide with a saddle point at X2 and a homoclinic bifurcation occurs, which leads to the response falling from the higher energy branch to the lower energy one at point X3.

k R  ).
Further increasing the load resistance pushes the Hopf bifurcation point toward lower frequencies. Figure 7 depicts the response of the system at load 5 k R  . As illustrated in this figure, stable limit cycles exist between the Hopf bifurcation point and a very close adjacent point X1, and at this point the limit cycles lose their stability, and therefore an instability region is formed between X1 and the bifurcation point SN1, where a stable equilibrium point is created. In this region, since there is no stable solution, the system response grows unbounded up to this order of approximation. However, in reality, the motion does not grow without bound. As the amplitude further increases, higher orders of nonlinearity play an important role and in the present analysis, accounting for nonlinearities up to the third order is not valid anymore. Furthermore, in addition to taking higher order terms into account, a method capable of treating strongly nonlinear problems must be employed. However, for the current case study this region must be avoided since a greater increase in the deflection amplitude causes the beam to damage and fracture. Further increasing the load resistance causes the Hop bifurcation point to move again toward high frequencies and disappear at a specific load.   -10 show the amount of extracted power with respect to the excitation frequency for the above-mentioned three cases. As expected, the curve representing the amount of output power of the system behaves similarly to the vibration amplitude curve and is a multivalued one too, so the correct value of output power depends on the trajectory of the excitation. Ramping up the frequency from low excitation frequencies pushes the harvester along the high-energy stable branch of the response curve, while ramping it down causes the harvester power to lie on the low-energy branch. To examine the validity of our analytical solution, the equations are solved numerically by employing the Runge-Kutta method. To achieve this goal, the solution is numerically calculated in both a forward (increasing frequency) and backward (decreasing frequency) sweep. The comparison is presented for all three cases. As can be seen in these figures, a very good agreement is observed between the results of both the analytical and numerical methods.
In Figure 8, it can be seen that sweeping the excitation frequency from low to high frequencies first leads to enhancement of the extracted power. At the bifurcation point corresponding to SN2 in Figure 5, a jump takes place and power falls to the lower energy branch. Reducing the excitation frequency first results in extracted power enhancement, at the SN1 point power jumps to the higher energy branch, and after that a further decrease in frequency leads to power reduction.   As illustrated in Figure 9, a significant difference exists in the system response in this figure compared to the response shown in Figure 8. In this case, increasing the excitation frequency increases the output power before the Hopf bifurcation takes place. After that, the creation of limit cycles leads to a decrease in extracted power. At a specific point corresponding to X1 in Figure 5, due to the collision of the limit cycles with the saddle point and occurrence of a homoclinic bifurcation, the power falls down to the lower branch. As shown in Figures 8 and 9, one of the most significant advantages of nonlinear harvesters in comparison to their linear counterparts is their improved bandwidths. While typical linear harvesters have extremely narrow bandwidths, the bandwidth of a nonlinear harvester is enhanced due to the bending of the response curves, which leads to the generation of considerable amounts of energy in wider excitation frequency intervals. Figure 11 shows the temporal response of output voltage on the high-energy branch for two excitation frequencies, one before the Hopf bifurcation point and one after but before the homoclinic bifurcation occurrence. As can be seen in this figure, the creation of limit cycles leads to the enhancement of the peak output voltage, which can be beneficial in many harvesters but especially MEMS-based types, which usually generate low output voltages due to the thinness of their piezoelectric layer. This is because the maximum piezoelectric thickness which can be deposited is generally limited to a few m  in conventional deposition processes [45].
Although at first it seems that the power should be enhanced due to the creation of new harmonics in the response, as mentioned in Figure 9, limiting cycles' formation leads to an output power decrease. The reason for this can be seen in Figure 12, which shows the frequency domain of the voltage response. As can be seen in this figure, while other harmonics are added to the response due to the formation of limit cycles, the amplitude of the main harmonic is decreased significantly, which leads to a decrease in total generated power. This is because as the response point   , a  , defining the amplitude and phase of the vibration, periodically moves on a limit cycle, it takes more time for the point to pass through the region near the saddle point and, since the vibration amplitude and consequently output voltage at the saddle point are less than the node at the center of the limit cycle, this leads to a considerable decrease in the main harmonic amplitude. Figure 10 shows the variation in extracted power with excitation frequency for the case of 5 k R  . As can be seen in this figure, as the excitation frequency is increased the output power increases before Hopf bifurcation takes place, and after that, similar to the previous case, the power is decreased due to the creation of limit cycles in a very narrow frequency band before the instability occurs. A further increase in the excitation frequency leads to the disappearance of stable limit cycles, and the system enters an instability region. In this region, as mentioned before, there is no stable point up to this order of approximation, and as the numerical solution confirms the multiple scale solution, the vibration amplitude grows unbounded up to this order of approximation. After the instability region finishes, any further increase in excitation frequency decreases the output power.   Figure 4. At a critical excitation, a Hopf bifurcation occurs at a critical load resistance and the system behavior changes to that shown in Figure 5. For the systems with lower mechanical quality factors, this Hopf bifurcation appears at higher excitations. Increasing the excitation levels leads to the appearance of this bifurcation in a wider interval of load resistances. In a sub-domain of this region, an instability interval starts to appear in the system response, as shown in Figure 7.
The variation in maximum extractable power at superharmonic resonance with respect to load resistance is illustrated in Figure 14. From Equation (47), the peak power occurs approximately at half the linear natural frequency E  . As mentioned in Equation (46), two harmonics are involved in the voltage response, the first with a frequency equal to the excitation frequency and the second with a frequency two times this. As this figure illustrates, the power due to the first harmonic is maximized at a lower load resistance in comparison to the second harmonic, and the optimum load opt R , at which the total extracted power reaches its maximum, falls between these two values. As can be seen in this figure, most of the power is generated by the second harmonic, which is created as the result of superharmonic resonance. Therefore, this secondary resonance, which makes it possible to extract considerable amounts of power at fractions of natural frequency, can be very beneficial, especially in MEMS-based devices which have high natural frequencies in comparison to typical environmental vibrations.  Consequently, it can be concluded that by increasing the excitation amplitude the effect of superharmonic resonance is amplified, and the contribution of the second harmonic to the generated power is increased. For intermediate excitation levels the optimum resistance should be selected with respect to the system mechanical quality factor.

Conclusions
Nonlinear dynamics of an electromagnetically coupled, piezoelectric unimorph beam, employed as a vibrational energy harvester during large amplitude motions, were studied. To investigate the harvester behavior in the neighborhood of first resonance frequency, the governing order-reduced nonlinear differential equations of the system are firstly presented. An analytical solution employing the perturbation method of multiple scales is presented for the equations at both primary and superharmonic resonances. The results indicate that, at relatively small excitation amplitudes, the system response involves a typical Duffing resonance, which is a common behavior of cantilever beams under large lateral vibrations due to the presence of curvature and inertia nonlinearities. As the excitation amplitude is increased, a supercritical Hopf bifurcation occurs on the high-energy branch of the response curve and stable limit cycles are created around the unstable nodes in a narrow frequency band. By increasing the frequency, these limit cycles grow and finally collide with the saddle equilibrium point of the system, leading to a homoclinic bifurcation which pushes the response to fall from the high-to the low-energy branch. These limit cycles lead to the enhancement of output voltage while decreasing the output power of the harvester. Increasing the load resistance at this excitation level moves this bifurcation point toward lower frequencies so that the created limit cycles will be unstable before a saddle point is created, and a narrow instability region is observed where there is no stable limit cycle or equilibrium point up to this order of approximation. Although, in reality, the response amplitude will be bounded considering higher orders of nonlinearity, the system operation at this region should be avoided, since it can cause the beam to fracture. Any further increase in the load resistance moves the Hopf bifurcation point toward higher frequencies again, and this region disappears. The analytical findings are verified using a numerical solution where the results are in a very good agreement. Due to electromagnetic coupling, a superharmonic resonance occurs at half the natural frequency, which is able to generate energy at smaller excitations compared to the ordinary one-third superharmonic resonance of the Duffing oscillator. This can be very beneficial, especially in MEMS-based devices which have high natural frequencies in comparison to the environmental vibration spectrum.