High-Level Vibration for Single-Frequency and Multi-Frequency Excitation in Macro-Composite Piezoelectric (MFC) Energy Harvesters, Nonlinearity, and Higher Harmonics

This paper presents an extensive experimental investigation to identify the influence of signal parameters on a piezoelectric harvester’s performance. A macro-fibre composite energy harvester was studied as an advanced, flexible, high-performance energy material. Gaussian white noise, and single-frequency and multi-frequency excitation were used to investigate nonlinearity and multiple-frequency interactions. Using single low and high frequencies, we identified the nonlinearity of the harvester’s vibration. Multi-frequency excitation with a series of low-to-high-frequency harmonics mimicked the practical vibration signal. Under such multi-frequency excitation, the harvester’s nonlinear behaviour was studied. Finally, the interaction effects among multiple frequencies were identified. The results show that under pure resonant excitation, high-level vibration led to high-level mechanical strain, which caused nonlinear vibration behaviour. Moreover, it was shown that the different harmonics excited the various structure bending modes, which caused the nonlinearity of multi-frequency excitation. The first four harmonics of the real-time signal were important. The experimental results emphasise the resonant nonlinearity and interactions of multi-frequency excitation effects.


Introduction
Piezoelectric materials with electrical-to-mechanical conversion abilities are alternatives to power electronic devices, as they utilise the often-wasted electrical energy. Piezoelectric Vibration Energy Harvesting (PVEH) has found industrial [1] and bio-related [2] applications. By moving toward self-powered electronics with PVEH, electronic devices and biomedical sensors/actuators can be installed in inaccessible areas, such as in vivo or remote areas, using the often-wasted available energy [3].
Practical vibration signals from machines or the environment are often random, multifrequency signals. Despite extensive PVEH from practical vibrations [4,5], a laboratorycontrolled sinusoidal input is typically used as an excitation signal, which is different from real-time vibration systems, as real-time vibration signals are stochastic. Real-time randomness may affect the efficacy of piezoelectric harvester optimisation practices [6], emphasising the importance of multi-frequency vibration analysis of piezoelectric energy harvesters (PEHs). Some studies have analysed PEHs using random vibration methods, such as white noise [7] and Gaussian-coloured noise [8]. Moreover, PVEH using wind force, as another familiar environmental energy source, entails randomness, which has been experimentally investigated [9]. These studies often focus on the PEH output without properly focusing on different excitation frequencies and multi-modal beam models. In PVEH, the frequency and amplitude of the vibration source are prominent parameters of the output power of a PEH [10,11]. The PEH maximum voltage is a vibration signal with an excitation frequency close to the natural frequency [12]. This process is called frequency

The Need for Multi-Frequency Harmonic Analysis
Practical vibration sources are random vibration signals comprising a series of harmonic signals. The Fourier Transform (FT) of a time-domain vibration signal demonstrates the frequency information of the vibration source. Therefore, a full-range frequency analysis is required for energy harvesting power estimation. As a practical example, the water pump acceleration signal of a running water pump with rotation speed Ω = 2970 rpm (Ω ≈ 49.5 Hz) is shown in Figure 1a [22]. This acceleration signal is influenced by the roller bearing elements and structural elements, which affect the FT peak amplitudes due to factors such as bearing ball diameter and number of balls [23,24]. This acceleration signal is not purely harmonic and consists of harmonic multipliers of Ω, according to the FT signal shown in Figure 1b. In the zoomed-in view of the FT signal, the four dominant frequencies are shown and marked with 1 × Ω, 2 × Ω, 3 × Ω, and 4 × Ω.
Micromachines 2022, 13, x FOR PEER REVIEW 3 of 17 to the FT signal shown in Figure 1b. In the zoomed-in view of the FT signal, the four dominant frequencies are shown and marked with 1 × Ω, 2 × Ω, 3 × Ω, and 4 × Ω. A common assumption in piezoelectric harvester analysis is the single-mode assumption, which simplifies modelling by only considering a frequency range around the fundamental harvester frequency. Analytical modelling studies are often simplified for single-frequency harmonic vibrations, but practical vibration signals are not single-frequency harmonics. High-frequency harmonics may excite a harvester's high-vibration modes. Since an electrode covers the whole piezoelectric surface, electrical charge cancellation may be experienced in multi-vibration modes. Therefore, a full-range frequency analysis is needed to assess the high-frequency harmonics and interaction effects among high frequencies.
Moreover, most analytical models present a linear model for energy harvesting power estimation. The linear model assumes a linear relationship between input acceleration and output voltage. Nevertheless, assessing linearity under conditions of multiple vibration frequencies is difficult due to different vibration mode shapes. A common assumption in piezoelectric harvester analysis is the single-mode assumption, which simplifies modelling by only considering a frequency range around the fundamental harvester frequency. Analytical modelling studies are often simplified for single-frequency harmonic vibrations, but practical vibration signals are not singlefrequency harmonics. High-frequency harmonics may excite a harvester's high-vibration modes. Since an electrode covers the whole piezoelectric surface, electrical charge cancellation may be experienced in multi-vibration modes. Therefore, a full-range frequency analysis is needed to assess the high-frequency harmonics and interaction effects among high frequencies.
Moreover, most analytical models present a linear model for energy harvesting power estimation. The linear model assumes a linear relationship between input acceleration and output voltage. Nevertheless, assessing linearity under conditions of multiple vibration frequencies is difficult due to different vibration mode shapes.

