Nonlinear structural dynamics of a new sliding-mode triboelectric energy harvester with multistability

A new sliding-mode triboelectric energy harvester in the form of a cantilever beam with a tip mass that is acted upon by both magnetic and friction forces is modelled and simulated. A numerical scheme based on the trapezoidal rule with the second-order backward difference formula (TR-BDF2) method is introduced to solve the combined non-smooth mechanical and stiff electrical system. This is the first study of the structural dynamics of the sliding-mode triboelectric energy harvesting; additionally, a magnetic field that induces multistability is present. A comparison between the coupled and uncoupled electromechanical models suggests that the electrostatic force between the electrodes can be ignored, which makes the uncoupled model preferable in the dynamical analysis. The influence of the non-conservative force (the friction force) on the multistability of the system is investigated. It is found that the distribution of the multistability on the parametric plane changes even when a small amount of friction is involved, and the areas of bistability and tristability shrink while that of the monostability expands. A comparison among these three types of stability reveals the superiority of invoking bistability as it facilitates broadband energy harvesting. The excitation level plays an important role in inducing the snap-through motion (the interwell oscillation) by enabling the crossing of the energy barriers between wells. The increase in the friction shrinks the frequency band of interwell oscillations from high frequencies down to low frequencies on the discrete frequency sweep. An analysis of the basins of attraction finds that at low frequencies the bistable system can undergo only interwell oscillations, while the tristable system can merely experience intrawell oscillations. The basins can intermingle with each other in both bistable and tristable systems. Finally, an increase in the excitation level can break the basins into discrete pieces and/or points.


Triboelectric energy harvesting
Harvesting vibrational energy from ambient vibration and converting it into useful electric energy is a topic that has received substantial attention. Piezoelectric [1] and electromagnetic [2] energy harvesters are the two most common types and have been studied extensively. The triboelectric nanogenerator (TENG), which was first introduced in 2012 [3], offers a new fashion of energy transduction. It works through a combination of contact electrification and electrostatic induction [4], both of which happen between two triboelectric materials of opposite tribo-polarities. (e.g., PTFE and aluminium) [4]. During interaction, materials of positive tribo-polarities tend to lose electrons while materials of negative tribo-polarities tend to gain electrons, and thus an electrical current.
Compared with piezoelectric and electromagnetic energy harvesters, triboelectric energy harvesters (TEHs) have several advantages: They can be constructed from a wide range of common materials that are much cheaper than piezoelectric crystals and ceramics and neodymium magnets, TEHs have a high conversion efficiency, are easy to fabricate, have widespread applications, and are cost-effective [5]. Moreover, the outputs of TEHs are usually high in voltage but low in current, while piezoelectric and electromagnetic energy harvesters are the opposite, which makes TEHs a better choice when the targeted sensors in application need high-voltage inputs. In addition, TEHs have four different working modes which offer more freedom when designing TEHs for a specific application, and they are the vertical contactseparation mode [6], the in-plane sliding mode [7], the single-electrode mode [8], and the freestanding triboelectric-layer mode [9]. Although the four modes differ from each other, the electrical outputs of TEHs based on them can all be described using the V-Qx relationship [4], where V, Q, and x are the voltage, the transferred charge, and the separation distance between electrodes, respectively.
A variety of TENG-based prototypes have been developed for practical applications such as, triboelectric motion sensing [10], human health monitoring [11], wind speed sensing [12], plasmonic ultraviolet detection [13], and virtual reality 3D-control sensing [14]. In these devices, the triboelectric materials are usually etched to have micro-or nanoscale surface patterns including nanowire arrays [7], nanopores [15], or pyramid and cubic patterns [16]. Patterned triboelectric materials have been demonstrated to be much more efficient than non-patterned ones [5]. In addition to surface modification, an electron blocking layer-based interfacial design has been found to dramatically increase the output of TENGs [17]. Moreover, instead of modifying a material's surface to improve its performance, some new triboelectric materials have been developed, such as a friction material [18] and a flexible fibre material [19].
Despite the development of diverse TENGs, research into their structural dynamics is rare, and a study investigating sliding-mode TEHs is still missing. Structural vibration is closely related to harvesting performance and thus should be considered. Piecewise stiffness may be introduced to widen the effective frequency bandwidth of TEHs that use the vertical contact-separation mode [20]. Investigations of vibroimpact dynamics benefit the design of harvesters involving repeated impact (or contact) and separation, such as in a contact freestanding mode TEHs [21]. Further, a velocity-dependent coefficient of restitution can be identified experimentally and used in the modelling of impact to better approximate a physical device [22]. The electrostatic force between two charged electrodes may be too weak to affect the structural dynamics and thus may be ignored in modelling, which then leads to decoupling of the mechanical system from the electrical system [22]. The frequency-shifting phenomenon, which often appears in piezoelectric energy harvesting due to the resistive shunt damping effect [1], may disappear in TEHs [22]. It has also been found that a TR-BDF2 numerical integration scheme can be effectively used when the electrical system is stiff [22].
In microelectromechanical systems (MEMS), electrostatic transducers (or generators) have been studied extensively [23]. However, the necessity of using a charge source along with switching losses reflects the limitations of the standard (electret-free) electrostatic generators [24]. Electret-based electrostatic generators can eliminate the need for an initial charge source and, therefore, are now prevalent [24,25]. Their working mechanisms are similar to the TEHs in the vertical contact-separation mode or the in-plane sliding mode, but usually without any contact between the electret and the counter-electrode [25]. The main difference is that a variable open-circuit electric potential difference (V oc ) is used in TEHs [4] while a fixed surface voltage (V) is used in electret-based electrostatic energy harvesters [26]. In terms of the topology of the harvesters themselves, most of the electrostatic converters are miniaturized and combshaped with designs that are usually derived from accelerometers [23][24][25][26], whereas TEHs can adopt various shapes and dimensions [4].

