Regular and chaotic vibration in a piezoelectric energy harvester

We examine regular and chaotic responses of a vibrational energy harvester composed of a vertical beam and a tip mass. The beam is excited horizontally by a harmonic inertial force while mechanical vibrational energy is converted to electrical power through a piezoelectric patch. The mechanical resonator can be described by single or double well potentials depending on the gravity force from the tip mass. By changing the tip mass we examine bifurcations from single well oscillations, to regular and chaotic vibrations between the potential wells. The appearance of chaotic responses in the energy harvesting system is illustrated by the bifurcation diagram, the corresponding Fourier spectra, the phase portraits, and is confirmed by the 0–1 test. The appearance of chaotic vibrations reduces the level of harvested energy.


Introduction
Broadband energy harvesting systems for many applications are often nonlinear, exhibiting such nonlinear phenomena as material nonlinearities, geometrical nonlinearities, multi-scale responses and the appearance of multiple solutions. These phenomena can be observed from the nonlinear time series analysis of simulated mathematical models or from measured system responses in experiments. It is well known that the efficiency of many engineered systems may be enhanced by operation in a nonlinear regime. Nonlinear vibrational energy harvesting shows such an advantage, making possible broadband frequency vibration energy accumulation and transduction into useful electrical power output.
A range of vibration energy harvesting devices have been proposed [1][2][3][4][5][6]. The key aspect of nonlinear harvesters is the use of a double well potential function, so that the device will have two equilibrium positions [7][8][9][10][11][12]. Gammaitoni et al. [8] and Masana and Daqaq [13] highlighted the advantages of a double well potential for energy harvesting, particularly when inter well dynamics were excited. The Duffing oscillator model has been used for many energy harvesting simulations, with the addition of electromechanical coupling for the harvesting circuit [14,15]. Electromagnetic harvesters with a cubic force nonlinearity have also been considered [16]. Litak et al. [17] and Ali et al. [18] investigated nonlinear piezomagnetoelastic energy harvesting under random broadband excitation. McInnes et al. [19] investigated the stochastic resonance phenomena for a nonlinear system with a double potential well. Gravitationally induced double potential wells in a system with a vertical elastic beam with a tip mass have been studied extensively [20,[23][24][25].
Recently, in the context of broad-band energy harvesting, bifurcations and chaotic vibrations have been studied in several papers. Cao et al. [26] studied chaos in the fractionally damped broadband piezoelectric energy generator in the system with additional magnets. Syta et al. [27] analysed the dynamic response of a piezoelectric material attached to a bistable laminate plate. It is worth to note that chaotic vibrations are, in most cases, characterized by moderate amplitude of vibrations and simultaneously give continuous spectrum of frequency, which can be useful to increase mechanical resonator durability.
In this article we discuss different solutions appearing in that system. Especially, intra-and interwell oscillations as well as periodic and chaotic vibrations lead to different efficiency in the energy harvesting. Therefore, we identify the properties of given solutions by using nonlinear methods.

Mathematical model and equations of motion
For nonlinear energy harvesting an inverted elastic beam is considered with a tip mass and the base is harmonically excited in the transverse direction. Only a summary of the equations is provided here; Friswell et al. [20] give a full derivation of the equations. Figure 1 shows the beam as a vertical cantilever of length L with harmonic base excitation zðtÞ ¼ z 0 cos xt. The beam carries a concentrated tip mass, M t , with moment of inertia I t , at the tip of the beam. The horizontal and vertical elastic displacements at the tip mass are v and u respectively, and s represents the distance along the neutral axis of the beam.
In the following analysis the beam is assumed to have uniform inertia and stiffness properties; a nonuniform beam is easily modeled by including the mechanical beam properties in the energy integrals. The beam has cross sectional area A, mass density q, equivalent Young's modulus E, and second moment of area I.
In this paper we assume that the tip mass is sufficiently large so that a single mode approximation of the beam deformation is sufficient. The displacement at any point in the beam is represented as a function of the tip mass displacement through a function for the beam deformation, wðsÞ, as v p ðs; tÞ ¼ v p ðL; tÞwðsÞ ¼ vðtÞwðsÞ: In this paper we will assume The equation of motion of the beam-mass system is derived in terms of the displacement of the tip mass using Lagrange's equations [20][21][22] as Fig. 1 Schematic of the inverted beam harvester. M t denotes the tip mass attached to the elastic beam, while v and u denote the horizontal and vertical displacements of the mass. P denotes an arbitrary point on the beam whose position is described by the coordinates s, v p , and u p . Piezoelectric patches are placed along the beam but are not shown here where the constants are given in the ''Appendix''. Damping may also be added to these equations of motion, for example viscous, material or aerodynamic damping.
The model of the piezoelectric patches is now considered. The mechanical stiffness and mass density of the piezoelectric layers may be included in the beam constants already derived. The electromechanical coupling constants are where L c is the active length of the piezoelectric material, which is assumed to start at the clamped end of the beam. For a unimorph configuration with excitation in the 31 mode, with thickness h c and width b c , the constant c c is where h is the thickness of the beam, d 31 is the piezoelectric constant, E c is the Young's modulus of the piezoelectric material and z is the effective neutral axis [28]. These expressions assume a monolithic piezoceramic actuator perfectly bonded to the beam; Bilgen et al. [29] considered the effect of the structure of a Macro-Fiber Composite (MFC) on the coupling coefficient, and also the effect of the bond and the insulating Kapton layers.
On the electrical side the piezoelectric patches may be considered as a capacitor, and the electrical circuit is represented by a resistive shunt connected across the piezoelectric patch. The electrical equation then becomes where R l is the load resistor and C p is the capacitance of the piezoelectric patch.

Equilibrium positions
The equilibrium positions with no forcing are obtained by setting the velocity and acceleration terms to zero in Eq. (3) to give, This equation has either one or three solutions, and v ¼ 0 is always a solution. Since where M tb is the tip mass so that the beam is about to buckle.

Numerical simulations and results
Following [20] we use the same set of system parameters, as given in Table 1. The beam-mass system is excited at the base with harmonic excitation. Figure 2 shows the equilibrium position of the tip mass, using the analysis described in Sect. 2.1, and shows that the post buckled response has two stable equilibrium positions above the critical tip mass M tb ¼ 10 g. In the vicinity of these equilibrium points small amplitude vibrations can appear provided that excitation is also present. To observe its influence on the system dynamical response we used the harmonic base excitation at a fixed amplitude, z 0 , and frequency, f ¼ x=ð2pÞ and varied the tip mass M t . The corresponding bifurcation results are presented in Fig. 3. As expected the dynamical bifurcation pattern follows the split in the equilibrium positions, highlighting the significant change above the critical value tip mass M tb ¼ 10 g. The character of the oscillations change by undergoing a series of bifurcations starting from the single frequency (in the limit for small M t ) to nonperiodic vibrations (for relatively large M t ), that correspond to single point and a widely distributed set of points in the bifurcation diagram, respectively. Looking more closely at the bifurcation diagram, one can identify periodic windows characterised by distributions of discrete points among the non-periodic solutions. It is also possible that in this regions several solutions are present. However we used the same initial conditions for every frequency, and thus our approach was not able to capture all of the solutions simultaneously.
Additionally, Fig. 4 shows the energy harvesting performance, i.e. the power output is plotted versus the tip mass. One can see that the power output increases above the value of M t corresponding to the critical buckling mass. Moreover the periodic windows visible in Fig. 3 give a higher power output.  Fig. 5a-f, respectively. One can easily see the evolution of the system response. Figure 5a shows a single frequency, while Fig. 5b shows symmetric oscillations with two frequencies indicating single well oscillations around the central equilibrium point (see Fig. 2). Figure 5c corresponds to a small amplitude single frequency non-symmetric solution indicating that the stable equilibrium point has already bifurcated. Finally, Fig. 5d-f show cross well oscillatory responses including periodic (Fig. 5e) and non-periodic cases (Fig. 5d, f). Furthermore, the displacement, velocity and voltage signals have the same generic nature. The non-periodicity in the kinematics are reflected in the voltage output (see Fig. 5d, f).
The cases b, d-f (cases numerated as in Fig. 5) are considered in more detail in the frequency domain by calculating the Fourier spectrum of the mechanical resonator displacement, as shown in Fig. 6. The corresponding phase portraits (the displacement-velocity phase portraits of the tip mass) are shown in Fig. 7. It should be noted that these qualitative measures clearly indicate the appearance of chaotic motion in cases d and f. Interestingly, the subharmonic visible in Fig. 5b corresponds to f 1 ¼ 0:167 Hz, and this is one third of the excitation frequency f ¼ 0:5 Hz. The nature of this solution is illustrated in Fig. 7b with three loops in the phase portrait and three isolated Poincaré points. The chaotic solutions are characterised by strange attractors with complex trajectories and the wide distributions of the Poincaré points. It should, however, be noted that these plots are projections onto the displacement-velocity plane (neglecting the voltage coordinate). The next step is to determine the quantitative parameter showing the chaotic character of the   solutions. This can be done by means of the 0-1 test [30,31] which is especially simple in systems with many degrees-of-freedom. The corresponding analysis and the estimation results will be presented in the next section.

The 0-1 test
The '0-1 test', invented by Gottwald and Melbourne [30,31], can be applied for any system of finite dimension to identify chaotic dynamics, but it is based on the statistical properties of a single coordinate only. Thus it is suitable to quantify the response where only one parameter is measured in time. The test is related to the universal properties of dynamical systems, such as spectral measures, and can therefore distinguish a chaotic system from a regular one using a single variable.
A particular advantage of the 0-1 test over the frequency spectrum is that it provides information regarding the dynamics in a single parameter value, similar to the Lyapunov exponent. However, the Lyapunov exponent can be difficult to estimate in any non-smooth simulated system or measured data [32]. Therefore the 0-1 test can provide a suitable algorithm to identify the chaotic solution [33,[36][37][38].
To start the analysis, we discretize the investigated time series vðtÞ ! vðiÞ using the characteristic delay time dt equal to one quarter of the excitation period 2p=x. This roughly indicates the vanishing of the mutual information [33,39]. Starting from one of the initial map coordinate v(i), for sampling points i ¼ 1; . . .; N, we define new coordinates p(n) and q(n) as ðvðjÞ À vÞ r v cosðjcÞ; qðnÞ ¼ X n j¼0 vðjÞ À vÞ r v sinðjcÞ; where v denotes the average value of v, r v the corresponding standard deviation, and c is a constant 2 ½0; p. Note that q(n) is a complementary coordinate in the two dimensional space. Furthermore, starting from bounded coordinate v(i) we build a new series p(n) which can be either bounded or unbounded depending on the dynamics of the examined process. Continuing the calculation procedure, the total mean square displacement is defined as The asymptotic growth of M c ðnÞ can be easily characterised by the corresponding ratio  In the limit as n ! 1 (in present calculations n ¼ n max ¼ 150 while N ¼ 1350) we obtain the corresponding values of K c for a chosen value of c. Note, our choice of n max and N limits (in Eqs. 4,5) is consistent with that proposed by Gottwald and Melbourne [34,35,40]. N; n max ! 1 but simultaneously n max should be about N/10.
It is important to note that the parameter c acts like a frequency in a spectral calculation. If c is badly chosen, it could resonate with the excitation frequency or its super-or sub-harmonics. In the 0-1 test, regular motion would yield an expanding behaviour in the (p, q)-plane [34] and the corresponding M c ðnÞ has an asymptotic growth rate even for a regular system. The disadvantage of the test, its strong dependence on the chosen parameter c, can be overcome by a proposed modification. Gottwald and Melbourne [34,36,37] suggested randomly chosen values of c are taken and the median of the corresponding K c -values are computed.
In the above, the covariance covðx; yÞ and variance varðxÞ, for arbitrary vectors x and y of n max elements, and the corresponding averages x and y respectively, are defined as covðx; yÞ ¼ 1 n max X n max n¼1 ðxðnÞ À xÞðyðnÞ À yÞ; Finally, the median is taken of the K c -values in Eq. (6), corresponding to 100 random values of c 2 ð0; pÞ. Such an average K-value can now be estimated for various tip masses, M t . Figure 8 shows hKi versus the tip mass. Here the chaotic solutions are clearly distinguished as hKi % 1, from any type regular solutions where hKi % 0. This is consistent with the previous qualitative methods based on the frequency spectrum (Fig. 6), phase portraits and the Poincaré maps (Fig. 7). It is worth noting that for all cases with maximal energy outputs (Fig. 4) their dynamical responses were periodic (Fig. 8).

Conclusions
The examined response vibrational energy harvester exhibited regular and chaotic oscillations. We have reported the corresponding transition in the system behaviour around the equilibrium bifurcation from a single point to two equilibrium points. We observe this transition by changing the tip mass. The system was excited by harmonic inertial force at a constant frequency, and the harvester responded at this frequency and also with additional subharmonics dependent on the value of the tip mass.
The mechanical vibrational energy was continuously converted to the electrical power through the piezoelectric patch. Interestingly, we observed the rise of the power output once the system changed from single well to cross barrier vibrations. Furthermore the appearance of the chaotic response lowered the power output considerably.
The bifurcation between regular and chaotic vibrations was reported using the bifurcation diagram, the corresponding Fourier spectra, and the phase portraits, and confirmed by the 0-1 test. We noticed that the appearance of chaotic vibrations significantly reduced the harvested energy.

Appendix
The constants N 1 to N 9 are given by