Closed-Form Solutions for Mechanical and Electrical Responses
a n are the sampled acceleration signal measurements of acceleration function a(t). To obtain the closed-form solutions under this general load, the acceleration is represented by a series of harmonic functions using the Fast Fourier Transform (FFT).
where A i is the external acceleration FFT and N is the number of sampled measurement data.
For the vibration solution of the harvester, the piezoelectric beam deformation is denoted by where φ r (x) is the displacement-dependent function (beam mode shapes) and η r (t) is the time-dependent mechanical deformation function.
The electromechanical equations of a piezoelectric unimorph without tip mass are given by [15] ..
where the parameters are given in Table 1. The unimorph parameters are shown in Figure 2.
Beam width (m)   Assuming a linear framework, mechanical displacement ) and output voltage ) are the summation of the outputs from each harmonic Ω , starting from = 0 to +∞. Individual harmonic components are denoted with ̅ , and , , respectively. Assuming a linear framework, mechanical displacement η r (t) and output voltage V R (t) are the summation of the outputs from each harmonic Ω i , starting from i = 0 to +∞. Individual harmonic components are denoted with η r, Ω i and V R,Ω i , respectively. Therefore, the overall mechanical vibration and electrical responses are where η r, Ω i is the mechanical response and V R,Ω i is the piezoelectric voltage response due to a nominal excitation harmonic with Ω i frequency. Note that the over-bar indicates the magnitude. For obtaining η r, Ω i and V R,Ω i , the harmonic solution analysis of the piezoelectric energy harvester differential is carried out [25], as the steady-state relationships can be given by Eliminating mechanical response η r, Ω i between Equations (6.a) and (6.b), the output voltage can be expressed as Finally, the total voltage generation is calculated as The output power is calculated as This power conversion term is frequency and load dependent, in addition to being dependent on the material and geometrical properties. Equation (9) contains essential information: (1) The total output power is calculated using two series over the frequency range and the modal mode shapes for general input acceleration. Therefore, parameters such as Y r and Λ r are mode-shape dependent, in addition to being dependent on the external frequency, which is Ω i . The interaction between the external excitation frequency and the mode shapes is a complex research object, specifically when an electrode covers the piezoelectric layer. (2) Moreover, there is a linear relationship between the output power and the external square acceleration amplitude. This linear dependency over different frequencies needs to be researched, since excitation and modal mode interactions can be nonlinear.
This paper tackles the above two research questions. Figure 3a shows the test rig for all the experimental measurements. Two aluminium base plates connected a B&K shaker to the piezoelectric sample, tightened with four bolts. The amplifier was controlled with a National Instruments NI 9263 module, which generated analogue voltage signals. A KEPCO BOP 100-10MG amplifier amplified the signals and powered the shaker. The piezoelectric harvester output wires were connected across a resistive load of resistance R. A Data Acquisition (DAQ) system, an 8-channel National Instruments NI 9201 module, was employed for reading the voltage across the resistive load (which was also the voltage across the harvester). The NI modules were placed in an NI cDAQ 9172 chassis connected with a USB cable to the computer.

