Nonlinear dynamics and triboelectric energy harvesting from a three-degree-of-freedom vibro-impact oscillator

A triboelectric energy harvester based on a three-degree-of-freedom vibro-impact oscillator is presented. Both the dynamic model of the oscillator and the theoretical model of the oscillator-based triboelectric energy harvester are established. The dynamic response and its effect on the electrical output are considered for various mass ratios and mass spacings. The study leads to the conclusions that the symmetric mass configurations of the oscillator are more beneficial to energy harvesting than the asymmetric cases. The extent of the initial spacing between the masses influences the dynamics of the system and the electrical output by triggering grazing bifurcation. High-order periodicity is found to accompany a reduction in the electrical power. An increase in mass ratio tends to increase the electrical output, and there may exist an optimal mass ratio at which the electrical output is maximized. Chatter and sticking motion can improve the output performance dramatically, while resonance, as usual, corresponds to large amplitude response, but these large amplitudes are not optimal for triboelectric energy harvesting, and thus the maximal output does not appear around resonance. This is different from other types of vibration-based energy harvesters, such as piezoelectric and magnetoelectric energy harvesters which are usually designed to operate at resonance. In addition, chatter is found to occur at low excitation frequencies, which can help harvest energy from low-frequency ambient vibration.


Mechanical energy harvesting
Among the various energy harvesting categories, harvesting energy from ambient vibration has received much attention over the past few decades. A key factor that has contributed to this is the dramatic increase in the use of sensors [1,2]. A conventional way of powering sensors is using batteries; however, batteries have to be monitored and replaced when necessary. This monitoring and replacement work can be onerous when the sensors are in remote and/or hazardous places. Energy harvesters that convert the ambient vibration energy into electricity enable the development of self-powered wireless sensors that can address these problems.
According to their basic working mechanisms, there are mainly three types of vibration-based energy harvesters: namely piezoelectric energy harvesters (PEHs), magnetoelectric energy harvesters (MEHs) and triboelectric energy harvesters (TEHs). Com-pared to PEH and MEH, TEH is relatively new and has been gaining research attention in recent years. The working mechanism of TEH involves a combination of triboelectrification and electrostatic induction through relative motion (e.g. contact-separation and sliding) between materials that have the opposite tribopolarities [1][2][3]. Moreover, TEH can work in four fundamental modes: namely, the vertical contactseparation mode, the in-plane sliding mode, the singleelectrode mode, and the free-standing triboelectriclayer mode [1][2][3]. So far, most research on TEHs is limited to experimental studies, especially involving the technology of material surface processing. This has resulted in some self-powered prototype sensors tailored to far-reaching applications including active alcohol breath analysers [4], 3D acceleration sensors [5], dynamic displacement monitoring systems [6], and wind-speed sensors [7]. The triboelectric energy harvester studied in this paper can be modified for use in a number of real applications. A notable example is to power wireless health monitoring sensors on structures such as bridges, wind turbines and undersea gas/oil pipe lines. Similar devices could also be used to harvest energy from human motion, and/or to partially replace conventional batteries so as to reduce pollution and cost [8].