Bistable energy harvesting
Bistability is often exploited in piezoelectric energy harvesters, though it can also be found in electromagnetic energy harvesters [27][28][29]. However, the application and study of bistability in triboelectric energy harvesting was only very recently reported in the context of a tuneable shock sensor [30]. Bistable energy harvesters can be categorized into three different mechanisms according to their manifestations of bistability [31]: magnetic attraction bistability [32][33][34], magnetic repulsion bistability [28,35], and mechanical bistability [30,36,37]. The former two mechanisms rely on a magnetic force while the latter depends on either a mechanical pre-load [30,36] or material anisotropy [37].
Bistable systems are well known for having a double-well potential energy curve which results in two distinct types of oscillation: low-energy intrawell and high-energy interwell oscillations, the latter of which has been found to be more efficient in energy harvesting applications [31]. Bistability, as a particular case of multistability, has been extensively studied in the context of symmetries, phase transitions, and hysteresis [38]. Numerous attractors can coexist for a fixed set of parameters of a bistable system. For instance, a small change in initial conditions can turn a steady-state intrawell oscillation into an interwell oscillation owing to their coexistence [39]. For such systems, basins of attraction are often investigated, such as in Refs. [40,41]. In addition, bistability has been shown to broaden the frequency range from which energy can be extracted, and to provide a highenergy orbit over a wide range of frequencies [42]two features that are highly desirable in energy harvesting applications. In addition to numerical and experimental studies, some approximate analytical methods have been used to investigate bistable harvesters, such as harmonic balance method [43], the method of multiple scales [44], and Melnikov's method [45]. By analysing the bifurcation and stability characteristics of the system, design and operation guidance can be provided. For example, the frequency span between a saddle-node and a period-doubling bifurcation might be used to approximate the bandwidth within which the harvester can achieve highenergy orbits and thus produce high-energy outputs [46]. In addition, multistability rather than bistability has also been utilized in energy harvesting, such as in Refs. [47,48]. Although various nonlinearities have been exploited in vibration energy harvesters, it should be noted that there still exists a notable knowledge gap in integrating nonlinear energy harvesters with effective nonlinear rectifying and power management circuits for practical applications [49], because using energy harvesters as the sources is quite different from using batteries [50]. However, this paper will not target on bridging this knowledge gap.

Aim and objectives of this study
A new triboelectric energy harvester that works in a sliding mode and includes bistability is presented. This study aims to investigate as well as optimize it from a structural dynamics perspective. The following contributions concerning TEHs are made: 1. A sliding-mode triboelectric energy harvester is modelled and studied from a structural dynamics perspective. 2. Both the coupled and uncoupled electromechanical models are established and compared to investigate the effect of electrostatic force between the electrodes. 3. The different types of stability resulting from a magnetic field are identified and compared for a TEH device. 4. The effect of friction in triboelectric energy harvesting is studied in the context of structural dynamics.
The outline for the rest of this paper is as follows: Section 2 introduces the design of the harvester and describes its working mechanism. Section 3 presents the mechanical modelling of the vibration of a cantilever beam when its tip is subject to both magnetic and frictional forces. The electrical modelling of the harvester's output is given in Sect. 4. The method used to solve the non-smooth mechanical dynamical system and the stiff electrical system is described in Sect. 5.1. Section 5.2 analyses the effect of the electrostatic force. Coupled and uncoupled electromechanical models are compared. The influence of the non-conservative force (the friction force) on the multistability of the system is investigated in Sect. 5.3. Section 5.4 makes the comparison among three types of stability of the system via discrete numerical frequency sweeps. Section 5.5 presents a study investigating the effect of the excitation level on the system response. The influence of the friction is discussed in Sect. 5.6. Afterwards, the basins of attraction are analysed in Sect. 5.7. The conclusions of this study are drawn in Sect. 6.