PVEH Device Optimum Characterisation
The piezoelectric sample was primarily analysed and characterised with the model and experimental tests presented in this subsection. The material properties for modelling were as follows: Young's moduli of MFC and aluminium layers were 15.85 [28] and 68.9 GPa, and the corresponding densities were 5540 and 2700 kg/m 3 . The damping coefficient was 5% [29]. The relative permittivity coefficient and piezoelectric constant d31 were 1800 and −170 × 10 −12 C/N, respectively.
For this energy harvester, the first undamped natural frequencies obtained with the current model were compared with those of the experimental tests, as reported in Table  2. The natural frequency difference was 1.3 Hz. The higher-mode natural frequencies are also given in Table 2 and were employed in the multi-modal analysis. This difference in the natural frequencies can be linked to the non-uniform piezoelectric MFC; in practice, the commercial MFC sample had an active area where the piezoelectric material was placed, and on the outer areas, there were only Kapton layers. Including this non-uniformity would have required advanced piezoelectric beam modelling, which was beyond the scope of this paper. For detailed demonstration and modelling of an MFC, one can The shaker input voltage signal, y(t), was varied in frequency and amplitude, and the piezoelectric voltage was measured. The shaker input voltage was a single-frequency or multiple-frequency harmonic.
The piezoelectric sample comprised a macro-fibre composite (MFC) of 0.3 mm in thickness, an aluminium substrate shim of 0.12 mm, and an epoxy rapid 332 bonding layer of~0.4 mm in thickness. The bonding layer has little effect on the natural frequency, but it is a significant source of damping [26]. Therefore, the bonding layer was neglected in the structural modelling of the natural frequency; however, its effect was considered by extracting the damping coefficient experimentally, so the bonding layer influence was observed. The bonding layer joined the MFC and substrate shim. The MFC was the M-8528-P2 type from Smart Material GmbH (Dresden, Germany) [27], a piezoelectric bending energy harvester. The MFC had seven sub-layers: two Kapton outer layers, two acrylic layers, two electrodes, and one central active piezoelectric layer. More information about the MFC can be obtained from Smart Material Inc. [27]. Figure 3b shows the MFC sample.

PVEH Device Optimum Characterisation
The piezoelectric sample was primarily analysed and characterised with the model and experimental tests presented in this subsection. The material properties for modelling were as follows: Young's moduli of MFC and aluminium layers were 15.85 [28] and 68.9 GPa, and the corresponding densities were 5540 and 2700 kg/m 3 . The damping coefficient was 5% [29]. The relative permittivity coefficient and piezoelectric constant d 31 were 1800 and −170 × 10 −12 C/N, respectively.
For this energy harvester, the first undamped natural frequencies obtained with the current model were compared with those of the experimental tests, as reported in Table 2. The natural frequency difference was 1.3 Hz. The higher-mode natural frequencies are also given in Table 2 and were employed in the multi-modal analysis. This difference in the natural frequencies can be linked to the non-uniform piezoelectric MFC; in practice, the commercial MFC sample had an active area where the piezoelectric material was placed, and on the outer areas, there were only Kapton layers. Including this non-uniformity would have required advanced piezoelectric beam modelling, which was beyond the scope of this paper. For detailed demonstration and modelling of an MFC, one can refer to [16]. Piezoelectric unimorph harvesters are characterised by a natural frequency and optimum power generation. The studied harvester was first evaluated over a frequency range to obtain these conditions. Figure 4a shows power with load resistances over frequency. The resonant frequency was constant for all the load resistances and was equal to ω r=1 = 20.4 Hz. This frequency was the first resonant frequency of the harvester. Moreover, there was another resonant frequency, ω pc = 31.0 Hz, associated with the piezoelectric structural effects, here called piezoelectric-coupled frequency. Second, output power versus load resistances were evaluated for optimum load resistance selection. Figure 4b shows output power versus load resistance at the resonant frequency. The optimal load resistance of 21.8 kΩ led to the highest power generation.  Piezoelectric unimorph harvesters are characterised by a natural frequency and optimum power generation. The studied harvester was first evaluated over a frequency range to obtain these conditions. Figure 4a shows power with load resistances over frequency. The resonant frequency was constant for all the load resistances and was equal to = 20.4 Hz. This frequency was the first resonant frequency of the harvester. Moreover, there was another resonant frequency, = 31.0 Hz, associated with the piezoelectric structural effects, here called piezoelectric-coupled frequency. Second, output power versus load resistances were evaluated for optimum load resistance selection. Figure 4b shows output power versus load resistance at the resonant frequency. The optimal load resistance of 21.8 kΩ led to the highest power generation.

