Based on the assumed numerical data characterizing the mathematical model of the analyzed energy-harvesting system, numerical experiments were performed to investigate the impact of the d shift on the efficiency of energy harvesting. In the first stage, the impact of the shift on the location of the zones of chaotic and periodic solutions was assessed. Among the periodic solutions, small and large orbits were found. Areas of chaotic solutions can be identified via different numerical tools, such as bifurcation diagrams [
21], 0–1 test [
19,
24] or by determining the maximal Lyapunov exponent [
25]. In our numerical simulations (performed in Mathematica), the areas of chaotic and periodic solutions were visualized in the form of multicolor maps showing the distribution of the maximal Lyapunov exponent (
Figure 2). According to the authors, this approach makes it possible to look at system dynamics with a wide range of variability of the parameters characterizing the source of mechanical vibration affecting the energy-harvesting system. The maximal Lyapunov exponent was estimated in a three-dimensional phase space following the numerical procedure proposed by Wolf et al. [
26]. The phase space was spanned by the computed coordinates (
x,
,
u). All the included two-dimensional, multicolor distribution maps were plotted for the assumed nodal initial conditions:
x(0) = 0,
(0) = 0 and
u(0) = 0. Furthermore, two control parameters, ω and p, were selected to describe changes in the external excitation source. Technically, the small distance of the initial conditions between the tested trajectory and the reference trajectory was assumed to be ε(0) = 10
−5. To obtain a satisfactory resolution, the ranges of the control parameters ω and
p were divided into 500 subintervals.
It is worth mentioning that the positive values of λ relate to the chaotic dynamic response of the system, otherwise (negative, λ < 0) the system response is regular with the corresponding phase trajectories tending to stable points or periodic orbits. However, when λ approaches values equal to zero, we are dealing with the so-called bifurcation points (or quasiperiodic solutions). In a low frequency range ω, narrow repetitive zones of chaotic solutions are arranged along the curves that can be approximated by functions with a fairly high exponential growth. Regardless of the parameter
d value, the largest area of chaotic solutions is located in the central parts of the multicolor maps showing the distribution of the maximal Lyapunov exponent. For the case of
d = 0, this region is located in the band ω [
1,
2] (see
Figure 2a). However, for
d = 0.3, this zone stretches towards higher values of the dimensionless frequency ω [1, 2.4] (see
Figure 2b). An indicator expressed as the effective value of the voltage induced on the piezoelectric electrodes was taken as a measure of energy-harvesting efficiency. To determine the effect of the solution on the efficiency of energy harvesting, the multicolor maps of the maximal Lyapunov exponent distribution were compared with the diagrams of effective voltage values. The results clearly indicate that when the solution changes from periodic to chaotic, the voltage induced on the piezoelectric electrodes is limited. Examples of these landmarks are marked in black. The RMS(
u) diagrams (RMS—root mean square) were identified in relation to every value of the dimensionless excitation frequency from the range ω [0,4], based on a time sequence of 50 periods of the mechanical vibration signal affecting the energy-harvesting system.
With regard to the zero initial conditions, it can be observed that—irrespective of the displacement value of the “cut” halves of the potential barrier d (
Figure 1b,c) and the amplitude of the excitation source p—periodic solutions with a periodicity of 1T (single excitation period T) dominate in the range of low, dimensionless excitation frequencies of
ω < 0.5. It should be noted that increasing the width of the hysteresis loop reduces the energy-harvesting capability. Consequently, the reduction in the voltage induced on the piezoelectric electrodes is particularly noticeable in the range of higher parameter
p values. In the range of low amplitudes of forced vibration,
p < 0.1, particularly in the zone of chaotic solutions, no significant decrease in effective voltage value is observed. As predicted by RMS(
u), the voltage value is directly proportional to the amplitude and average kinetic energy of cantilever beam vibration.
3.1. Influence of the Parameter d on the Geometrical Structure of a Chaotic Attractor
The below graphs (
Figure 3) show examples of attractors that occur in the regions of chaotic solutions. It is worth mentioning that in the classic visualization of the Poincaré cross-section, numerical results are presented as a set of points located on the phase plane. Much more information about the geometric structure of chaotic attractors can be obtained by analyzing the density of the distribution of the points of intersection of the trajectory and the control plane. In this way, one obtains information from the areas of the phase plane where the trajectory most often intersects with the control plane. The graphical representations of the Poincaré cross-sections are then plotted against bifurcation diagrams. From a theoretical point of view, bifurcation diagrams can be drawn based on the following: Poincaré cross-section, phase trajectory and time sequence. In our case, these diagrams were plotted based on the local minima (points marked in red) and maxima (points marked in blue) of the analyzed time series (
Figure 3).
The plotted Poincaré cross-sections show different geometrical structures of attractors depending on the external excitation frequency value affecting the energy-harvesting system. To quantify the geometric structure of the chaotic attractors, their correlation dimension was identified, and the effective voltage values induced on the piezoelectric electrodes were also given. Based on the model tests, it was found that in the case of the “developed” attractors (a developed attractor is understood as a structure consisting of the Poincaré cross-section points forming an even distribution along the geometric structure), the correlation dimension of the attractor of the system with a displaced characteristic (
ω = 1.7
Figure 3b) is comparable to the attractor plotted for the classical potential function (
ω = 1.7
Figure 3a). In conclusion, the parameter
d does not significantly affect the position of the band of chaotic solutions. Both for
d = 0 and
d = 0.3, unpredictable solutions are located in comparable ranges of the variability
ω ϵ [1.5, 1.9]. Similar values of the correlation dimension (DC) also result from the similarity of the geometric structure of the chaotic attractors and the density distribution of the points of intersection between the trajectory and the control plane.
The most important element distinguishing the two discussed cases is the parameter representing the effective value of the voltage induced on the piezoelectric electrodes. In the case of the “split” and shifted characteristic (
ω = 1.7 in
Figure 3b), the value is lower by approximately 0.08. This RMS voltage value is directly related to the lower amplitude of the displacement of the free end of the flexible beam I (
Figure 1). The influence of the displacement of the “cut” halves of the potential barrier is best visible in the bifurcation diagrams. In particular, this is very visible in the range of low values of
ω, where increasing the parameter
d results in narrowing the zones of chaotic solutions. It is also worth noting that in the range of low values of ω, the zones of unpredictable solutions are additionally shifted, as a result of which it is impossible to directly refer to the Poincaré cross-sections that were plotted for the cases
d = 0 and
d = 0.3. The numerical modelling results presented in this section were obtained for zero initial conditions. This approach however does not fully reflect the capability of harvesting energy from vibrating mechanical devices. Therefore, the numerical simulations focused on the assessment of initial conditions from the point of view of coexisting solutions.
With reference to the presented Poincaré sections, Fourier amplitude–frequency spectra were plotted. They provide information on the dominant harmonics in the time series, which constitutes the formal basis for depicting Poincaré maps. Based on the spectra, it is possible to conclude that in the case of developed attractors, for which the correlation dimension DC > 1.5, harmonic components representing the frequency of the source of excitation dominate in the Fourier spectra. We deal with such geometrical structures of Poincaré cross-sections in the widest zone of chaotic solutions, located in the band ω ϵ [1.4, 1.9]. However, in the range of low ω values, the correlation dimension of the Poincaré cross-sections DC <1.5. At the same time, in the amplitude–frequency spectra, it is represented by the domination of harmonic components that are a multiple of the frequency of mechanical vibrations affecting the energy-harvesting system. We deal with such a sequence of dominant harmonic components of the Fourier spectrum when the correlation dimension of Poincaré cross-sections is in the DC range of [1.1, 1.5]. At this point, it is worth noting that such spectra occur both in the case of continuous and smooth (p = 1, d = 0, ω = 0.42) and discontinuous (p = 1, d = 0.3, ω = 0.26) characteristics representing the potential barrier. On the basis of the presented results of numerical experiments, it is also possible to state that for correlation dimensions DC < 1.1, in the amplitude–frequency spectra, the dominant harmonic components are the multiples of the combination of the two fundamental frequencies ω and ω1. As in the previous example, this is the case with both continuous and smooth (p = 1, d = 0, ω = 0.21) and discontinuous (p = 1, d = 0.3, ω = 0.48) characteristics representing the barrier potentials.
The diagrams (the bottom panels
Figure 3a,b) show a direct comparison of Feigenbaum’s steady state bifurcation diagrams with the diagrams of RMS values of the voltage induced on the piezoelectric electrodes. In the wide range of variability of the dimensionless excitation frequency
ω ϵ [0, 4], bifurcation diagrams are dominated by one relatively wide band of chaotic solutions, within the
ω range of [1.4, 1.9]. On the other hand, in the band
ω < 1, there are many bands corresponding to unpredictable solutions, separated by pieriodic areas and bifurcation zones. The individual bands characterizing the dynamics of the tested energy-harvesting system were distinguished by colors: the bands of chaotic solutions were depicted with a light cyan color, and the bifurcation zones with a light magenta color. Based on the presented graphs, it is difficult to unambiguously characterize the dynamics of the system, because both in the bands of chaotic solutions and bifurcation zones we deal with an increase and a decrease in the effective voltage induced on piezoelectric electrodes. We are also dealing with both the decrease and increase in the voltage induced on piezoelectric electrodes in the areas of periodic solutions. For example, the voltage drop induced on piezoelectric electrodes in the area of periodic solutions is in the ω band of [0.225, 0.245] (
Figure 3c). The same is the case when the halves of the “cut” potential barrier overlap with
ω in [0.23, 0.255] (
Figure 3d), and the range of variation of the dimensionless excitation frequency is similar. In the remaining bands of periodic solutions, for example
ω ϵ [0.337, 0.395] (
Figure 3c) and
ω ϵ [0.355, 0.42] (
Figure 3d), we deal with the opposite situation, i.e., an increase in the voltage induced on the electrodes is observed. While any detailed inference regarding the correlation of the bifurcation diagram with the diagram of rms voltage values in the bands of chaotic solutions and bifurcation zones is difficult due to the large number of points, in the case of periodic solutions it is possible to formulate a probable diagnosis.
It is possible to formulate the hypothesis that if the mean of the slope coefficients approximating the branches of the bifurcation diagram in the area of periodic solutions is positive, then we are dealing with an increase in the effective value of the voltage induced on the piezoelectric electrodes. On the other hand, when the average value of the slope approximating branches of bifurcation diagrams takes negative values, a decrease in the RMS voltage is observed on the electrodes. If the mean value of the slope of the branch takes values in the vicinity of zero, then the RMS voltage diagram assumes a constant value.
3.2. Identification of Multiple Solutions
It should be emphasized that the results presented in the previous section pertain to selected cases of system dynamics for zero initial conditions. However, one of the most important properties of nonlinear dynamic systems is related to the coexistence of different solutions depending on the initial conditions. For the defined system and excitation parameters, one has to find the system’s evolution with various initial conditions. This problem comes down to the study of phase plane orbit topologies [
22], the origins of which are located in different places. In the case of systems with a greater number of degrees of freedom, coexisting periodic and chaotic solutions are identified in a multidimensional phase space. The space dimension is a multiple of the number of degrees of freedom. Considering the research problem formulated in this way, it is possible to conduct complex numerical calculations, the results of which can be illustrated in the form of a diagram of solutions (DS) [
27,
28]. In this approach, the information about the efficiency of harvesting energy from coexisting solutions is neglected. For this reason, this work contributes to the state of the art by proposing a different approach, the essence of which boils down to multiple plotting of voltage effective value diagrams. Every diagram is plotted for randomly chosen initial conditions of the phase space. Such presentation of computer simulation results provides information about the number of coexisting solutions. However, it does not provide information about their periodicity. For this reason, additional detailed computer simulations were performed to identify the periodicity of individual branches appearing in the diagrams of effective voltage induced on the piezoelectric electrodes (
Figure 4). The following convention was adopted to define periodicity: the digit before the letter T indicates the periodicity of the solution, while the number of solutions with a given periodicity is denoted by the right subscript.
It should be noted that both the plotted DS diagrams and the diagrams of effective values plotted for randomly selected RMS(
u) initial conditions were obtained in a simplified manner. As a result, their accuracy is a compromise between the precision of numerical calculations and the time of computer simulation. Below are some examples of diagrams showing the effective values of the voltage induced on the piezoelectric electrodes. On their basis it was possible to determine the ranges of variability of the dimensionless frequency ω in which energy harvesting is most effective. Examples of numerical results showing the effects of the dimensionless excitation amplitude p and the shift d of the left and right potential halves overlap in the barrier zone are given in the graphs (
Figure 4).
Given the identified branches on the plotted diagrams of the RMS values of the voltage induced on the piezoelectric electrodes, additional numerical simulations were performed to identify phase trajectories of the coexisting solutions. The periodicity of individual branches was identified for the frequency
ω > 0.4. Irrespective of the
d shift value of the “cut” potential barrier, the highest energy-harvesting efficiency was obtained for 1T-periodic solutions in the band
ω < 2. A comparison of the numerical results presented in this section with the RMS(
u) diagrams (
Figure 2) reveals that the highest energy-harvesting efficiency in this band is achieved by assuming zero initial conditions. However, in the zone
ω > 2, energy-efficient solutions can be obtained too, for example by initiating appropriate initial conditions. It is also worth noting that the parameter
d does not affect the nature (periodicity) of a given solution. Its influence is mainly visible in the shift of individual branches in relation to the frequency axis ω. However, with increasing the parameter
d value, this shift is directed towards higher values of the dimensionless excitation frequency ω. Moreover, as the parameter
d is increased, the efficiency of energy harvesting is reduced. A detailed examination of the single points and branches in the diagrams (
Figure 4) shows that they represent transient solutions which finally become attracted to permanent periodic solutions over long enough time spans. This is the case with, e.g., the 4T
2 branch that appears in the diagrams identified for
p = 2.
For clarity, the selected cases were studied with the help of phase plane trajectories (
Figure 5). The images are shown against the background of the potential barrier. It should be noted that the multiple solutions, which are marked in different colors depending on the solution periodicity, occur with the evolution of the frequency
ω. The Poincaré points are also identified, which confirms the periodicity of the solutions obtained. The multiple solutions of the same period are plotted using the same color. Usually, they are trapped in the left and right potential wells.
The trajectories, which were depicted against the background of the “cut” and shifted halves of the potential barrier, illustrate the possible cases of coexisting solutions.
Figure 6 presents selected examples of coexisting solutions. The examples were selected to illustrate the effect of the superimposition of the “cut” halves of the potential barrier, with respect to large (
Figure 6a) and small (
Figure 6d) orbits of coexisting solutions and unpredictable responses (
Figure 6b). The following convention was adopted during the graphical visualization of the results of computer simulations: the colors of the coexisting solutions were assigned to the individual branches of the diagrams of the RMS values of the voltage induced on the piezoelectric electrodes (
Figure 4), while the coexisting solutions are plotted with dashed lines.
Based on the results of computer simulations, it can be concluded that the increase in the value of the parameter
d does not change the nature of the solution. This behavior of the system is observed both in the case of periodic and chaotic solutions. It is also worth noting that the vibrations of the elastic cantilever beam for periodic solutions do not undergo a phase shift with the increase in the parameter
d (
Figure 6a,c). In the case of large and medium orbits, only a reduction in the amplitude of vibrations is observed. However, amplitude limitations were not observed for solutions whose orbits are located inside the potential well (
Figure 6d). If the response of the system is given in the form of a chaotic solution (
Figure 6b), the change of parameter
d does not cause a significant deformation of the geometric structures of the Poincaré sections, and their similarity is preserved. The differences in the plotted Poincaré cross-sections become apparent when the correlation dimensions are identified. With an increase in the value of the parameter
d, a decrease in the value of the correlation dimension
Dc is observed. On the other hand, in graphical images of amplitude–frequency spectra, it is manifested by limiting the amplitude of the excited harmonics, with a simultaneous slight shift towards higher values.