Vibro-impact systems in energy harvesting
Vibro-impact systems are usually characterized by repeated vibration and impact between the components in the systems. Vibro-impact appears in many engineering applications, such as hand-held percussion machines, pile driving machines, periodic rubbing between the rotor blades and the stators in turbomachinery, braking systems in automobiles and cutting and grinding machines [9][10][11][12]. Vibro-impact systems have been studied extensively in different contexts. For example, Cao et al. [13] studied the bifurcations and penetrating rate of percussive drilling, while others have focused on suppressing the vibration of civil engineering structures by using SMAs (shape memory alloys) [14,15]. Other studies have focused on the finite codimension bifurcation and coexisting attractors in forming machines [16,17], the numerical and experimental bifurcation of a forced impacting beam [18,19], the dynamic analysis of an impact force generator in a heat exchanger [20], and vibro-impact systems under random excitation [21,22]. Interestingly, some systems purposely exploit the vibro-impact, while others try to avoid it. For example, percussion and pile driving machines employ vibro-impact motion, while turbomachinery and braking system are designed to prevent in vibro-impact conditions. However, examples of vibroimpact systems can be found in almost all types of vibration-based energy harvesters.
Vibro-impact is often used in PEH systems since impact improves harvesting efficiency through frequency up-conversion [23][24][25][26], which helps cope with the frequency incompatibility between the ambient vibration and the harvester [23]. The vibro-impact in PEH systems usually happens between flexible beams and relatively rigid stoppers. In one of the two scenarios, stoppers vibrate and trigger the oscillation of the beams [23], which are mostly cantilevers. Alternatively, the beams oscillate and hit the stoppers [27]. Either way, the low-frequency ambient vibration can be up-converted to the higher frequency vibration of the beams, which then improves harvesting efficiency. Frequency up-conversion can also be achieved without any physical impact or contact in PEH systems, for instance, through magnetic plucking [28][29][30] or using a cantilever shell structure [31]. Both of these configurations can generate an impact-like force or impulsive excitation which then increases the dominant response frequency of the cantilevers. Alternatively, some PEH designs use techniques other than frequency up-conversion to harvest energy across relatively wider frequency ranges. Examples of these techniques include the use of variable-width bistable piezoelectric cantilever beams [32] and the use of selfresonating behaviour to create a passive self-tuning PEH [33].
The most common configuration of MEH uses a tube-structure to house magnets of opposing polarity [34]. In the modelling of such systems, vibro-impact motion between the magnets is generally not considered since impact only occurs at large accelerations. However, some studies [35][36][37][38][39] have employed vibroimpact beam systems to enhance the harvesting efficiency of MEH designs using frequency up-conversion. In these cases, magnets [35,36] or coils [37][38][39] are mounted on the free end of a cantilever beam as well as on a stopper or another beam. The theoretical modelling of such vibro-impact system is mostly done by simplifying the continuous beam system into a corresponding discrete mass-spring-damper system and modelling a non-smooth impact into a bilinear or non-smooth stiffness and damping, such as in Refs. [36,37].
Among the four fundamental triboelectric energy harvesting modes, the vertical contact-separation mode and the in-plane contact-sliding mode are the two most commonly used modes. Generally speaking, the vertical contact-separation mode relies on vibroimpact. The TEHs which work on the vertical contactseparation mode take on simple structural forms, such as the most common experimental TEH system [40] which consists of two plates (assumed rigid) connected by springs. Thus far, TEH research is mostly experimental in nature and is mainly focused on the state-of-the-art surface processing technology of triboelectric materials, such as the nanopore-, nanowireor nanotube-based surface modifications [40,41]. A more complicated configuration is the triple-cantilever TEH [8,42] which may have rich dynamic responses, but the corresponding theoretical model has not yet been studied. The contact-mode freestanding TEH [43,44] seems to have higher working efficiency, and the relevant research shows some applications, such as quantitative measurements of vibration amplitude [43] and the portable power-supplying system [44].
In addition to the experimental research on TEHs, theoretical modelling of the power generation (but not the structural dynamics) of the contact-mode as well as other working modes has been explored as well. The V-Q-x relationships (where V is the voltage between electrodes, Q is the amount of transferred charges between the electrodes, and x is the separation distance in between) for contact-mode TEHs have been developed [2,[45][46][47]. The theoretical simulation for the equivalent circuit model of the TEH in the SPICE software has also been investigated [48], but structural dynamic aspects (such as impact) are not involved in the simulation. In these studies, only the constant and harmonic variations of the separation distance are used as the dynamic input to describe the relative motion between electrodes or plates, and then the equivalent capacitance and voltage are used in SPICE. However, if the dynamic input is more complicated, SPICE may not adequately model the variable capacitor and voltage source. Therefore, there is a need for more robust circuit simulators for vibro-impact systems.
Although vibro-impact systems have been employed and studied extensively both in PEH and MEH devices, a theoretical study of the dynamics of a vibro-impact TEH system has yet to be reported. Looking beyond energy harvesting, there are few reports that study three-degree-of-freedom vibro-impact systems where any two adjacent masses can come into impact. In the study of the impact among rotating shrouded blades [49][50][51], each blade is usually simplified as a cantilever beam with a tip mass (the shroud). Its impact with the adjacent shrouds is modelled using the impact between the tip mass and the two-sided weightless rigid bodies that are supported by springs. Therefore, the impact is modelled as a bilinear or non-smooth contact stiffness. A more sophisticated recent model [52] treats one active blade and two passive blades as cantilever Euler-Bernoulli beams. A tip mass representing the shroud is at each free end and any two adjacent shrouds can impact. In addition to the above, a friction wedge damper in a train suspension system was investigated both experimentally [53] and theoretically [54]. The system was simplified to a three-degreeof-freedom vibro-impact model in the latter work. Both impact (which might cause chatter) and friction (which might induce stick-slip motion) between the wedges and the bolster were considered. The theoretical model developed here is distinct from the aforementioned vibro-impact systems as it considers a three-degree-offreedom vibro-impact system in the context of TEH.
The rest of this paper is structured as follows: Sect. 2 presents the mechanical model and formulates the dynamical equations. The conservation of momentum and the coefficient of restitution are used to calculate the just-after-impact velocities. The contact forces during sticking motion are derived and applied to prevent interpenetration and keep the equilibrium of forces. Section 3 gives the theoretical electrical model of the oscillator-based triboelectric generator and shows the Simulink simulator. Then, Sect. 4 illustrates the study of the impact characteristics, the effects of initial mass spacing, mass ratio, and the use of chatter, respectively. Section 5 presents the conclusions that can be drawn in this study.

Mechanical model and dynamic equations
The configuration of the three-degree-of-freedom vibro-impact system representing the dynamic elements of a TEH is shown in Fig. 1. This model could be used to approximate a physical system in which three beams or flat plates are stacked. The three masses are initially separated by two gaps-δ 1 between m 1 and m 2 , Fig. 1 The three-degree-of-freedom vibro-impact system and δ 2 between m 2 and m 3 . The masses are connected to the sinusoidally driven base via springs and dampers.
The equations of motion away from any two consecutive impacts of the system are given as where y = A sin (ωt), and the prime denotes a derivative with respect to time, t. The separation distances, between m 1 and m 2 , and between m 2 and m 3 are represented by d 1 and d 2 , respectively. These distances are given by , and x c = δ 1 . The non-dimensional equations are then where The over-dot denotes a derivative with respect to τ , and the separation distances become When impact occurs, i.e. D 1 = 0 and/or D 2 = 0, by employing the principle of conservation of momentum and the coefficient of restitution, one can find the velocities just after the impact aṡ where μ is the coefficient of restitution, and the subscript '-' represents the state just before impact, while '+' represents the state just after impact. Sticking motion between any two adjacent masses can happen since the impact is assumed to be partially elastic. During sticking, contact force must maintain the dynamic equilibrium and prevent interpenetration between two sticking masses. Let N 12 and N 23 represent the contact forces between m 1 and m 2 , and between m 2 and m 3 , respectively. Then, the sticking conditions between m 1 and m 2 are D 1 = 0, X 1 =Ẋ 2 and N 12 > 0, while those between m 2 and m 3 are D 2 = 0,Ẋ 2 =Ẋ 3 and N 23 > 0, and the expressions of the contact forces are given as 3 Theoretical model of triboelectric generator and the Simulink simulator

Theoretical model of the oscillator-based triboelectric generator
To simulate the harvester's electrical output and assess its performance, the theoretical model of the triboelectric energy harvester based on the three-degreeof-freedom vibro-impact oscillator is established. The middle metal layer (aluminium layer) is attached to the middle mass m 2 and is assumed to have negligible mass. This layer is regarded as a freestanding layer [2,47], which works only as a charge inductor. This kind of triboelectric generator is called the contactmode freestanding triboelectric generator [1,3,47]. The working principle of the vertical contactseparation triboelectric generator involves a combination of contact electrification and electrostatic induction [1][2][3]. The charge transfer process is shown in Fig. 2. Initially, there is no charge induced, hence no electric potential difference exists between the two electrodes. As the freestanding metal layer (i.e. the middle mass m 2 ) moves upward, it comes into contact with the top dielectric (e.g. PTFE film) which is assumed massless and is attached to the top electrode (which is adhered to m 1 ). Here, charge transfer takes place at the contact area as a result of the triboelectric effect [1,2]. The aluminium tends to lose electrons, while the PTFE gains electrons in accordance with the triboelectric series [2]. This charge transfer process then results in net positive charge on the aluminium layer and equivalent net negative charge on the dielectric film (Fig. 2a). When the middle metal layer separates from the top dielectric and moves downward, an electric potential difference is simultaneously formed. The electrons in the bottom electrode are then driven to flow to the top electrode, which results in an instantaneous current that flows in the opposite direction of the electron flow (Fig. 2b). After the middle metal layer contacts with the bottom PTFE film, all induced charges are neutralized and rebalanced between those two contact surfaces (Fig. 2c). Hereafter, they separate and form the electric potential difference, and the electrons are forced to flow from the top electrode to the bottom electrode, which results in a reversed current flow (Fig. 2d). Finally, the middle metal layer comes into contact with the top PTFE film again and the cycle begins anew.
The model equation of the triboelectric generator [2,[45][46][47][48] can be given as where V is the electric potential difference between electrodes, C is the equivalent capacitance, Q is the amount of transferred charges between electrodes, and V oc is the open-circuit voltage. The equivalent capacitance and the open-circuit voltage are given as [46,47] where ε 0 is the electric constant, S is the contact area, σ is the tribo-charge surface density; and d 0 = d 1 ε r 1 + d 2 ε r 2 , where d 1 and d 2 are the thicknesses of the top dielectric and the bottom dielectric, respectively, and ε r 1 and ε r 2 are the dielectric constants of the materials of the top and bottom dielectric layers, respectively (the two dielectric layers are both assumed to be PTFE films with dielectric constants of ε r 1 = ε r 2 = 2.0). D 1 and D 2 are the separation distances between the top dielectric and the middle metal layer, and between the middle metal layer and the bottom dielectric, respectively. Assumed values for the physical parameters are S = 0.01m 2 , σ = 15 µC/m 2 and d 1 = d 2 = 0.1 × 10 −3 m. By combining Eqs. (11)-(13), one can obtain a charge transfer equation similar to the one presented as Eq. (9) in Ref. [55], the latter of which has been validated by experiment [55], providing a degree of confidence in the present modelling approach. Moreover, the results of this study agree with those in Ref. [55] in the sense that both studies conclude that higher impact or contact velocity corresponds to higher electrical output. See Sect. 4.3.2 for more discussion of this.

The Simulink simulator
The equivalent circuit model [2,46,48] for the triboelectric generator can be modelled by using the equivalent capacitance and the open-circuit voltage to represent the electrical properties of the triboelectric generator in the circuit as shown in Fig. 3.
By using Simulink software package, the output performance of the oscillator-based triboelectric generator can be simulated and analysed conveniently. The Simulink model of the equivalent circuit is shown in Fig. 4, in which the input signals, i.e. the variable capacitor and the controlled voltage source, are defined by Eqs. (12) and (13). These are combined with the dynamic responses, i.e. the separations D 1 and D 2 , of the three-degree-of-freedom vibro-impact oscillator to To assess the performance of the oscillator-based triboelectric energy harvester, the average output power per forcing cycle is introduced: where R is the resistance used in circuit and taken to be R = 100 kOhm, T is the dimensional excitation period, N is the number of dimensional excitation periods, and I (t) is the current across the resistor.

Impact characteristics
The mass distribution of the three-degree-of-freedom vibro-impact system influences the way in which the three masses interact with each other. To determine designs that work best for energy harvesting, four con-  1, 2, 3), if not otherwise specified. The fourth-order Runge-Kutta method is used to solve the differential equations. The tolerances used to detect impact or sticking between two masses are both set to 10 −6 .

Different impact characteristics
It has been found that different mass distributions can lead to different impact positions. In symmetric case  Fig. 5a. Interestingly, the opposite is true in symmetric case (2), as shown in Fig. 5b. In the third configuration, impact only appears below the static equilibrium of each mass, even for m 2 (see Fig. 5c), and it turns out to be the opposite in case (4), as shown in Fig. 5d. From the energy harvesting perspective, it is preferable for m 2 to encounter impact both above and below its static equilibrium position. Therefore, the two symmetric configurations are of greater interest and will be compared in the following study.

Comparison between the two symmetric configurations
It is not easy to make general observations regarding the two symmetric configurations since altering some of the specified parameters (e.g. the mass ratio or excitation frequency) might lead to different conclusions.
In an attempt to make a broad comparison, both low and high mass ratios of each configuration are considered, and the average power output of each configuration at each mass ratio is calculated both at lower and higher excitation frequencies.
In the first symmetric configuration, the mass distribution is configured to m 2 = 120×10 The average power P a versus the number of excitation periods N p is calculated for different excitation frequencies and mass ratios for both configurations. The results are shown in Fig. 6. It can be seen that the first symmetric configuration has overall better electrical output performance than the second. Therefore, the first symmetric configuration will be studied further.

The effect of mass spacing
The initial mass spacing, δ 1 and δ 2 , could be symmetric or asymmetric, and the choice will influence the dynamics of the system deeply because the spacing determines the impact criteria. Additionally, the mass

Asymmetric mass spacing
Let R δ = δ 2 δ 1 be the mass spacing ratio in the asymmetric spacing case. The bifurcation diagram of X 2 versus R δ and its local enlargement are calculated and shown in Fig. 7a, b, respectively, while Fig. 7c, d show the corresponding "bifurcation diagrams of the velocities just prior to impact", in which the values of the quantities concerned, for example, relative velocities between m 1 and m 2 , (i.e. V 12 ), and between m 2 and m 3 , (i.e. V 23 ), are only sampled at impact events within each forc-ing cycle. This is one way of displaying the dynamic behaviour of vibro-impact and the present bifurcation idea is similar to that in Ref. [56], which is quite different from the conventional bifurcation concept. Therefore, the countable dots at each value of R δ represent the number of impacts, rather than the order of periodicity, in one vibration period of periodic motions, while the uncountable dots indicate the impacts of non-periodic vibrations. The parameter values used in simulation are δ 1 = 20 mm, = 0.70, m 1 = m 3 = 90 × 10 −3 kg and m 2 = 120 × 10 −3 kg.
From the bifurcation diagrams of X 2 , i.e. Fig. 7a, b, it could be seen that the system undertakes complicated dynamic variation with the change of R δ . The corresponding bifurcation diagrams of the just before impact Grazing bifurcation occurs at R δ = 1.419, and it is obvious in Fig. 7b that the dynamics of the system becomes much more complicated due to grazing bifur-  Fig. 9 to verify and show the details, and it is clear that grazing impact appears only between m 1 and m 2 . Besides, the vibration responses are asymmetric owing to the asymmetric mass spacing.
Since there occurs period-doubling bifurcation with non-periodic windows in Fig. 7b, and both periodic and non-periodic motions appear with the variation of R δ , a comparison of the electrical output among different motion characteristics is convenient. The PSD (power spectral density) and Poincaré section of X 2 , the average harvested power versus the number of excitation periods are shown in Fig. 10 for different types of response. The average power output for these different responses is shown in Fig. 11. It is concluded that high-order periodic motion has lower average power output relative to simple periodic motion. It appears that power decreases as the order of the periodicity increases. However, none of them seem to generate more than a few micro-watts of power.

Symmetric mass spacing
In the consideration of symmetric mass spacing, the structural parameters (i.e. mass, stiffness, damping, mass spacing) are also all assumed symmetric. Consequently, the dynamic response of the system is mostly symmetric, though symmetry-breaking [57] of the dynamic response can exist.
In the modelling of the non-dimensional system, the mass spacing, δ 1 , is used as a normalization parameter (i.e. the non-dimensional displacements are X i = 1, 2, 3)). Therefore, for a given dimensional displacement, a smaller mass spacing would result in a larger non-dimensional displacement. To circumvent this dependence, the bifurcation diagram in Fig. 12 is shown in terms of δ X 2 versus δ rather than X 2 versus δ. Here, the excitation frequency is held fixed at = 0.75. The diagram shows that vibration amplitude tends to grow with an increase in mass spacing and that the response can be complicated at smaller mass spacing. It is noted that there exists Hopf bifurcation, which might be of great importance in some dynamic systems, such as in railway vehicle systems, where Hopf bifurcation can lead to unstable motions, such as hunting [58].