Results and Discussion
The result section provides a comprehensive result set from the experimental tests. Through the experimental tests, power generation with = 21.8 kΩ was studied with various excitation signals. Note that the load resistance of = 21.8 kΩ was employed

Results and Discussion
The result section provides a comprehensive result set from the experimental tests. Through the experimental tests, power generation with R opt = 21.8 kΩ was studied with various excitation signals. Note that the load resistance of R opt = 21.8 kΩ was employed in all experimental results presented in this section. Excitation signals with Added White Gaussian Noise (AWGN), single-frequency harmonics, and multi-frequency harmonics were studied. Moreover, the effect of the vibration level from low-level to high-level vibrations was studied with harmonics with single frequency and multi-frequencies. Because in piezoelectric energy harvesting applications a harvester's natural frequency is matched to the dominant acceleration frequency, here, Ω = ω r=1 was assumed.

Effects of Adding White Noise to the Excitation Signal
Pure analytical solutions often consider a single harmonic signal without noise as an excitation signal. In contrast, practical vibration sources often have added noise. Here, the effect of adding white noise was studied by considering Added White Gaussian Noise (AWGN) for a single harmonic signal with matched frequency, i.e., Ω = ω r=1 . The shaker vibration signal was y(t) = 0.1(1 + AWGN) sin(2πΩt) (V), where "AWGN" is the white noise percentage, 1% or 2%. The shaker vibration signals with and without AWGN are shown in Figure 5a. The output power results showed that adding white noise increased the power slightly (see Figure 5b). White noise, i.e., random noise with a uniform frequency domain value, caused additional vibration on the piezoelectric harvester, so power generation slightly increased. An important conclusion is that white noise in the practical vibration data did not show a reduction in power influence. shown in Figure 5a. The output power results showed that adding white noise increased the power slightly (see Figure 5b). White noise, i.e., random noise with a uniform frequency domain value, caused additional vibration on the piezoelectric harvester, so power generation slightly increased. An important conclusion is that white noise in the practical vibration data did not show a reduction in power influence.

Effects of Increasing the Vibration Level with Single-Frequency Harmonics
A single-frequency harmonic excited the harvester, and the vibration-level effect on the output power was investigated by increasing the harmonic amplitude. This experimental setup tested the linearity assumption of power generation versus input square acceleration amplitude. Note that a single frequency does not imply constant frequency. Various single-frequency harmonics with different frequencies were studied. Considering the matched frequency (Ω = ω ) up to six harmonics, the natural multiples of the fundamental frequency were studied. Therefore, Ω = × ω , = 1, … , 6.
Shaker input signals are symbolised with ) = ), where is the magnitude of the shaker excitation signal and is the driving frequency multiplier. Subscript denotes different excitation magnitude levels for one driving frequency. By changing , the vibration level moves from low-level to high-level vibrations; in other words, the strain on the piezoelectric material changes from low levels to high levels.