Design
The configuration of the proposed sliding-mode triboelectric energy harvester is shown in Fig. 1a. All the components are integrated onto a base which receives a sinusoidal excitation. A cantilever beam is clamped onto the base and its free end is attached with a slider in which a cylindrical magnet is embedded. Additionally, two identical magnets are fixed to the base and symmetrically located to attract the oscillating magnet. As shown in Fig. 1b, the slider has a circular contact area and slides over two specially made patches. The slider is attached with an electrode (known as the top electrode) and a dielectric (such as PTFE). The base is composed by an electrically active patch (EAP, aluminium) and an electrically inert patch (EIP, PTEG). The circular electrically active area works as the bottom electrode and has the same circular area as the top electrode. Both the top and bottom electrodes are wired into an outside circuit to form a loop (not shown in Fig. 1). Under the base excitation, the slider can oscillate, and its motion will be affected by the magnetic field.
During sliding, triboelectrification only happens in the contact area between the slider and the electrically active patch (EAP). The charge transfer process [51] of the TEH is depicted in Fig. 2, where a resistor denoted by R is wired with the harvester and the current in the circuit is represented by I (note that the electrically inert patch (EIP) is not shown in Fig. 2). When the top dielectric fully overlaps the bottom electrode (the EAP), charge transfer occurs at the contact area owing to the triboelectric effect, and the aluminium loses electrons, while the PTFE gains electrons. The charge transfer process at this stage will result in net positive charges on the EAP and equivalent net negative charges on the dielectric, as shown in Fig. 2a. As the slider slides out of the EAP, an electric potential difference between the two electrodes is established simultaneously, which drives the electrons in the top electrode to flow to the bottom electrode. The electron flow induces an instantaneous current that flows in an opposite direction as shown in Fig. 2b. When the slider fully slides outside the EAP, the electrons will be balanced between the top dielectric and top electrode, as shown in Fig. 2c. Once the slider starts sliding back, as illustrated in Fig. 2d, the balance gets broken and the current will flow from the top electrode to the bottom electrode until the slider fully overlaps the EAP again. Thereafter, the process is repeated as slider oscillates. Here, a relatively simple or a purely resistive circuit rather than a more complicated circuit is used, because this study is mainly focused on the structural dynamics rather than the electrical dynamics of the proposed vibration energy harvester. This is commonly seen in the study of vibration energy harvesters from the perspective of structural dynamics. Nevertheless, the use of more complicated circuits can also be found in the literature. In the study of a piezoelectric energy harvester, a resistor-inductor resonant circuit, which is  Fig. 1 Configurations of a the sliding-mode triboelectric energy harvester and b the slider on patches in serial connection with the harvester, has been demonstrated to be able to enhance the energy harvesting efficiency [52]. Thus, it would be interesting to investigate the influence of such circuits on the performance of triboelectric energy harvesters, and this might be carried out in a future study.
Note that although the configuration is similar to some bistable piezoelectric energy harvesters [39], the energy harvesting mechanism is completely new and the system is no longer smooth due to the involvement of friction between the slider and the patches. Nevertheless, the presented triboelectric energy harvester can be combined with a bistable piezoelectric energy harvester to form a hybrid energy harvester, in which case the harvesting performance may be enhanced. As for the practical applications, the proposed vibration energy harvester can be used to power sensors for structural health monitoring, such as those installed on bridges, wind turbines, and oil pipelines. The integrated network of the harvesters and sensors can be wireless as well as self-powered, which no longer needs the use of conventional batteries and thus avoids the associated motoring, replacement and recycling costs.