The effect of mass ratio
While the size of the mass spacing mostly affects the location at which the masses come into contact, the mass ratio can influence the way the masses interact with each other. If the mass ratio equals one, i.e. R m = m 2 m 1 = m 2 m 3 = 1, all three masses are equal, there is no impact among them after the system reaches its steady state (all other parameters, i.e. the mass spacings, the stiffness and the dampings, are assumed to be  Fig. 11 The comparison of P a between different types of motions Fig. 12 The bifurcation diagram of δ X 2 versus the symmetric mass spacing δ symmetric as well). This is due to the absolute symmetry of the structural parameters.

Mass ratio's effect on the number of impacts
To observe the dynamic response and electric output at mass ratios other than one, the time histories of D 1 and the generated current I and the diagram of the average power P a are shown in Fig. 13 for different mass ratios. It is be observed that the number of impacts between m 1 and m 2 in one excitation period increases from one to three as the mass ratio R m increases. This results in an increase in electrical output (note that impacts between m 2 and m 3 show the same effect).

Optimal mass ratio
The bifurcation diagram of X 2 versus R m when = 0.72 is shown in Fig. 14a, while Fig. 14b shows the cor-responding "bifurcation diagram of V 12 against R m ". Here, it can be seen that there is an optimal mass ratio at which the impact velocities reach a maximum at each impact. In Fig. 14b, the number of the countable dots at each R m represents the number of impacts that happen in one excitation period for periodic motions; for instance, there are three impacts between m 1 and m 2 in one forcing cycle of the periodic motion at R m = 4, while the number of the 'line segments' at each R m (of the values approximately between R m = 3.06 and R m = 3.64) corresponds to the number of impacts that happen in one excitation period for non-periodic motions; for example, three impacts happen between m 1 and m 2 at the non-periodic motion of R m = 3.5. Because higher impact velocity corresponds to higher electrical output, there exists an optimal R m at which the system has the best electrical output performance. Figure 15 shows the average power outputs of three different mass ratios for comparison and verification. It is clear that the average power output reaches the maximum around R m = 3.5. The parameters used are the same as before except δ 1 = δ 2 = 20 mm and m 1 = m 3 = 50×10 −3 kg. It is again noted that a higher impact or contact velocity induces a higher electrical output agrees with the conclusion drawn in Ref. [55].
The optimal mass ratios, R mo , of the frequencies between = 0.55 and = 0.95 are obtained and shown in Fig. 16, and it is clear that the optimal mass ratio decreases with the increase of excitation frequency. The fitted polynomial is found to be R mo = 11.43 2 − 28.46 + 18.03, ∈ [0.55 0.95]; However, there would be no impact owing to the absolute symmetry of structural parameters as R mo gradually approaches one. The bifurcation diagrams of X 2 and V 12 of two higher excitation frequencies, = 1.30 and = 2.00, are calculated and given in Fig. 17 as examples to show the variation of impact velocity at higher excitation frequencies. It can be observed that the magnitude of impact velocity would gradually approach a stable value as the mass ratio grows bigger, which means there is not an obvious optimal mass ratio for higher excitation frequencies.
For lower excitation frequencies, the response of the system can become quite complicated with the variation of mass ratio. The bifurcation diagram of X 2 versus R m for = 0.30 and its local enlargement are given in Fig. 18 as an example, and there appears a period-doubling cascade with non-periodic windows in Fig. 18b. Figure 19 shows the Poincaré sections of