Effects of Increasing the Vibration Level with Single-Frequency Harmonics
A single-frequency harmonic excited the harvester, and the vibration-level effect on the output power was investigated by increasing the harmonic amplitude. This experimental setup tested the linearity assumption of power generation versus input square acceleration amplitude. Note that a single frequency does not imply constant frequency. Various singlefrequency harmonics with different frequencies were studied. Considering the matched frequency (Ω 1 = ω r=1 ) up to six harmonics, the natural multiples of the fundamental frequency were studied. Therefore, Ω i = i × ω r=1 , i = 1, . . . , 6.
Shaker input signals are symbolised with y J K (t) = Y J sin(Kω r t), where Y is the magnitude of the shaker excitation signal and K is the driving frequency multiplier. Subscript J denotes different excitation magnitude levels for one driving frequency. By changing Y, the vibration level moves from low-level to high-level vibrations; in other words, the strain on the piezoelectric material changes from low levels to high levels.
The output power with R opt was recorded during three independent runs. Table 3 shows the output power, and average and experimental errors for all the harmonics and vibration levels. An experimental error of less than 4% implies experimental repeatability. The sensitivity comparisons showed the dramatic power increase obtained by magnifying the excitation amplitude of all excitation harmonics; however, the trends were not the same for all harmonics. Output power (P) versus square shaker signal amplitude (Y 2 ) for Ω 1 to Ω 6 frequencies are plotted in Figure 6a. According to the analytical model in Equation (9), the relation between power and square input amplitude is linear for all frequencies; however, Figure 6a demonstrates a separate pattern for Ω 1 excitation. This apparent pattern indicates that the resonant excitation frequency differed from the other frequencies.
Further analyses of the power pattern were conducted with the dimensionless study of power versus square vibration amplitude. The dimensionless study divided the parameters into the lowest-level vibration amplitude and power. The dimensionless power-vibration values for Ω 1 excitation frequency and for Ω 2 to Ω 6 excitation frequencies are plotted in Figure 6b and c, respectively.
As shown in Figure 6c, the slope of the fitted line for Ω 2 to Ω 6 harmonic excitation frequencies was approximately one; power linearly changed by the square of input acceleration at these frequencies. This conclusion is in line with the analytic model in Equation (9). In contrast, for Ω 1 frequency, where the excitation frequency was the fundamental natural frequency (Figure 6b), the slope of the fitted line was approximately three. It has been demonstrated that nonlinearity exists even in a typical no-added-tip-mass energy harvester. In many typical energy harvesters, an added mass adjusts the fundamental frequency and increases the power [30], inducing physical deformation, i.e., they become noticeably enlarged, at the acting frequency [31]. As observed in Figure 6b, increasing physical deformation creates nonlinearity, and the added-tip-mass effect is expected to enlarge nonlinearity. The tip mass affects the first vibration mode more substantially than other vibration modes (because of the larger first-mode-shape magnitude); therefore, the added-tip-mass effect on the nonlinear effects of resonant excitation is expected to be considerable.

Effects of Increasing the Vibration Level with Multi-Frequency Harmonics
As demonstrated in Section 2, practical vibration energy sources only have singlefrequency signals instead of a series of harmonic multipliers. Thus, a vibration-level study was carried out on a multi-harmonic signal using vibration-level control. The standard baseline excitation signal was the summation of harmonics with Ω to Ω frequencies, denoted by ); it excited the piezoelectric harvester, and the corresponding RMS According to the linear bending theory, the axial strain depends on the beam curvature [15], e.g., ε xx ∝ ∂ 2 w ∂x 2 . By increasing the vibration level, the strain on a piezoelectric material also increases. Assuming the matched frequency (Ω 1 = ω r=1 ), Ω 1 -frequency vibration excitation creates resonance deformation on the piezoelectric beam, and the resonance deformation is large. Therefore, a high vibration level is expected to create elevated levels of strain, which goes beyond the linear assumption between the axial strain and beam curvature. However, the experimental dimensionless study showed that for higher harmonics, e.g., Ω 2 and beyond, a high vibration level does not cause a high strain level; therefore, the linear analytical model is valid.
It has been demonstrated that nonlinearity exists even in a typical no-added-tip-mass energy harvester. In many typical energy harvesters, an added mass adjusts the fundamental frequency and increases the power [30], inducing physical deformation, i.e., they become noticeably enlarged, at the acting frequency [31]. As observed in Figure 6b, increasing physical deformation creates nonlinearity, and the added-tip-mass effect is expected to enlarge nonlinearity. The tip mass affects the first vibration mode more substantially than other vibration modes (because of the larger first-mode-shape magnitude); therefore, the added-tip-mass effect on the nonlinear effects of resonant excitation is expected to be considerable.