Modelling of the mechanical system
The equations modelling the mechanical system can be derived using the Lagrange-d ' Alembert principle which is an extended Hamilton's principle in the presence of external forces, such as friction force. The corresponding variational equation is where T and U are the kinetic and the potential energies, W is the work done by the external force, i.e., the friction force (kinetic) F f , and W ¼ Þ is the Dirac's delta function and w rel x; t ð Þ is the transverse displacement relative to the base, d represents an infinitesimal variation, and t 1 and t 2 are the start and the end times.
The kinetic energy of the cantilever beam can be given as where q b and A b are the cantilever beam's density and cross-sectional area, and gðtÞ is the base displacement. The kinetic energy of the slider is where m s and J s are the mass and the second moment of inertia of the slider, where the latter is considered negligible for the small tip mass used in this harvester. The elastic potential energy of the cantilever beam is where E is the Young's modulus of the beam and I is the second moment of area of the beam's cross section. The relative locations of the magnets are shown in Fig. 3. Magnetic dipoles are used to represent them and the modelling work below follows Refs [42,53].
---------+ + + + + + + + + + + + + + + + + + ------------------+ + + + + + + + + + + + + + Attached to slider EAP Bottom electrode (Al) Fig. 2 The charge transfer process of the presented sliding-mode triboelectric energy harvester The elastic potential energy provided by the magnetic field is given as where n is the number of fixed magnets (two in this case), B i is the magnetic field generated by the ith fixed magnet at the location of the tip magnet (in the slider), l s is the magnetic moment vector of the tip magnet and is expressed as where M s is the magnitude of the magnetization of the magnet at the cantilever tip and can be estimated using the residual flux density B r as M s ¼ B r =l 0 , where l 0 is the vacuum permeability; V s is the magnet's volume, h is the slope of the beam at the tip,ê x andê y are the unit vectors along the x-and the y-axis, respectively. Similarly, the magnetic moment vectors for the fixed magnets are given as where M i and V i are the magnitude of the magnetization and the volume of the ith fixed magnet, respectively. The separation vectors from the positions of the fixed magnets to that of the tip magnet are expressed as where d x is the longitudinal displacement of the tip and it can be derived as Then, the magnetic field generated by the ith fixed magnet upon the tip magnet is of the form where Á k k 2 denotes the Euclidean norm and r represents the gradient operator.
Assuming the cantilever beam is a uniform Euler-Bernoulli beam and supposing the first mode is dominant, the transverse displacement of the beam can be written in the form where q t ð Þ is the temporal modal coordinate, /ðxÞ is the mode shape function where r and b are determined from the associated eigensystem, C can be derived from the corresponding orthonormality conditions (the corresponding formulas are not given here). Substituting Eq. (11) into Eqs.
(2) to (5), one can get the expressions of kinetic energies of the beam and the slider and potential energies of the beam and the magnetic field in the form of the modal coordinate and mode shape function as where the over-dot and the prime denote a derivative with respect to time, t, and location, x, respectively, r i are the Euclidean norms of vectors r i , and cosh and sinh can be obtained through geometrical approximation. These terms are given by The Lagrangian of the mechanical system is expressed as Based on the Lagrange-d ' Alembert principle, which is equivalent to the Euler-Lagrange equations with external forces, the governing equation of the mechanical system can be obtained as d dt where F f is the kinetic friction force exerting on the slider and will be detailed in the following section. Incorporating Eqs. (11), (13) and (14) into Eq. (15) and applying the orthonormality conditions and including damping yield where f is the damping ratio, x is the undamped natural frequency of the first mode of the beam without the tip mass, € g ¼Âsinu t ð Þ is the base excitation, whereÂ is the acceleration amplitude and _ u t ð Þ=2p ¼ f is the excitation frequency in Hz, and F m is the magnetic force exerted by the magnetic field at the tip. The expression for F m and those of a and x are given as Owing to the rotation of the cantilever tip mass, a circular contact area between the slider and the patches is more convenient for modelling the charge transfer process and calculating the friction force in between. As shown in Fig. 1b, w rel is the displacement of the slider, A 1 and A 2 are the contact areas between the slider and the EAP, and between the slider and the EIP, respectively. The slider has the same circular area as the EAP, and its radius is r and hence According to Coulomb's law of friction, the kinetic friction force on the slider can be given as where l k1 and l k2 are the kinetic coefficients of friction between the slider and the EAP and between the slider and the EIP, respectively; sgn Á ð Þ is the signum function, A is the total contact area, and N is the normal force are expressed as where C r 2 0 1 ½ is the contact force ratio describing the contact situation between the slider and patches and g is the gravitational acceleration.
Without loss of generality, stick-slip may happen between the slider and the patches during oscillation. During sticking, the beam slider system is in a dynamically balanced state and the slider has no motion relative to the moving base. The static friction force between the slider and the patches is given by where m e is the equivalent mass of the cantilever beam and can be expressed as m e ¼ 33 140 q b A b L b . The maximal static friction force is as follows where l s1 and l s2 are the static coefficients of friction between the slider and the EAP and between the slider and the EIP, respectively. The conditions for sticking are then 4 Modelling of the electrical system According to Ref. [51], during slipping and when w rel L b ; t ð Þ j j \2r, the equivalent capacitance, i.e., C e , and the open-circuit voltage, i.e., V oc , for the presented triboelectric energy harvester can be derived as where e 0 and e r are the vacuum permittivity and the relative permittivity of the top dielectric (PTFE), respectively, t d is the thickness of the top dielectric, and r is the tribo-charge surface density. The relationship between the voltage across a resistor R in an outer circuit, i.e., V, and the amount of transferred charges between electrodes, i.e., Q, of the harvester can be given by [4,51] Applying Ohm's law, i.e., V ¼ RI ¼ R dQ dt , and substituting Eqs. (21) and (22) into Eq. (23) yield the differential equation However, a singularity will occur in Eq. (24) when the top electrode (attached to the slider) fully slides out of the bottom one (the EAP), i.e., the denominators of the coefficients of the second and the third terms become zero when w rel L b ; t ð Þ!2r, and it will result in a zero equivalent capacitance and an infinitely large open-circuit voltage. Nevertheless, the transferred charges between the electrodes will saturate when the top electrode fully slides out of the bottom one [7,51], and the open-circuit voltage V oc has been measured to be finite in experiment when the two electrodes fully slide out of each other [7]. To avoid this singularity, it is assumed that the top electrode will not fully slide out of the bottom electrode during oscillation. This will be monitored during the numerical simulations.
When sticking, the transferred charges between electrodes can be given by where q s is the temporal modal coordinate when sticking motion starts to take place and t s is the corresponding time instant.