Use of chatter
Chatter is a challenge in machining processes including turning, milling and drilling. It can cause problems like poor surface finish, excessive noise, breakage of machine tools, decreased tool life, and low productivity [59,60]. Chatter is characterized by many consecutive impacts and might therefore be exploited in triboelectric energy harvesting from vibro-impact systems. In the study of a multi-degree-of-freedom vibroimpact system, Wagg [56,61] discussed the chattering and sticking behaviour which appeared at low forcing frequencies, and showed that a variety of non-smooth events could possibly occur. Here, a numerical simulation is conducted to determine whether chatter can appear in the present system at low forcing frequencies.
The bifurcation diagram of X 2 versus , which includes a chatter range, is shown in Fig. 20. The parameters used are m 1 = 50 × 10 −3 kg, R m = 3 and δ = 5 mm. To observe the details of the chatter range as well as the resonant range, local views are shown in Fig. 21(a1), (b1), respectively, while Fig. 21(a2), (b2) gives the relevant "bifurcation diagrams of V 12 ", in which the number of the dots at each represents the number of impacts that happen in one excitation period, and the y-values of the dots represent the relative impact velocities. It is clear that there are plenty of consecutive impacts when chatter takes place while only three impacts exist at the excitation frequencies around resonance. Besides, chatter does not appear around a resonance but at some lower frequencies. Resonance, as usual, can have the strongest vibration and impact, but its high vibration amplitude would results in large separation distances between plates, which might not be beneficial to triboelectric energy harvesting. Therefore, there is the possibility of using chatter to harvest more energy than can be harvested by operating at resonance. It can be seen, from Fig. 22a2, that the middle mass m 2 encounters chatter on both sides. Hence, there are many impact events. However, only six impacts in total would occur in the resonant case in Fig. 22b. On the other hand, most of the impacts in chatter are quite small compared with those at resonance. Therefore, to determine which one is superior in terms of harvesting energy, the average power output both around chatter zone and resonant zone are given in Fig. 23 to compare the capability of efficient energy harvesting, in which it is clear that chatter frequencies perform much better than resonant frequencies in producing higher electrical output. Besides, since the chatter zone appears in a relatively lower frequency range compared with the resonant zone under the present parameters, it will be quite beneficial to harvesting energy from lowfrequency ambient vibrations.
To consider the details of both chatter and resonance cases and understand why forcing the system in the chatter regime results in such relatively high power output, some more detailed results of = 0.600 and = 0.765 are shown in Figs. 24 and 25. The nondimensional time and dimensional time in these figures correspond to each other. It can be seen that the equivalent capacitance gets bigger when (D 1 + D 2 ) becomes smaller, and the peaks of the equivalent capacitance correspond to the peaks of current. Although the non-chatter case has higher equivalent voltage, i.e. at = 0.765, the equivalent capacitance in the chatter case is much bigger, actually two orders of magnitude bigger, than in the non-chatter case. The relationship between the equivalent capacitance and the generated current indicates that the equivalent capaci- Fig. 21 The local enlargements of Fig.18 for (a1) chatter range and (b1) resonant range, and their corresponding bifurcation diagrams of V 12 , i.e. (a2) and (b2), respectively tance's variation seems to play a bigger role than that of the equivalent voltage. Besides, from Fig. 24a, b, it can be observed that the two chatter events, respectively, between m 1 and m 2 , and m 2 and m 3 appear just one after the other, which means that the middle mass m 2 gets into chatter motion with one of the side masses, i.e. m 1 or m 3 , firstly and then gets into chatter motion with the other one. More specifically, the relative motion between the middle mass m 2 and one of the side masses firstly undergoes chattering until they split apart from sticking, just after which the other side mass impacts the middle mass (Fig. 24d, e) pushing it into an additional impact with the mass from which it had just separated. After this, the original side mass flies away, while the newcomer (the other side mass) begins chattering with the middle mass. This is similar to the motion of a Newton's cradle. Interestingly, it seems that the small impact events during chatter are not directly responsible for the spikes of the electric output. Rather, it is the chatter-induced sticking motion that provides a chance for these three masses to come into impact nearly at the same time. When this occurs, the capacitance between the top and bottom electrodes spikes up because of the great decrease in the separation distance between them (the capacitance between two electrodes is inversely proportional to the separation distance in between, see Eq. (12)).