Effects of Increasing the Vibration Level with Multi-Frequency Harmonics
As demonstrated in Section 2, practical vibration energy sources only have singlefrequency signals instead of a series of harmonic multipliers. Thus, a vibration-level study was carried out on a multi-harmonic signal using vibration-level control. The standard baseline excitation signal was the summation of harmonics with Ω 1 to Ω 6 frequencies, denoted by Z M (t); it excited the piezoelectric harvester, and the corresponding RMS power with optimum load was P Z M . Z M (t), as given in Equation (10) where M is the signal magnification factor. Note that Ω i = i × ω r=1 , i = 1, . . . , 6. In practice, four magnification factors, M = 1, 1.2, 1.5, 1.75, and 2.0, were assigned. The associated acceleration signals from Z 1 (t) to Z M=2 (t) are shown in Figure 7b.
Micromachines 2022, 13, x FOR PEER REVIEW 12 of 17 power with optimum load was . ), as given in Equation (10) and shown in Figure   7a.
where is the signal magnification factor. Note that Ω = × ω , = 1, … , 6. In practice, four magnification factors, =1, 1.2, 1.5, 1.75, and 2.0, were assigned. The associated acceleration signals from ) to ) are shown in Figure 7b. The optimum power with 21.8 kΩ for different magnification factors is shown in Figure 8a. As expected, the power increased with the increase in the excitation vibration level; however, a nonlinear variation was observed. Further analyses were conducted on dimensionless experimental power versus square magnification factor (Figure 8b). Figure 8b shows that the experimental correlation between and was not linear, while according to the analytical modelling of ∝ , the empirical correlation between power and the magnification factor is given by  The optimum power with 21.8 kΩ for different magnification factors is shown in Figure 8a. As expected, the power increased with the increase in the excitation vibration level; however, a nonlinear variation was observed. Further analyses were conducted on dimensionless experimental power versus square magnification factor (Figure 8b). Figure 8b shows that the experimental correlation between P and M 2 was not linear, while according to the analytical modelling of P∝M 2 , the empirical correlation between power and the magnification factor is given by The output power obtained in the experiments was smaller than that obtained by employing the linear analytical theory. For explaining this statement, the interaction effects between the harmonics needed to be studied.