Numerical scheme
Since the order of magnitude of the vacuum permittivity e 0 is -12 and the load resistance R can vary in a wide range, the ODE of the electrical output, i.e., Equation (24), is likely to be stiff. The TR-BDF2 (trapezoidal-backward differentiation formula of order two) method, which is an implicit single-step method of second-order accuracy, is prevalent in circuit and semiconductor simulations and has demonstrated good performance among various one-step methods [54][55][56][57] and, therefore, is used to solve Eq. (24). The method is composed of two substeps: a trapezoidal step from t n to t nþc and a second-order backward difference step from t nþc to t nþ1 , where t nþc ¼ t n þ cDt n , c 2 ð0; 1Þ and t nþ1 ¼ t n þ Dt n , and Dt n is the time step. For the integration of a smooth system described by dq dt ¼ hðqÞ, the solution is firstly advanced from t n to t nþc using the trapezoidal rule as In the second sub-step, the solution is then advanced from t nþc to t nþ1 using the second-order backwards difference rule as The choice of c can affect the solution in many ways [54,58], and it is suggested to take c ¼ 2 À ffiffi ffi 2 p because it results in the least truncation error [54], the same Jacobian matrix for both sub-steps (the Newton-Raphson method is used to solve the implicit difference equations), and the largest linearized stability region [58].
However, owing to the non-smoothness of the mechanical system, i.e., the sign change of the friction force and the stick-slip motion, the TR-BDF2 method can only be directly applied to steps without any nonsmooth events. For a time step containing non-smooth events, only the trapezoidal rule or the TR method is used for the whole step. The derivation of the corresponding difference equations is given in ''Appendix A''.
The values of the main parameters used in simulation are given in Table 1.

The effect of the electrostatic force
Theoretically, there exists an attractive electrostatic force between the slider and the EAP because the dielectric on the slider, i.e., the PTFE film, and the EAP (Al) are oppositely charged. The electrostatic force can be expressed as And the voltage across a capacitor, i.e., V c , can be given by Substituting Eqs. (21) and (29) into Eq. (28), the electrostatic force can be rewritten as Since the harvester is of the sliding mode, the electrostatic force will act as part of the normal force exerted on the slider, and the new normal force is then The electrostatic force, therefore, may affect the system's behaviour through friction, which is thus worth investigating. It is noted that the mechanical system described by Eq. (16) and the electrical system modelled by Eq. (24) are coupled when the electrostatic force is considered, in which case, the Jacobian matrix should be updated accordingly, and the new matrix elements are given in the appendix.
To simplify the comparison, the two fixed magnets are removed and the slider no longer experiences the magnetic force. The comparison of the tip displacements and the RMS voltage outputs under discrete numerical frequency sweeps (from 5 to 15 Hz with a step of 0.05 Hz) with and without the consideration of the electrostatic force (EF), is shown in Fig. 4. If the EF is adequately big, the model that includes the EF should produce a smaller amplitude response relative to the model without the EF because of the increase in the friction force caused by the EF. However, it can be seen in Fig. 4 that both the mechanical and electrical responses are almost identical in the two cases. Therefore, it can be concluded that the electrostatic force between the two electrodes (or between the slider and the EAP) for this TEH can be neglected and the EF will be not be modelled in the subsequent results. Nevertheless, the feedback of the EF into the electromechanical system may be taken into consideration in some miniaturized/MEMS-based electrostatic transducers, such as in Ref. [23].