Conclusions
A dynamic model of a three-degree-of-freedom vibroimpact oscillator is developed. In this system, every mass can impact its neighbour(s), which produces very complicated vibro-impact motions. Then the theoretical model of the oscillator-based triboelectric energy harvester is constructed, and the equivalent circuit (1) Different mass distributions could lead to very distinct impact characteristics. Among those four possible configurations, i.e. two symmetric configurations: m 1 = m 3 < m 2 and m 1 = m 3 > m 2 , and two asymmetric configurations: m 1 > m 2 > m 3 and m 1 < m 2 < m 3 , symmetric cases tend to be more beneficial to energy harvesting than asym-metric cases, and the first symmetric configuration has better harvesting performance than the second. (2) The initial mass spacing between the three masses affect the displacement conditions of impact directly, thus grazing bifurcation or impact is likely to appear, which could influence the dynamics of the system and the electrical output deeply. Besides, it has been found that period-doubling motion could somewhat reduce the electrical output, and the higher the order of the periodicity, the more the reduction.
The mass ratio determines how the masses impact and interact with each other. The number of impacts in one excitation period increases with the growth of the mass ratio, which then leads to the improvement of electrical output. There may also exist an optimal mass ratio at which the electrical output reaches the maximum, and it decreases with the increase in excitation frequency. (4) The electrical output increases greatly when chatter and sticking occur in motion. Resonance, as usual, can have the highest vibration amplitudes and impact velocities, but this can be detrimental to electrical output and thus the maximal electric output does not appear at around resonance, which is different from other types of vibration-based energy harvesters, such as piezoelectric energy harvesters and magnetoelectric energy harvesters, both of which usually have the highest electrical output at resonance. In addition, the appearance of chatter and sticking in the low excitation frequency range can help harvest energy from low-frequency ambient vibration, which should be leveraged in the design of future systems.