Interaction Effects between Different Excitation Harmonics
The interaction effects were studied by applying an excitation signal comprising six harmonics and controlling the harmonic amplitude. The excitation signal is denoted by Z(t) = α 1 sin(Ω 1 t) + α 2 sin(Ω 2 t) + α 3 sin(Ω 3 t) + α 4 sin(Ω 4 t) + α 5 sin(Ω 5 t) + α 6 sin(Ω 6 t) where α 1 to α 6 can be controllably changed. Note that Ω i = i × ω r=1 , i = 1, . . . , 6. shows that the experimental correlation between and was not linear, while according to the analytical modelling of ∝ , the empirical correlation between power and the magnification factor is given by The output power obtained in the experiments was smaller than that obtained by employing the linear analytical theory. For explaining this statement, the interaction effects between the harmonics needed to be studied. A three-level experimental design was proposed so that nonlinear relationships could also be captured and the overlap of the interaction effects could be avoided. In addition, to avoid overlapping the interaction effects, two levels for α 1 and three levels for α 2 to α 6 were employed. A fractional factorial design using the orthogonal method of resolution III was used with 36 runs of experiments [32]. Table 4 shows the orthogonal array with thirty-six runs. Moreover, three duplications were performed. Table 4. Orthogonal design of experiments on interaction effects of harmonics.  Table 5 shows the Analysis of Variance (ANOVA) table with which the meaningful variations were interpreted. Parameter F is the mean square of responses within the same treatment factor divided by the mean square of responses in all experimental runs. F > 1 indicates meaningful variation, and a greater F demonstrates a more sensitive treatment factor. The ANOVA table shows that variables α 1 to α 4 have F > 1, meaning that their influence on the output power was substantial. Table 5 shows that the influence of the α 1 factor (with Ω 1 frequency) was larger than that of the other factors and that the excitation frequencies of Ω 1 to Ω 4 had meaningful effects on the output power. This conclusion implies that harmonics in a practical vibration source up to the fourth harmonic should be considered.  Figure 9a-e demonstrate the interaction effects of different harmonics. The interaction effects of the Ω 1 -Ω 2 , Ω 1 -Ω 3 , and Ω 1 -Ω 4 harmonics were substantial, while the Ω 5 and Ω 6 harmonics had negligible effects. The charge cancellation effect of the Ω 1 -Ω 2 harmonics was visible, i.e., increasing α 2 (with Ω 2 frequency) adversely affected α 1 (with Ω 1 frequency) and increased the effect on power generation. In Ω 1 -Ω 3 and Ω 1 -Ω 4 harmonics, the interaction effects were positive. In other words, increasing α 3 (with Ω 3 frequency) or α 4 (with Ω 4 frequency) did not change the positive effect of the first harmonic.  Figure 9a-e demonstrate the interaction effects of different harmonics. The interaction effects of the Ω -Ω , Ω -Ω , and Ω -Ω harmonics were substantial, while the Ω and Ω harmonics had negligible effects. The charge cancellation effect of the Ω -Ω harmonics was visible, i.e., increasing (with 2 frequency) adversely affected (with 1 frequency) and increased the effect on power generation. In Ω -Ω and Ω -Ω harmonics, the interaction effects were positive. In other words, increasing (with 3 frequency) or (with 4 frequency) did not change the positive effect of the first harmonic. The interaction study showed that the negative interaction between the harmonic and the harmonic led to smaller power generation than when applying the linear theory, which is observed in Figure 8b. This observation can be linked to the link between the resonant harvester frequencies and the different harmonics in the excitation signal. As The interaction study showed that the negative interaction between the Ω 1 harmonic and the Ω 2 harmonic led to smaller power generation than when applying the linear theory, which is observed in Figure 8b. This observation can be linked to the link between the resonant harvester frequencies and the different harmonics in the excitation signal. As demonstrated in Table 6, the first resonant frequency and the piezoelectric-coupled resonant frequency (ω r=1 and ω pc ) were in the range of the first and second excitation harmonics (Ω i=1 and Ω i=2 ). Thus, these two excitation harmonics are expected to have considerable effects compared with the effects of high harmonics.

Proposing a Nonlinear Model
Large-deformation strain and nonlinear stress-strain constitutive equations are proposed. Base excitation a(t) deforms the piezoelectric beam with the deformation shape of w(x, t). In a nonlinear framework, the axial strain is given by at level z from the neutral axis. Moreover, the nonlinear stress-strain constitutive equations can be given by [33] T T s xx = Y s ε xx + where T xx is the axial stress, E z is the z-component electrical field, and β is the nonlinear coefficient. Thus, the internal bending moment can be calculated with Consequently, the beam stiffness (YI) and piezoelectric coupling coefficients (Y r and Λ r ) are obtained for the nonlinear model.
The electromechanical voltage equation for the nonlinear model is given by More details about the nonlinear model can be found in ref. [34].

Concluding Remarks and Future Works
This paper presents extensive experimental works aimed at practical PEH performance characterisation. Single-frequency harmonics with resonant and off-resonant excitation frequencies, and multi-frequency signals as practical vibration signals were studied as PVEH excitation signals. The investigations demonstrated the high-strain-level nonlinearities at resonance, while off-resonant harmonics followed the linear theories. Furthermore, with multi-frequency signals, resonance nonlinearity and different frequency interactions caused nonlinear performance in the PEH. This paper contributes to the realistic evaluation of piezoelectric energy harvesters. Nonlinear energy harvester models with more experimental investigation are proposed for future work. Moreover, further studies on the parameters of the flexibility characteristics of energy harvesters, such as thickness and width, are needed.