Response Identification in a Vibration Energy-Harvesting System with Quasi-Zero Stiffness and Two Potential Wells

Abstract

In this paper, the frequency broadband effect in vibration energy harvesting was studied numerically using a quasi-zero stiffness resonator with two potential wells and piezoelectric transducers. Corresponding solutions were investigated for system excitation harmonics at various frequencies. Solutions for the higher voltage output were collected in specific branches of the power output diagram. Both the resonant solution synchronized with excitation and the frequency responses of the subharmonic spectra were found. The selected cases were illustrated and classified using a phase portrait, a Poincaré section, and recurrence plot (RP) approaches. Select recurrence quantification analysis (RQA) measures were used to characterize the discussed solutions.

1. Introduction

The area of piezoelectric energy harvesting from sources of vibration has experienced intensive development in the recent years [1]. Many researchers have discussed various configurations of energy-harvesting systems that provide nonlinear effects on energy conversion [2,3,4,5,6]. The concept of ambient vibration energy harvesting is very useful in the case of power autonomous sensors, portable devices, or recharging corresponding batteries. Piezoelectric vibration energy-harvesters were proposed due to the fairly high density of energy used in the process of mechanical-to-electrical energy conversion [4]. Such devices consist of a mechanical resonator (most frequently a cantilever beam), an energy transducer (based on a piezoelectric layer), and an electrical circuit with the load resistor [7,8]. At the start of piezoelectric harvester development, linear resonators were proposed [7]. Unfortunately, such simple structures proved to be too sensitive to resonance conditions, which was very important in the context of ambient vibration source applications. Certain progress was made by including multiple-degree-of-freedom systems [9] to resonators, but major improvement resulted from the replacement of linear resonators with nonlinear resonators [10,11].

In contrast to linear harvesters, nonlinear harvesters are not fixed to a single resonant frequency. Due to the additional solutions, they can work efficiently in some frequency ranges with a fairly high-power output response. Nonlinear resonators modify the resonance characteristics and introduce additional coexisting solutions with changing frequencies. The main disadvantage rests in the fact that the optimal solution (leading to higher output) may not be the most probable one. Such nonlinear energy-harvesting systems can be modelled by various interactions including impacts and magnetic or pre-stressed mechanical systems with multiple stable equilibria [10,11,12,13,14,15]. Starting with small, non-resonant solutions around the stable equilibria, the resonator system can reach a higher response level for a number of resonance frequencies and finally overcome the potential barriers of having a larger orbit. Such orbits can appear for different frequencies including those belonging to subharmonic solutions [16,17,18,19,20].

In this paper, we continue the discussion on nonlinear energy-harvesting systems by analyzing the system with quasi-zero stiffness and two additional potential wells. Alongside the graphical recurrence plot (RP) representation [21], we propose an identification schema using recurrence quantification analysis (RQA).

2. The Model

The system for energy harvesting from vibrating mechanical sources is presented in Figure 1a. The system consists of a mechanical resonator with a cantilever beam, a piezoelectric transducer, and an electrical circuit. As a result of the piezoelectric deformation caused by the inertial force, the electromotive force is created inside the electrodes. It is worth noting that the electrodes are placed along the piezoelectric patches. Consequently, the electrical charge is moving through the electrical circuit, resulting in the voltage and power outputs on the load resistor. The special configuration of permanent magnets models the shape of the potential, forming a quasi-zero stiffness of the cantilever and restoring force characteristics with additional symmetric potential wells in the specific bent positions of the cantilever (Figure 1b).

Following the standard analysis, we focus on the single deflection mode of the cantilever beam (see Reference [6]). The mathematical model in the dimensionless form can be expressed by the following set of equations:

$$\{\begin{array}{l}\ddot{x}+2\zeta \dot{x}-\frac{1}{2}{x}^3\left(1-x^2\right)-\chi u=p{\omega }^{2}cos\left(\omega t\right),\hfill \\ \dot{u}+\lambda u+\kappa \dot{x}=0\hfill \end{array}$$

(1)

In the above equations, x and u denote the tip point of beam displacement (dimensionless) and the corresponding voltage output, respectively. Note that the corresponding dimensional form of the energy-harvesting device prototype can be revealed from the size of the resonator, which is related to h (Figure 1a) of the physical model. t denotes the dimensionless time defined as t = Ω_nt′, where t′ is the original time in seconds, while Ω_n is the natural free frequency of small oscillations around both of the stable equilibrium positions (x = ±1). By obtaining the excitation frequency (Ω), we define the dimensionless excitation frequency (ω), scaled by the natural frequency (Ω_n as ω= Ω/Ω_n). The remaining dimensionless parameters (Equation (1)) denote the mechanical damping (ζ), the amplitude of kinematic excitation (p), and the piezoelectric couplings (κ and χ), while λ denotes the reciprocal of the time constant of the electrical circuit.

The formulated mathematical model constitutes a formal basis for conducting studies in the area of energy harvesting. The model-testing results presented in this paper were obtained under the assumption of the mathematical model (Equation (1)) coefficients—ζ = 0.01, λ = 0.01, κ = 0.5, and χ = 0.05—and the harmonic excitation parameters—kinematic (dimensionless) forcing amplitude (p = 0.183) and forcing frequency (ω) swept in the interval [0, 4].

3. Results of Simulations and Recurrence Analysis

Using the set of dynamic equations of Equation (1) with the application of the Runge–Kutta fourth-order method (usually used in the dimensionless models) implemented in the MATLAB environment, a number of solutions were obtained. In particular, the frequency broadband effect was studied. Namely, the investigation into the dependence between the dimensionless excitation frequency (ω) and the nature of the dynamic system response was our main goal at this stage of the research. For the purpose of simplicity and to discuss various types of potentially coexisting solutions, the nodal initial condition was assumed for the voltage and minimal for the displacement of the tip point of the cantilever (x₀ = 1, see Figure 1b), while tip point velocity was changed randomly in the interval [−1.5, 1.5] to test the various levels of kinetic energy. In Figure 2, the bifurcation diagram determined for the voltage output is presented. On this diagram, a number of branches that formed above and below 0.2 can be observed. These branches make it possible to roughly distinguish between the intra and inter-well solutions (small and large orbits). This result is in contrast to the linear system in which only one resonant solution at ω = 1 is present. In the case of a nonlinear system, the additional transformation of the frequencies is a side effect of dynamic evolution. Consequently, we may expect not only the periodical solutions (related to a small amplitude of vibrations) but also an appearance of sub-harmonic and super-harmonic leading components in the solutions and non-periodic (chaotic) solutions with a higher amplitude of vibrations. Furthermore, as the amplitude of the solution is larger for the inter-well solutions, the most complex solutions have such a character. In the limited cases, such solutions can appear just below the amount of energy sufficient to overcome the potential barrier. This effect can be seen as an additional structure in the lower branches of a bifurcation diagram (Figure 2).

3.1. Analyses Using Phase Space Trajectories and Poincaré Maps

To clarify the nature of the solutions, 16 cases were selected for detailed study (see the points marked by letters A–P in Figure 2). The results of the further investigations are collected in Figure 3. In the left, center, and right panels of Figure 3, time series, phase portrait with Poincaré points, and corresponding voltage recurrence plots are presented.

The recurrence plots, which are presented in the right column of Figure 3, were estimated for the output voltage time histories (cases A–P). It should be stressed that the considered voltage time histories (Figure 3, left column) were determined for the dimensionless time (t), defined as t = Ωnt’, where t’ is the original time in seconds. Therefore, the estimated recurrence plots indicate the system recurrent properties for consecutive samples (i) of the output voltage dimensionless time histories.

The points A and B noticeably represent single-well solutions. Case A is characterized by a small orbit of the resonator (period one solution represented by one Poincaré point per one loop in the corresponding phase portrait). The small vibration amplitude indicates the non-resonance vibrations. Case B is a little more complicated and contains the super-harmonic component with a larger amplitude close to a ½ period resonance (one Poincaré point for two loops in the corresponding phase portrait). In a similar way to case A, the solution of case B is an intra-well solution.

The solution for case C has a fairly large amplitude. Case C is an inter-well solution of the period one (represented by one Poincaré point for one loop in the corresponding phase portrait). The large amplitude results in deformation of the phase portrait. Conversely, solution D is non-periodic (chaotic, as indicated by the wide distribution of the Poincaré points). Case D is a large orbit inter-well solution but the hopping between the potential wells is irregular in time. The solution for case E is similar to the case C solution. In fact, they are placed on the same branch in the bifurcation diagram (Figure 2). While increasing the frequency above ω = 2, we observed that the chaotic branch is still present. Consequently, the case F solution is chaotic again. Interestingly, case G corresponds to the intra-well solution represented by three Poincaré points per four loops in the phase portrait. Presumably, this is a resonant solution of a ¾ period. Case H is a larger orbit solution represented by six Poincaré points per two loops. This is effectively a dominating sub-harmonic component solution similar to the period three solution with a small deformation in each second loop due to a period doubling bifurcation.

Case I belongs to the chaotic branch of solutions. It is similar to the case D and F solutions. Interestingly, case J is a large orbit sub-harmonic-like solution with four Poincaré points per two loops in the phase portrait. This solution breaks the symmetry through the central vertical axis. Note that this solution effectively coincides with two Poincaré points per one orbit (with a period doubling bifurcation). Case K is again an intra-well solution, although it is characterized by four Poincaré points per four loops in the phase portrait. Note that the loops are of various sizes. Conversely, cases L and M, representing the large orbits, are similar to previously discussed cases H and J, respectively. They are related by the period doubling bifurcation.

Evidently, cases H and L and, separately, cases J and M, lie along the corresponding branches for sub-harmonic (period two and period three solutions, respectively) larger orbit solutions (see Figure 2). The branch family of period three (with cases H and L) is completed by case 0 represented by nine Poincaré points per three loops in the phase portrait. Case N is the simple, small orbit harmonic non-resonant solution, while case P is represented by nine Poincaré points per three loops. The case P solution differs from the case J and M solutions by the attractor represented by the specific distribution of the Poincaré points, and the voltage output is visible in Figure 2 (see an additional branch). Classification of the analyzed cases A–P is gathered in

Table 1. Classification of cases A–P. Large orbits consistently have an inter-well oscillation character, while small and medium orbits are located in one of the potential wells (intra-well oscillations) and differ by the amplitude of oscillations.

Case Annotation	Classification
A	small orbit intra-well solution
B	medium orbit intra-well solution
C	large orbit inter-well solution
D	large orbit non-periodic (chaotic solution)
E	large orbit inter-well solution
F	large orbit non-periodic (chaotic solution)
G	medium orbit intra-well solution
H	large orbit sub-harmonic solution
I	large orbit non-periodic (chaotic solution)
J	large orbit sub-harmonic solution
K	medium orbit intra-well solution
L	large orbit sub-harmonic solution
M	large orbit sub-harmonic solution
N	small orbit harmonic non-resonant solution
O	large orbit sub-harmonic solution
P	large orbit sub-harmonic solution

.

3.2. Recurrence Plots (RPs) and Recurrence Quantification Analysis (RQA)

In order to identify the dynamic system properties, the RP method was used for analyzing the simulated time series. This approach was introduced by Eckmann et al. [22] for qualitative studies of recurring patterns, non-stationarities, and bifurcations of the dynamic attractor (represented also by the Poincaré section). Significant development of the RP method was provided by Webber et al. [23] and Marwan et al. [24] by introducing statistics of diagonal and vertical lines from the RP diagrams. In their studies, this method is known as recurrence quantification analysis (RQA) [24,25,26]. The recurrence quantifiers illustrate the quantitative information contained in the recurrence plots. As the signal is represented by the states defined in the previously reconstructed embedded phase space, the RP-based methods reveal all instances (recurrences) when the phase space trajectory of the dynamic system coincides at various times. Consequently, two points on the system trajectory are considered as neighbors if they are close enough to each other. The corresponding distance matrix (R) includes phase states for all sampling times represented by the matrix element, $R_{ij}^{\epsilon }$:

$$R_{i,j}^{\epsilon }=H\left(\epsilon -\left|\left|u_i-u_j\right|\right|\right), i,j=1,\dots , N.$$

(2)

where N denotes the number of sampled time events (voltage, u_i) while ε communicates the threshold distance. The above measures are used in a so-called embedding space obtained with the help of delayed coordinates using a time delay (τ) obtained from a reduction in the mutual information and dimension (m) estimated from the disappearance of the false nearest neighbors.

In reference to the results in Figure 3 (on the left panel), we observe that the RPs are sensitive to changes in the system dynamic response. Cases A and N are characterized by the strait diagonal lines that correspond to the harmonic responses represented by a circular shape of the phase portraits. The period of response vibrations coincides with the space between the lines, and this period is constant. Note that the information about the character of lines, including its continuity, straightness (or their corrugations), and distance between the local diagonal lines (shorter and longer), are key to identifying the solution. The extreme examples are the cases with chaotic responses (cases D, F, and I). These cases demonstrate fluctuations in the line lengths and properties such as thickness and straightness. In addition, some of the lines change in direction from diagonal to vertical or horizontal. Apart from the chaos, this could be also the effect of modulation and non-stationarities. The cases presented above are stationary, but some features are visible in chaotic motion, which breaks the linear superposition rule of solution combinations, favoring intermittency with a local switch between solutions of characteristic periodicities together with their short transients. Conversely, the local corrugations (deviations from the straitens) can be identified in the solution with orbits (defined by portraits in the phase plane) far from the circular solutions. This feature can be noticed in all large orbit solutions, especially in those with dominating sub-harmonic components (see cases H, J, L, M, O, and P).

The leading frequency coincides with the distance between the corresponding diagonal lines. Note that such lines are not always continuous (see cases C, E, and P). This effect is strongly dependent on the value of ε but could still be included into the statistics. Another class of solutions are resonant intra-well solutions wherein the continuous lines coexist with the dashed lines (see cases B, G, and K). Note that, in cases G and K, the tendencies of perpendicular diagonal lines are present as well. Such features can signal super-harmonic components in the solutions.

As mentioned previously, more information can be obtained from the statistical approaches using the recurrence points, namely RQA [24,25,26,27]. In the following analysis, we focused on select parameters including the average length of the diagonal lines excluding central diagonal and singular points in RP (La), the average and maximum lengths of the vertical lines (TT and V_max), the average period (T₂) and the corresponding entropy of the period distributions, clustering (informing about the coordination number of the black points with respect to their neighbors, Clust), and transitivity (informing about the average number of the paths based on the black points across the RP, Trans). The last two parameters are topological ones. They have meanings closely related to the dimensions of the dynamic attractor. Obtained results of the RQA analysis are gathered in the

Table 2. The recurrence statistics (RQA) for ε = 0.1. For comparison, in the first two columns, we included information about the number of loops and the main periodicity obtained from the phase portraits and Poincaré points of the mechanical resonator response (see the central panel in Figure 3). P_n denotes the period of the solution as the multiple (n) of the excitation period, and L_n and sL_n denote the number (n) of large or small loops respectively. The calculations of the RQA parameters were performed in the CRP toolbox for MATLAB [28].

Point	Number of Loops	Periodicity	La	TT	Vmax	T2	RTE	Clust	Trans
A	sL1	P1	38.56	9.68	11	96.91	0.127	0.749	0.747
B	sL2	P1	15.97	6.53	11	63.01	0.655	0.665	0.663
C	L1	P1	107.20	9.69	11	96.93	0.100	0.740	0.737
D	L > 10	P > 10	28.53	29.57	206	295.80	0.850	0.625	0.624
E	L1	P1	106.55	9.68	11	96.93	0.100	0.740	0.738
F	L > 10	P > 10	35.28	27.20	115	272.07	0.777	0.647	0.645
G	sL3	P3	60.22	17.58	32	172.63	0.673	0.682	0.680
H	L2	P6	164.01	27.25	34	273.13	0.434	0.730	0.728
I	L > 10	P > 10	30.05	24.21	177	242.29	0.712	0.637	0.635
J	L2	P4	137.43	18.68	25	187.27	0.431	0.740	0.738
K	L4	P1	56.97	19.87	39	197.75	0.741	0.667	0.665
L	L1	P3	287.12	27.22	31	272.84	0.243	0.739	0.737
M	L1	P2	152.85	18.66	21	187.17	0.220	0.741	0.739
N	sL1	P1	186.87	9.67	11	96.90	0.126	0.737	0.735
O	L6	P2	40.72	17.72	35	176.35	0.747	0.701	0.699
P	L9	P3	104.23	27.31	51	273.65	0.480	0.647	0.645

.

Based on the statistical results for ε = 0.1, chaos (cases D, F, and I) is indicated by the relatively small value of La and fairly large values of TT and Vmax (see also the small values of Clust and Trans and the large values of RTE and T₂). There are some changes between cases but the values of several RQA measures have to be checked to identify such solutions correctly. In contrast, in the case of small orbits, the values of La, TT, and V_max are small.

For better clarity, we also studied the RQA parameters estimated for a smaller ε. In

Table 3. The recurrence statistics. The same RQA parameters as in Table 2 but for ε = 0.01.

Point	La	TT	Vmax	T2	RTE	Clust	Trans
A	1687.31	0	0	103.76	0.025	0.993	0.962
B	870.25	2.00	2	102.87	0.029	0.990	0.961
C	694.40	2.00	2	102.89	0.030	0.986	0.957
D	8.82	6.52	31	660.66	0.707	0.600	0.582
E	694.80	2.00	2	102.89	0.030	0.986	0.957
F	10.35	5.68	28	573.32	0.659	0.580	0.563
G	120.38	3.02	4	274.39	0.124	0.735	0.713
H	43.77	3.29	7	297.80	0.409	0.675	0.654
I	9.54	5.57	31	562.57	0.671	0.594	0.576
J	52.90	3.78	5	355.33	0.140	0.655	0.635
K	38.65	3.71	5	339.36	0.210	0.652	0.633
L	518.36	2.96	3	275.73	0.077	0.750	0.726
M	41.27	2.05	3	199.35	0.101	0.339	0.328
N	430.37	2.00	2	102.89	0.034	0.971	0.942
O	20.20	3.65	7	330.62	0.527	0.667	0.647
P	94.61	6.63	11	685.11	0.223	0.716	0.693

, we present the corresponding results. Note that V_max worked very well regarding the detection of chaotic solutions for its large values. T₂ also indicated chaos for large values but not as precisely as V_max. It is also suitable to distinguish small orbits in most cases.

Conversely, entropy in RTE indicated more complex solutions with sub- and super-harmonic components apart from the excitation frequency. Clust and Trans worked fairly well, demonstrating small values for more complex solutions of multiple loops including chaotic solutions. La also worked properly, indicating chaotic and multiple-loops solutions for its smaller values in most cases.

4. Summary and Conclusions

In this study, we reported the frequency broadband effect on nonlinearity of vibration energy harvesting. We used a model with two potential wells and quasi-zero stiffness appearing as a flattening effect in the potential with respect to the central position of the resonator beam. The quasi-zero stiffness effect enables the resonator to produce a larger amplitude and a higher voltage output due to the broader distribution of the potential well and to reduce the potential barrier in comparison to the classical Duffing model [20].

Namely, in our treatment, the potential wells are separated by the region of almost zero potential force interaction (Figure 1b) that favors acceleration in the kinematic mechanism of excitation. Note that the polynomial approximation is the simplest possible one.

The simulation results are based on the polynomial expansion of the potential. We used the lower-order approximation to the potential consistent with the quasi-zero stiffness at the center of the displacement axis coinciding with the unbent beam and two symmetric potential wells located on the sides. The proposed system can be made with the help of magnets distributed in the structure. Similar structures of the potential could be observed in the levitating magnets approach with combined magnetic and mechanical springs [29].

The linear energy-harvesting systems are ruled by the impedance matching principle [30,31]. According to this principle, optimal conditions for energy harvesting are reached when the amount of electrical dissipation energy (dissipation through the electrical load coinciding with transduced energy) is equal to the mechanical dissipation energy. In nonlinear systems where multiple solutions are present, this principle does not operate properly.

The bifurcation diagram (Figure 2), which indicates many branches of different solutions, is presented as our main contribution. We proposed to identify solutions using the phase portraits and Poincaré points of the resonator motion. To complete the investigation, we adopted the approach of recurrence plots and RQA. The advantage of the recurrence analysis rests in the fact that knowledge about the resonator motion is not required as the investigation is based on the voltage output time history only. In our study, RQA was used in order to support the interpretation of the informative content (patterns) of the estimated RPs. In practical applications, if the system is not expected to exhibit chaotic or transient behavior, the use of five basic, so-called “classical”, RQA measures is sufficient [32]. In the case of nonlinear systems, the RQA measures based on the vertical and horizontal elements of the RPs are also used [33]. Interestingly, the complete identification is not possible based on a single parameter [34,35,36]. However, we offer the argument that multi-parameter identification is possible.

In the course of the RQA performed for varying threshold ε values, the multiscale system properties [37] were revealed. During our analysis, two different thresholds of ε were used (0.1 and 0.01). On the basis of the obtained results, it can be stated that the RQA method works similarly for both thresholds. It should be stressed that a low threshold value would be necessary if the embedding is fixed for all cases; this is unlike our study as we estimated the individual embedding for each case separately.