Analysis of multistability
At first, friction is not considered. The total potential energy of the system, U, consists of two parts: the strain energy of the cantilever beam, U b , and the magnetic potential energy, U m . The variation of the strain energy of the cantilever beam only depends on the deflection of the beam. The magnetic potential energy, however, not only is related to the distance between the magnets but also the alignment of the magnets within the field. Therefore, the horizontal and vertical distances between magnets, i.e., d x and d y , can notably affect the total potential energy, and thus the behaviour of the system. When the horizontal distance is fixed at d x ¼ 17 mm, the variation of the total potential energy with the vertical distance and the tip displacement is shown in Fig. 5a, and the corresponding projections on the Uw rel plane are given in Fig. 5b. It can be seen that the system is monostable when the two fixed magnets are close (i.e., d y is small). In this case, the two fixed magnets act as one stronger magnet. As the two fixed magnets are placed farther away from each other, the system exhibits monostability, bistability, and weak tristability.
The system has different stabilities at different combinations of d x and d y , and the possible scenarios are shown in Fig. 6 on the parametric plane of d x -d y which is divided into three areas according to their stabilities. Therefore, the multistability of the frictionfree (conservative) system can be determined analytically through the potential energy function, i.e., the sum of U b and U m in Eq. (13). Apparently, there is no potential energy function for a non-conservative force or system. Thus, the multistability of the system with the presence of friction can no longer be revealed by Eq. (13). However, it is interesting to see how those three areas will evolve on the parametric plane with the involvement of friction, and such a study of the effect of non-conservative forces on the multistability of a system is rare in the open literature.
Using numerical simulations, the system involving friction at each combination of d x and d y (a grid of 101 Â 101 is used) is integrated at different excitation frequencies with different initial conditions. The algorithm used for the categorization of different types of stability is shown in a flowchart in Fig. 14 in ''Appendix B''. Even though the multistability of the system is not dependent on the level of the external excitation, a proper selection of the excitation frequency and amplitude can improve the efficiency of the categorization. This is because in multistable systems, large excitations can easily result in interwell oscillations and small excitations will mostly lead to intrawell oscillations. In practice, this means that more tip displacement time history samples (or more sets of initial conditions) are needed to confidently determine the type of the multistability at a given combination of d x and d y .
In simulation, the excitation amplitudeÂ is fixed at 0.15 g, and two excitation frequencies, 3 Hz and 13 Hz, are used. Generally, it is easier for an oscillator to settle down into a well if the initial displacement is within the potential well (assuming the initial velocity is zero). Therefore, 15 samples of the initial tip displacement linearly spaced between -0.025 m and -0.001 m are used with each of the two excitation frequencies (initial velocities and transferred charge are zero). The categorized results with three different contact force ratios are shown in Fig. 7. It can be seen that the categorization scheme performs quite well, and the performance at boundaries can be enhanced by using more sets of initial conditions and/or excitation frequencies. In comparison with Fig. 6, the bistable and tristable areas have shrunk while the monostable area has enlarged. These differences between Figs. 6 and 7 seem to require only a small amount of friction to materialize. With the increase in the contact force ratio or the friction force, there is no obvious change of the areas of various stabilities.

