Bifurcation analysis and nonlinear dynamics of a capacitive energy harvester in the vicinity of the primary and secondary resonances

The objective of the present study is to examine the effect of nonlinearity on the efficiency enhancement of a capacitive energy harvester. The model consists of a cantilever microbeam underneath which there is an electret layer with a surface voltage, which is responsible for the driving energy. The packaged device is exposed to unwanted harmonic mechanical excitation. The microbeam undergoes mechanical vibration, and accordingly, the energy is harvested throughout the output electric circuit. The dynamic formulation accounts for nonlinear curvature, inertia, and nonlinear electrostatic force. The efficiency of the device in the vicinity of the primary and super-harmonic resonances is examined, and accordingly, the output power is evaluated. Bifurcation analysis is carried out on the dynamics of the system by detecting the bifurcations in the frequency domain and diagnosing their respective types. One of the challenging issues in the design and analysis of energy-harvesting devices is to broaden the bandwidth so that more frequencies are potentially accomodated within the amplification region. In this study, the effect of the nonlinearity on the bandwidth broadening, as well as efficiency improvement of the device, are examined. It is seen that as the base excitation amplitude increases, the vibration amplitude does also increase and accordingly the nonlinearity dominates. The super-harmonic resonance regions emerge and get bigger as the vibration amplitude increases, and pull-in gaps appear in the frequency response curves.