Comparisons between various types of stability
The dynamics of the system can be highly dependent on the magnetic field. The magnetic field is mostly influenced by the relative locations between the magnets, i.e., d x and d y . As shown in the last section, the relative locations of magnets result in different types of stability. Besides, the system is non-smooth owing to the presence of friction. Thus, the steady state of the system can be dependent on the initial conditions [41], and various dynamic regimes can theoretically coexist, though only one of them can practically appear under certain conditions. Moreover, broad harvesting frequency span and large vibration amplitude are often beneficial to vibration energy harvesting. Therefore, discrete numerical frequency sweeps can be used to study the coexisting regimes (mainly the coexistence of intrawell and interwell oscillations) and compare between various types of stability.
To generate the discrete numerical frequency sweeps, 50 pairs of initial conditions w 0 ¼ w rel L; 0 ð Þ; _ w rel L; 0 ð Þ; 0 ð Þ are randomly selected from the set w rel L; 0 ð Þ; _ w rel L; 0 ð Þ ð Þ : w rel L; 0 ð Þ f 2 À0:025 0:025 ½ ; _ w rel L; 0 ð Þ 2 À1:5 1:5 ½ g where both w rel L; 0 ð Þ and _ w rel L; 0 ð Þ follow uniform distributions. The tip displacement w rel L; t ð Þ with these initial conditions under the discrete numerical frequency sweeps are shown in Fig. 8. The maximal and minimal vibration amplitudes in steady states at each excitation frequency are marked and connected through vertical line in this figure (blue circles for interwell oscillations and red asterisk for intrawell oscillations except the results of the monostable system shown in Fig. 8a), and the response of the latest instant at one frequency is used as the initial conditions for the integration of the system under the next frequency. Obviously, there are no coexisting behaviours in the monostable system. However, intrawell and interwell oscillations can coexist in the bistable system. The situation is similar in the tristable system. The bistable system has the broadest frequency span over which the oscillation is at high amplitude. Therefore, bistability is superior to both monostability and tristability in terms of broadband energy harvesting.
It is clear that there exist sudden jumps of the tip displacement in Fig. 8 in all the three kinds of systems. For a better explanation, the tip vibration amplitudes under the single forward (blue circles) and backward (red asterisks) frequency sweeps of the monostable system are shown in Fig. 9. It can be seen that there exist hysteresis and stiffness hardening, which are believed to be caused by the magnetic force.

The effect of the excitation level
In multistable systems, the excitation level often plays an important role in exciting high-energy orbit oscillations, especially when it comes to high and wide potential barriers between wells. The tip vibration amplitudes under discrete numerical frequency sweeps of the bistable system at different excitation levels are shown in Fig. 10. In combination with Fig. 9b, it can be observed that the frequency bandwidths at which interwell oscillations occur, shrinks considerably with the decrease in the excitation level.
For reduction of excitation level from 0.5 to 0.3 g, there still exists a considerable frequency span over which the system performs interwell oscillations, though intrawell oscillations completely become dominant at high frequencies. A further decrease to 0.1 g results in a significant reduction of the frequency span for interwell oscillations, and a much smaller excitation level ceases to excite any interwell oscillations, see Fig. 10c.

The effect of the friction
Triboelectrification is largely influenced by friction and can be enhanced through repeated rubbing [58]. Thus friction is essential in sliding-mode TEHs. On the other hand, friction is well known for dissipating energy by means of generating waves, atomic motions, and heat [59], and hence friction is also a limitation with regards to maximal energy harvesting. Therefore, the relationship between friction and energy harvesting (or electrical output) becomes even more complicated in TEHs and a study investigating this relationship is still missing. The mechanism by which friction affects triboelectrification is beyond the scope of this paper. Nevertheless, it is intended to study the effect of friction on the performance of the current TEH by investigating its influence on the system dynamics.
There are various combinations of triboelectric materials and each combination has its own values of coefficients of friction (both static and kinetic). Besides, the contact between the slider and the patches  Fig. 10a, it can be seen that the distribution of the intrawell oscillation spreads from high frequencies down to relatively low frequencies with the increase in the contact force ratio. In addition, the jump frequency (due to the hardening effect) shifts towards lower frequencies. In comparison, low frequencies are unlikely to result in intrawell oscillations.
Since the base excitation amplitude is fixed in terms of the applied acceleration, the resultant base displacement is inversely proportional to the forcing frequency square. Thus, the displacement response tends to be small at high frequencies.

Basins of attraction
Since the dynamics of the system can be sensitively dependent on the initial conditions, basins of attraction are computed to investigate the eventual dynamic  The steady-state of the system at each sampling point is categorized according to the attractor onto which the system eventually settles down. The attractors are categorized into period-1 intrawell oscillations (shown in red), non-periodic or high-periodicity intrawell oscillations (in black), period-1 interwell oscillations (in blue), and non-periodic or highperiodicity interwell oscillations (in green).
The generated basins of attraction of the bistable system under the excitation ofÂ ¼ 0:3 g and A ¼ 0:5 g at two different excitation frequencies are shown in Fig. 12. It can be seen that the increase in the excitation level has profound effects on the evolution of the basins of attraction. At f ¼ 8Hz, the spiral basins both for the period-1 interwell oscillations and the non-periodic or high-periodicity interwell oscillations become scattered with the increase in the excitation level, and the basin of the former becomes dominant, see Fig. 12a. At f ¼ 13Hz, the weak spiral structure atÂ ¼ 0:3 g disappears due to the rise of the excitation level toÂ ¼ 0:5 g, and the basins of the nonperiodic or high-periodicity intrawell oscillations and the period-1 interwell oscillations seem to intermingle in the off-centre area atÂ ¼ 0:5 g. Additionally, it has been found that the system undergoes purely interwell oscillations at lower frequencies (such as at f ¼ 2 Hz) even for smaller excitation levels. The corresponding basins of attraction are purely green and are not presented here.
Similarly, the basins of attraction of the tristable system under the excitation ofÂ ¼ 0:3 g andÂ ¼ 0:5 g at these two excitation frequencies are obtained and shown in Fig. 13. Interesting structures appear at f ¼ 8 Hz, see Fig. 13a. Basins of nonperiodic or high-periodicity intrawell oscillations arise around the two side stable equilibria while the basin of period-1 intrawell oscillations emerges in the middle and has a fractal boundary. In addition, the tail of the spiral basin of the period-1 interwell oscillation is separated by an intermingled spiral basin of the period-1 and non-periodic or high-periodicity intrawell oscillations, and the increase in the excitation level thickens the tail of the basin of the period-1 interwell oscillation. At the relatively higher frequency of f ¼ 13 Hz, the system merely performs intrawell oscillations at a relatively low excitation level, and basins of interwell oscillations start to emerge in a fractal structure with the increase in the excitation level; see Fig. 13b. The tristable system only performs intrawell oscillations at low frequencies (such as at f ¼ 2 Hz) even when the excitation level is increased (the corresponding basins of attraction are not given here), and relative to the bistable system, the tristable system always settles down onto attractors of intrawell oscillations at low frequencies, which is not beneficial for energy harvesting.

Conclusions
This paper presents the investigation of a new TEH which works in a sliding mode and involves three kinds of stabilities. The design concept is first introduced. The modelling of the mechanical vibration of a cantilever beam with a tip mass subjected to both magnetic and frictional forces under base excitation is then presented. The modelling of the electrical dynamics of the sliding-mode TEH based on the established mechanical system is presented accordingly. A numerical scheme, which is based on the TR-BDF2 method, for solving the system comprised of the non-smooth mechanical and the stiff electrical systems is discussed. Detailed numerical simulations are then carried out. This paper includes the first studies of both the sliding-mode TEH and magnetic multistability in TEHs from the perspectives of structural dynamics. The main conclusions are drawn as follows: 1. The comparison between the coupled and uncoupled electro-mechanical models (with and without the consideration of the electrostatic force between the electrodes, respectively) suggests that the electrostatic force can be ignored, which makes the uncoupled model preferable in the dynamical analysis. 2. The effect of the non-conservative force (the friction force) on the multistability of the system is investigated. It is found that the distribution of the multistability on the d x -d y (two distances describing the relative locations of two magnets fixed to the base in relation to one magnet embedded in the slider) parametric plane changes when friction is involved, and the proportions of occurrence of the bistability and tristability decrease while that of the monostability increases. 3. Among the three types of stability, bistability is found to provide the broadest frequency band over which the system oscillates on high-energy orbits or in interwell oscillations. Additionally, the stiffness-hardening effect caused by the magnetic force appears in systems of all three kinds of stabilities. 4. The base excitation level plays an important role in bistable or multistable system for broadband energy harvesting. Low excitation levels will not enable the systems to cross the energy barriers between the potential wells and achieve interwell motions. In addition, the increase in friction will transfer intrawell oscillations from high frequencies down to low frequencies and shrink the frequency band over which energy harvesting is effective. 5. The study of the basins of attraction reveals that at low frequencies the bistable system can undergo purely interwell oscillations while the tristable system will undergo intrawell oscillations. Intermingled basins can appear in both bistable and tristable systems. An increase in the excitation level can break the basins into discrete pieces and/ or points. ðA:4Þ

Appendix B
The algorithm used for the categorization of different types of stability is shown in a flowchart in Fig. 14, where f is the excitation frequency array, w rel L; 0 ð Þthe initial tip displacement array, t s the time within steady state, and w mean h ð Þ the mean of w max2 h ð Þ and w min2 h ð Þ.