Introduction
Wireless sensor networks (WSNs) have been focused in recent years due to their numerous advantages such as ease of installation, reduction in cost and weight, and elimination of wire connections [1]. Providing the power supply source is a challenge issue in these systems. Batteries provide an impermanent power supply solution which it takes lots of time and cost to replace them. An interesting and long-term solution to extend the operational time of WSNs is vibration energy-harvesting devices which convert the ambient vibrations into electrical energy to charge the batteries or replacing them. In the past few years, many researchers have focused on the dynamics of energyharvesting devices so that to enhance the efficiency of the energy harvester [1][2][3]. Among the numerous publications reported so far, some of them have focused on harvesting energy based on electromagnetic effect [4][5][6][7][8][9], and some have focused on piezoelectric-based ones [1,2,[10][11][12][13][14]. Capacitive energyharvesting devices have also been focused in the literature thanks to their ease of fabrication; however, as a major challenge, they require a voltage source to induce electric charges on the electrodes of the variable capacitor and consequently start the conversion [3,[15][16][17][18]. The required voltage can be provided by means of using dielectric, capable of keeping internal charges for many years. In variable gap capacitors, the capacitance value can be changed by different physical quantities such as pressure or acceleration. This change is measured by an electronic circuitry. Due to the electrostatic nonlinearity, the working range of devices using electrostatic force is limited. Some researchers have explored the possibility of amplitude enhancement thanks to the internal resonance [14,19,20]. A very important research subject which has received a great research interest is the effect of nonlinearity on the performance and the efficiency of an energy-harvesting devices. The findings of the most important of these are reviewed here. Ngan Tran et al. [21] reported a review on nonlinear techniques for performance enhancement of ambient vibration of energy harvesters. They reported the most important publications devoted to the effect of stochastic loading, internal resonance, being multidegree of freedom, mechanical stoppers, and parametric excitation, which lead to nonlinear behavior and accordingly enhance the efficiency of the energyharvesting device. Lu et al. [22] explored the connection between the resonance response interaction and bubble-shaped response curve that can occur in the forced response of a nonlinear magnetoelectric coupled system. Hao Wu et al. [23] proposed a nonlinear two-degree-of-freedom piezoelectric energy harvester based on the similar linear model proposed by Wu et al. [24]; they achieved broader operating bandwidth in mono-stable condition by taking benefit of the nonlinearity and also tuning the resonance response peaks in the frequency domain. A Cammarano et al. [25] improved the bandwidth of an energy harvester exploiting nonlinearities. They reported that the stronger is the nonlinearity, the broader is the bandwidth and accordingly the better is the performance; this is due to the fact that more frequency contents can be accommodated in the bandwidth.
Malaji and Ali [26] investigated the energy harvesting from a two-degree-of-freedom near periodic structure composed of two pendulums connected using a linear spring. They studied the mutual effect of mistuning due to changing the length of pendulums and the amount of harvested energy on each other. A bistable energy harvester subject to periodic and quasi-periodic excitations has been examined by Rajagopal et al. [27]. They investigated the nonlinear behavior of such system by exploring different dynamical properties. Using the bifurcation diagrams of the harvester under periodic excitation, they showed that hidden attractors coexist. They also used a 2D lattice array of the harvesters to study their collective behavior. Kumar et al. [28] considered the energy harvested from the vibration of a flexible cantilevered flapper induced by vortex in the wake of a rigid circular cylinder. They carried out experimental analysis in the wind tunnel, in addition to numerical simulations, to investigate how the amount of harvested energy is affected by the gap between the cylinder and the cantilever. They showed that the location of vortex shedding with respect to the flapper will vary as the flow speed changes, and this will affect the dynamics of the flapper and consequently the energy harvesting.
Based on the literature, one of the main issues with the vibrating energy harvesters is the narrow bandwidth in case of high-quality factor resonators [21]; the impetus of the current research is to broaden the bandwidth of a capacitive energy-harvesting device exposed to unwanted harmonic external excitations, thanks to the existence of nonlinearity. We also want to determine the low-frequency resonances which are activated due to the nonlinearity of the response and mainly known as super-harmonic resonance regions. Since the mechanical excitations are mostly of low frequency, in order to harvest the energy of lowfrequency ambient mechanical vibrations, activation of super-harmonic resonances is a necessity; besides the endeavor made to broaden the bandwidth, we have examined the nonlinear resonance regions to maximize the energy harvesting in the super-harmonic resonance domains. The effect of nonlinearities including geometric, inertia, and nonlinear displacement-dependent electrostatic force on the efficiency of the device is taken into account. The study encompasses the response in the vicinity of both primary and super-harmonic resonances. The frequency response curves are computed by means of shooting technique, and the stability of the periodic solutions is ascertained based on Floquet theory. The bifurcations in the frequency domain are determined, and their types are examined based on the loci of the Floquet exponents with respect to the unit circle.

Model description and mathematical modeling
As depicted in Fig. 1, the proposed energy harvester is composed of a cantilever beam exposed to underneath capacitance. Between the lower electrode and the cantilever beam exists an electret layer with a surface voltage of V s , thickness of h e , and relative permittivity of e e , which is responsible for providing the required voltage of the device.
In order to model the mechanical behavior of the micro-cantilever, an Euler-Bernoulli beam theory is utilized. The coordinate system (x-y-z) is attached to the fixed end of the cantilever beam with its origin on the neutral axis and the x-axis along the beam and the y-axis in the direction of the base excitation. Assuming the energy harvester operating as a high-quality factor resonator and accordingly undergoing large deformations, the geometrical and inertial nonlinearities due to inextensionality conditions get more pronounced and the linear theory does not hold true anymore. Therefore, a full nonlinear beam model is employed which takes the curvature nonlinearities and shortening effects into account. The nonlinear equation of motion for the considered micro-cantilever based on the Euler-Bernoulli beam theory is given by [29], qA € w þ EI w 000 þ w 0 w 002 þ w 000 w 02 À Á 0 ÀJ € w 00 where w denotes the transverse displacement and the superscript ( 0 ) represents the derivative with respect to the s which is the undeformed length measured from the root of the beam to the reference point, q is the mass density, A is the cross-sectional area, and E is the modulus of elasticity. Here, I denotes the area moment of inertia with respect to the neutral axis, J is the rotary inertia, w b is the base displacement due to ambient vibrations, and F E is the electrostatic force per unit length acting on the micro-cantilever. The first and third terms on the left-hand side of Eq. (1) account for translational and rotary inertias, respectively. The first part of the second term in the left-hand side accounts for linear stiffness, and the last part is due to curvature nonlinearities. Also, the last integral term in left-hand side accounts for shortening effect due to inertial nonlinearities. The last term in right-hand side accounts for translational inertia due to base excitation and is considered as a harmonic displacement in the form of: The electrostatic force between the electrodes of a parallel plate capacitor is given by [30]: where U E is the electric potential energy stored in the capacitor and is given by [30]: where V s is the surface voltage of the electret layer, and C denotes the equivalent capacitance of the Fig. 1 Schematics of the proposed energy harvester system which consists of two series capacitances including the electret layer capacitance (C e ) and the variable capacitance (C 1 ). To obtain the equivalent capacitance of the device, the capacitor is considered as an infinite number of differential parallel capacitors, Each of them is in series with the electret layer capacitance.
Neglecting the fringing effects and assuming complete overlap between the micro-cantilever and the substrate, the equivalent differential capacitance of the system is obtained as follows: Integrating Eq. (7) over the length and substituting the resultant in Eq. (4), and using Eq. (3), the electrostatic force per unit length acting on the micro-cantilever is obtained as follows: Neglecting the effects of parasitic capacitances, the equivalent electric circuit of the energy harvester is modeled as Fig. 2.
Applying the Kirchhoff's voltage law in the equivalent circuit of the device, the electrical governing equation of the system reduces to [3]: Equation (1) along with Eq. (9) forms the nonlinear electromechanical governing equation of the system subjected to the following boundary conditions: For convenience, the governing equations of the system are non-dimensionalized using following dimensionless parameters: where Q e is the electrical charge stored in electret layer and is defined as: Introducing the dimensionless parameters given in Eq. (11) and the electrostatic force in Eq. (8) to Eq. (1) and Eq. (9) and removing the hat notation, the dimensionless governing equations reduce to: where and the boundary conditions in the non-dimensional form reduce to: The Galerkin decomposition method is employed to eliminate the spatial dependence in Eqs. (13) and (14). To this end, the transverse deflection of the micro-cantilever is represented as a series expansion in terms of the eigenfunctions of the micro-cantilever, i.e., where U i t ð Þ is the ith generalized coordinate and u i s ð Þ is the ith linear undamped mode shape of the microcantilever. Based on Galerkin method, the reducedorder model of the system is obtained as given in appendix. The mean output power of the energy harvester P is obtained based on the amplitude of the steady-state response as: Due to the surface voltage of the electret layer, it is necessary to calculate the static pull-in voltage of the micro-cantilever. To this end, the time-dependent terms in Eq. (13) are eliminated and the static behavior of the micro-cantilever is given by: Using the Galerkin method, the static deflection and the pull-in voltage of the energy harvester are determined.

Energy conservation
To clarify the problem and examine the energy conservation, we have developed the following equivalent electromechanical circuit for the reduced-order model which accounts for the output energy through the electret layer and the pumped energy throughout the base excitation. A part of this energy is harvested in the resistance, and the remaining is stored in terms of kinetic and potential energies in the vibration system and the potential energy in the variable gap capacitor.
Equation (20) represents the energy conservation where the left-hand side corresponds to input energies and the right-hand side represents the harvested, mechanical, and dissipated energies. where: In Eq. (20), W ex is the input energy due to the base excitation, QV s represents the output energy from the electret layer, E m is the mechanical energy of the vibratory system, and the last term in Eq. (20) denotes the dissipated energy due to mechanical damping. M, K s , B d stand for the equivalent mass, stiffness, and damping coefficients of the vibratory system.

Stability of the fixed points
In order to determine the fixed points in the absence of base excitation, the motion equations (Eq. (13)) are transformed into the following phase space form: The fixed points are defined by the vanishing of the vector field; say F S; t ð Þ ¼ 0. The place in the phase space where this condition is satisfied is known as stationary or equilibrium point. Equilibrium points are usually available in the absence of time-varying forces, and the stability of any fixed point is examined by the determination of the eigenvalues of the socalled Jacobian matrix A, at the corresponding fixed point as follows: . . .
where F 1 S 1 ; S 2 ; . . .; S 2n ð Þ ; . . .; F 2n S 1 ; S 2 ; . . .; S 2n ð Þare the components of the vector field F. When the eigenvalues of the Jacobian matrix at a particular fixed point have negative real parts, then the corresponding fixed point is asymptotically stable.

Stability of the periodic orbits
We have applied shooting method for autonomous systems to capture the periodic motions and determined the stability of the periodic orbits by means of Floquet theory [2].

Results and discussion
The mechanical, electrical, and geometrical properties of the harvester are given in Table 1. The values given in Table 1 have been taken from [3]. The harvester is assumed to be made of silicon and CYTOP electret layer with breakdown voltage of 90 kV/mm [31]. Figure 4 illustrates the dimensionless tip deflection of the cantilever beam versus the applied electret surface voltage. It is shown that for a definite surface voltage, there exist either one or three equilibrium positions. In case merely one equilibrium position is available, the eigenvalues of the Jacobian matrix evaluated at the equilibrium point are negative and accordingly the fixed point is stable. However, once there are three equilibrium positions, two of them are stable (only one of them is of physical interest as the second one lays outside the gap and hence is not physically possible) and the other one is unstable. As the surface voltage is increased, the stable and unstable positions approach each other and they coincide on the so-called pull-in voltage where the system undergoes a saddle node bifurcation and beyond it, both of the equilibrium positions disappear.
As it is shown in Fig. 4, contribution of the first two modes in Eq. (17) provides a reasonable convergence and accordingly in the rest of the study the effect of the third and the higher modes in the system response is neglected. Figure 5 depicts the time history and the phase diagram of the dynamic response of the energy harvester in the absence of the base excitation but exposed to initial disturbance. Here, the impetus is to examine the energy conservation by means of comparing the input energy with that captured throughout the output circuit. As depicted, once the initial disturbance is applied, the output circuit initiates transforming the mechanical energy into electrical and accordingly the motion amplitude dissipates [2].  Relative permittivity of electret layer,e e 2 Permittivity of free space, e 0 ( C 2 N:m 2 ) 8:85 Â 10 À12 Fig. 3 The equivalent electromechanical circuit for energy conservation examination In the absence of base excitation, the energy conservation is examined by comparing the consumed and harvested energies throughout the output circuit. Once the circuit is closed due to the electrostatic force, the microbeam bends toward the electret layer and the electrical charge on both sides of the capacitor increases and accordingly an electric current flows in the circuit until the capacitor is charged. During this time span, a part of the energy is harvested throughout the resistance. Table 2 gives the amount of consumed, harvested, and stored energy in the capacitor for various electret layer voltages. Figure 6 depicts the time response, phase plane, and charge distribution on the capacitor plate for different electret layer voltages and R ¼ 50MX.
Considering the base excitation, Fig. 7 illustrates the harvested power for different load resistances with 250V electret surface voltage. As we increase the electrical resistance, the maximum harvested power increases up to a particular resistance (R ¼ 15MX) beyond which the harvested power decreases; this is because the effect of the higher electrical resistance on the harvested power is counterbalanced with the current reduction in the circuit and accordingly the harvested power exhibits a Lorenzian-type behavior [2]. The solid and dashed lines represent the stable and unstable limit cycles. Capturing the unstable periodic solutions usually demands cumbersome iteration process as they have no attraction basins; some of the periodic solutions on the unstable branch have not been captured since they require a lot of iteration process.
For a given R as the frequency is swept forward, the system response undergoes a jump to the higher (a) (b) Fig. 5 The tip response, subjected to initial disturbance, in the absence of base excitation for V s =250 (V) a: time response, b: phase plane Table 2 The input and harvested energy for various output resistances and initial disturbances in the absence of base excitation Rdt (Harvested energy) (pJ) 1  amplitude periodic attractor in the cyclic fold bifurcation point where the branch of stable and unstable solutions collides and beyond it both of the branches disappear. As we sweep the frequency in the backward direction, the response undergoes another jump to the lower amplitude periodic attractor in the period doubling bifurcation point where the Floquet multipliers leave the unit circle through -1. The difference between the forward and backward frequency sweeps is regarded as hysteresis which is very common in the literature of nonlinear systems. In the region of multi-response solutions where the system exhibits more than one stable limit cycle, it is suggested that we push the system to vibrate on the highest amplitude limit cycle so that to harvest more energy; this is usually possible either by forward/ backward frequency sweep or by disturbing the mass by an appropriate electrostatic initial voltage which pushes the system to the basin of the attraction of the desired limit cycle. Figure 8 depicts the jumps in the time response which occur in the bifurcation points in both forward and backward frequency sweeps.

(a) (b)
(c) Fig. 6 The tip response, subjected to initial disturbance, in the absence of base excitation, a: time response, b: phase plane, c: charge distribution on the capacitor (parameters are in non-dimensional form) Figure 9 represents the frequency response curves for the same parameters as of Fig. 7. The results show that the system experiences higher damping effects for higher values of resistance. Therefore, by increasing the resistance, the nonlinear behavior of the beam is reduced and the amplitude-frequency diagram of the beam straightens back to the right.
The more is the harvested power, the less is the motion amplitude since more mechanical energy is extracted from the vibratory system. As the nonlinearity dominates the response, the bandwidth broadens and accordingly the harvested power within the bandwidth increases. Figure 10 depicts the frequency response curves and the corresponding harvested powers for R ¼ 15MX; V s ¼ 250V, and different base excitation amplitudes.
As the base excitation amplitude increases, the amplitude of the super-harmonic resonance of order two surges and accordingly the harvested power in this region rises. Further increase in the amplitude of the base excitation leads to further energy harvesting as a result of the amplification of the motion amplitude in the super-harmonic resonance. A significant observation is that by increasing the level of base excitation, the amplitude of the response increases until a separation happens in the frequency response curve and a band gap emerges. A band gap is emerged in the vicinity of each resonant frequency (e.g., one gap emerges at X ¼ 2:5 for 0:3\Y b \0:6, and the other is generated at X % 1:25 for Y b % 1:3). For Y b ¼ 1:2, the frequency response in the super-harmonic resonance region splits into two distinct parts and accordingly two more bifurcations including one cyclic fold (point A) and one period doubling bifurcation (point B) manifest. As we further increase the base excitation amplitude to Y b % 1:3, band gap appears at X % 1:25 in the frequency response curve where pull-in occurs for any initial condition. After the band gap has emerged, no significant increase is observed in the maximum amplitude of the vibration in response to further increase in base excitation amplitude. Instead, the emerged band gaps grow as the level of base excitation is amplified. This urges one to compromise between the benefit of a high amplitude vibration and the disadvantage of a wide band gap due to large base excitations. For the energy-harvesting application, once the device undergoes the excitation frequency in the vicinity of the super-harmonic resonance, it is preferable for the system to be excited with the highest possible amplitude before the system undergoes band gap where pull-in occurs.
The results of Fig. 10 illustrate which regions are robust for energy harvesting against the variation of base excitation amplitude. For instance, the vibration amplitude is relatively high in the vicinity of excitation frequency X ¼ 3 for lower levels of base excitation (e.g., Y b ¼ 0:3). Although it gets even larger as the base excitation is amplified, it is observed that the stable response at X ¼ 3 disappears as the band gap grows for higher values of Y b . It is also observed that the secondary resonance is not significantly contributing to the response of the system for lower levels of excitation. Increasing the amplitude of base motion results in a larger amplitude of secondary resonance at lower frequencies. That is, the energy harvesting can be made at lower frequencies for higher values of base excitation. However, the growth of pull-in band gaps due to the increase in the excitation amplitude should be taken into consideration. Figure 11 illustrates the frequency response curve for R ¼ 15MX and Y b ¼ 0:3. The corresponding force response curves for three different frequencies X = 2.5, 2.6, and 2.7 are illustrated in this figure. The number of stable and unstable periodic solutions depicted in force response curves are in accordance with those in the frequency response curve. As the motion amplitude increases (Fig. 11b), the stable manifold 1 (SM1) and unstable manifold 1 (UM1) on the force response curves get closer and they eventually intersect at the cyclic fold bifurcation point (CF) beyond which both solutions disappear. However, the stable manifold 2 (SM2) and unstable manifold 1 (UM1) emerge at the other cyclic fold bifurcation (d) Fig. 10 The frequency response in the vicinity of the primary resonance for different base excitations at R ¼ 15MX (solid and dot  dash lines show the  stable and  unstable branches,  respectively) point (CF) and the corresponding manifolds diverge as the motion amplitude increases.

Conclusion
We proposed a nonlinear capacitive energy-harvesting device which harvests the energy of ambient mechanical vibrations. The equations of motion were derived based the Euler-Bernoulli beam theory and then discretized to a reduced-order model by means of Galerkin method. The minimum required mode shapes to gain a reasonable converged response were determined. The periodic solutions were computed by shooting technique, and their stabilities were examined based on the values of the corresponding Floquet multipliers. We inspected the energy-harvesting enhancement in the multi-response region which revealed a considerable increase in the ratio of the harvested energy to the input mechanical energy. The bifurcation types were studied on the frequency response curves, and upward and downward jump in the cyclic fold and period-doubling bifurcation points in the forward and backward frequency sweeps were reported, respectively. We have observed that as the amplitude of the base excitation increases, the superharmonic resonance of order two gets more pronounced and the harvested power increases considerably within resonance region. As we further increases the base excitation amplitude, the amplitude of the limit cycles in the resonance region increased, but a band gap immerged in the frequency response curves where we could not capture any stable periodic attractors and accordingly any excitation in the band led to pull-in. We have demonstrated that the nonlinearity of the response broadens the bandwidth which in turn enhances the energy harvesting.

Declarations
Conflict of interest: The authors declare that they have no conflict of interest.
Data availability Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.  Fig. 11 a Frequency response curve for Y b ¼ 0:3 and at ¼ 15MX; force response curves for b :X ¼ 2:5; c:X ¼ 2:6; and d: X ¼ 2:7