Global Dynamics and Optimal Design of a New Highly Efficient Nonlinear Energy Sink Based on Magnetic-Elastic Impacts with Negative Stiffness


 A new highly efficient elastic-impact bistable nonlinear energy sink (EI-BNES) based on magnetic-elastic impacts with negative stiffness and bistability is proposed and optimized through global dynamical analysis. The EI-BNES has better robustness and higher energy dissipation rates with nearly more than 96.5\% for broadband impulsive excitations than the traditional cubic NESs and single bistable NESs. The structure of negative stiffness impacts is realized by reasonable layout of permanent ring magnets and springs. A two-degree-of-freedom (two-DOF) elastic-impact system is established to describe the coupled nonlinear interaction between the main structure and the attached EI-BNES. A global Melnikov reduction analysis (GMRA) is proposed to study global dynamics and homoclinic bifurcations of the reduced two-dimensional subsystem, which is used to explain the mechanism of nonlinear targeted energy transfer (TET) and detect the threshold of impulsive amplitudes of EI-BNES for in-well and compound motions between in-well and cross-well resonance responses. A special type of saddle-center equilibrium points is also found in the non-smooth system of the EI-BNES and can be used to effectively increase the energy dissipation rates. The optimal design criterion of the tuned EI-BNES for better dissipation performance is also first discussed based on the GMRA and numerical techniques for calculating the Melnikov function of the non-smooth systems. The effectiveness of the analytical GMRA is also verified by numerical simulations.


Introduction
Structural vibration problems are common in mechanics, aerospace, civil engineering and other fields. They generally have negative effects on mechanical operations and structural responses. For example, chatter instability which may occur in the process of tool turning often results in the decline of machining quality [1], the engines fixed on the wings of an aircraft are prone to produce unsatisfactory vibrations and noise during operation, which leads to the reduction of flight comfort [2], and excessive vibrations of building structures threaten people's life and property safety under the actions of wind, earthquake or other dynamic loads [3]. Hence, how to absorb vibrations efficiently and dissipate vibration energy quickly to reduce the disadvantageous effects of vibrations is a hot issue in the field of vibration control.
In order to solve the problems of structural vibrations effectively, Frahm [4] has first proposed the concept of dynamic vibration absorbers, it is also known as a tuned vibration absorber (TVA). The dynamic vibration absorber connects a small mass to a main structure by a linear spring with appropriate stiffness. Furthermore, the other relevant linear dissipation devices have also been studied in [5][6][7][8]. The linear vibration absorbers are demonstrated to be effective for vibration reduction only at a priori natural frequency of the main structure. Hence, how to employ nonlinearity in dissipation devices to expand vibration suppression bandwidth of frequencies and improve robustness has always attracted the attention of researchers.
In the mid-20th century, Roberson [9] has connected a small mass by springs with linear and cubic nonlinear stiffness to a single-degree-of-freedom (SDOF) main structure under sinusoidal excitations and found the nonlinear vibration absorber has a wider vibration reduction frequency bandwidth than the TVA. Furthermore, Gendelman et al. [10] and Vakakis et al. [11] have proposed a class of passive dissipation structures with cubic nonlinear stiffness and linear damping and found the impulsive excitation energy of the main structure can be transferred to the additional light weight structure in one direction and locally dissipated by the damping. This phenomenon is called TET and these nonlinear attachments are named as NESs. Since then, more and more scholars have carried out a lot of studies on designs, dynamical analysis, experiments and applications of NESs, and ones can refer to the monograph [12] and the recent comprehensive reviews [13][14] for the latest research progress.
It is well known that the cubic NESs have low energy dissipation rates under ultra-low and ultra-high input impulsive excitations [15][16]. In order to further improve the vibration reduction performances of NESs and broaden impulse magnitude, many scholars had begun to discuss NESs with the characteristic of bistable since 2014. By compressing linear springs attached to the main structure to realize the negative linear and nonlinear stiffness, AL-Shudeifat [17] has proposed a bistable NES and showed that the bistable NES has more efficient vibration absorption performances than the cubic NESs under impulsive loads. In addition, the detailed dynamics of bistable NESs with negative linear stiffness and cubic stiffness under impulse excitations have been respectively studied by an analytical method in [18] and numerical techniques employed to study the strongly modulated and chaotic cross-well oscillations in [19]. On this basis, Romeo et al. [20] have furthermore studied the transient and chaotic low-energy transfers in a bistable NES coupled to a linear oscillator. Habib and Romeo [21] have established a tuned bistable NES (TBNES) connected to one of two coupled symmetric linear oscillators and shown that the absorption capability of the TBNES is higher than the TMD and the purely cubic NES for a wide range of impulsive loads. Qiu et al. [22] have employed the complexification-averaging method to study the TET and optimal design of a bistable NES with negative linear stiffness and cubic stiffness under periodic excitations. Yang et al. [23] have studied the enhanced TET for adaptive vibration suppression of fluid conveying pipeline. In additional, Yao et al. [24][25] have respectively studied the NESs for rotor systems with bi-stable and multi-stable characteristics.
Many scholars have also considered the advantage of fast energy consumption under impacts and obtained many studies about NESs with impact characteristic. In the study of structural protections in earthquake, a vibro-impact nonlinear energy sink (VI-NES) has been proposed in [26] and verified that it can transfer the seismic energy of the main structure passively to the additional nonlinear attachment. In addition, Lee et al. [27] have studied the periodic orbits and frequency energy orbits of elastic and inelastic VI-NESs and shown that the VI-NESs have the advantage of being able to absorb and dissipate most of the energy in the main structure in a short enough time scale. The dynamics of VI-NESs have also been studied analytically by means of the multi-scale method and power series expansion in [28]. AL-Shudiefat et al. [29] have proposed an asymmetric non-smooth NES and verified the highefficiency performances of unilateral VI-NES under impulse excitations. Furthermore, Gendelman et al. [30] have studied the dynamics and energy transfer of a periodically forced VI-NES and shown that the strongly modulated response contains the random distribution of resonant periodic and non-resonant motions. Gourc et al. [31] have built a new type of VI-NES and obtained that the periodic forced VI-NES is effective for TET. A two-DOF parallel VI-NES under periodic and transient excitations has been studied in [32] to show that the vibration control is more effective when the VI-NES is activated with two impacts per-cycle. A new VI-NES combined with cubic stiffness and bilateral barriers has been analyzed by harmonic balance method in [33]. Qiu et al. [34] have presented a VI-NES optimal criteria for vibration control under periodic and transient excitations and verified that the parallel VI-NES designed by this standard can achieve the best energy transfer in a certain frequency range.
Considering the current developments of NESs, it is possible to combine the double advantages of bistability and impact to design a new highly efficient NES. Furthermore, how to develop analytical methods for optimal design of NESs and realization of highly efficient TET has also important value of engineering applications. So, a new more efficient EI-BNES combining the characteristics of impacts and bistability is first considered to design by magneticelastic negative impacts in this paper. Based on the geometric explanation of the Hamiltonian dynamics in terms of the slow manifolds of the dynamics, a new analytical GMRA to enhance TET in the damped systems will be further presented.
The highlights are listed as follows. Firstly, the EI-BNES has better robustness and higher energy dissipation rates than the traditional cubic NESs and single bistable NESs. Furthermore, the structure of negative stiffness impacts is realized by reasonable layout of permanent ring magnets and springs and the main process and results are shown in detail in the Theorem 1. In addition, the GMRA is proposed to study global dynamics and homoclinic bifurcations of the reduced two-dimensional subsystem. The special saddle-center equilibrium points which can be used to effectively increase the energy dissipation rates are also found in the non-smooth impacts system of the EI-BNES. Finally, Melnikov function is calculated by numerical method and first used to discuss the optimization of the EI-BNES system. The effectiveness of the analytical GMRA is also verified by numerical simulations.
The structures of this paper are presented as follows. In the second section, the new highly efficient EI-BNES is proposed. Furthermore, the energy dissipation rates and dynamic analysis of the EI-BNES and the realization of the negative stiffness impacts are discussed. In the third section, the GMRA is proposed to study the mechanism of TET and homoclinic bifurcations of the reduced two-dimensional subsystem. In the final section, the Melnikov function is solved by the combination of implicit functions and numerical methods and used to detect the chaotic thresholds of the homoclinic bifurcations of the system. Furthermore, the Melnikov function and the chaotic threshold curve are used to optimize the EI-BNES model and the effectiveness of the GMRA for optimal designs is verified by numerical simulations.

Model of the EI-BNES
In this section, the new highly efficient EI-BNES is inspired by a bilateral elastic impacts oscillator studied in [35]. The model under investigation consists of a linear oscillator (LO) with mass m 1 , a c 1 -viscous damper, a k p -spring and a mass m 2 attached to the mass m 1 as the EI-BNES. The mass m 2 is vertically connected with a pair of symmetric linear k 1 -springs and slides on a thin rod which is fixed on the main structure m 1 and has the viscous damping with coefficient c 2 . When the mass m 2 is in the center of the rod, the k 1 -springs are compressed and its length is l at this time. Furthermore, when the elastic impacts between the m 2 and the k 1 -springs happen instantaneously, the k 1 -springs are in the original state without compressed or stretched, and their original length of the k 1 -springs is denoted as L. The two ends of the thin bar are respectively wound with k 2 -linear symmetric springs. The model is shown in Figure  1. x for m 2 2 x for m When an impulsive excitation X is acted on the main oscillator m 1 , the mass m 2 will move from the left stable state. If the moving relative distance of the main oscillator m 1 and the mass m 2 reaches a certain value a or −a, the elastic impact happens between the mass m 2 and the k 2 -spring at this time. When the mass m 2 moves to the center of the rod, the mass m 2 is in an unstable state because of the k 1 -springs are in compression state. In additional, the bistability will also appear because of the influence of the springs with stiffness coefficients k 1 and k 2 . In this process, it is assumed that the mass m 2 will not impact with the leftmost or rightmost ends of the thin bar, namely, there will be no rigid impacts. When the mass m 2 slides on the thin rod, the viscous damping with the coefficient c 2 is also applied here.
The coordinate origin is selected at the center of the EI-BNES frame shown in the Fig. 1. According to the Newton's second laws and the geometric configuration of the k 1 -springs, the dynamical equations of the mass m 1 and the lightweight m 2 with a piecewise-smooth restoring force are described as follows: (1) where x 1 and x 2 are the displacements of the main oscillator m 1 and the mass m 2 , a = √ L 2 − l 2 , the dots denote derivations with respect to the time t.

Energy dissipation rates of the EI-BNES
For impulsive excitations applied to the LO with different magnitude X, namely, the initial velocity of the LO iṡ x 1 (0) = X, the energy dissipation rates E NES (t) of the EI-BNES can be analytically obtained on the time internal [0, t] in the following form It is well known that achieving higher dissipation rates in a short time under broadband impulsive excitations is an important index to evaluate whether NESs have better vibration reduction effects. Therefore, based on energy absorption efficiency, energy absorption speed and range of impulsive excitations, the vibration reduction effect of the EI-BNES proposed in this paper is evaluated and the parameters k p = 1, c 1 = 0.001, c 2 = 0.01, k 1 = 0.2, L = 1, l = 0.9, m 1 = 1, m 2 = 0.05 are chosen in the following numerical simulations.
The energy dissipation rates E NES (t) of the EI-BNES under different impulsive inputs X applied to the LO for t = 100s and t = 200s are first studied. The energy dissipation rates E NES (t) of the EI-BNES will be discussed when the stiffness coefficient k 2 is respectively selected as positive, zero and negative. Where k 2 = −0.05, k 2 = 0 and k 2 = 0.05 are respectively chosen to study the energy dissipation rates E NES (t) within 100s and 200s under different energy impulsive excitations X, the results of energy dissipation are shown in Figure 2(a) and Figure 2(b), respectively.
When k 2 = −0.05, the energy dissipation rates of the EI-BNES are significantly higher than the cases k 2 = 0 and k 2 = 0.05 under different energy impulsive excitations X in the Fig. 2. When the different impulsive excitations X of the LO are in low level inputs, namely, X = 0 ∼ 0.4, the negative k 2 can maintain a higher dissipation rates than the positive stiffness and zero stiffness. The maximum difference of the dissipation rates among negative stiffness, positive stiffness and zero stiffness is even as high as about 40%. When the impulsive excitations X of the LO are in medium inputs, namely, X = 0.4 ∼ 0.7, the dissipation rates curves are almost the same under the three different stiffness conditions. However, the dissipation rates curve under negative k 2 is still slightly better than the other two cases. When the different impulsive excitations X of the LO are in high energy level, namely, X = 0.7 ∼ 1, the dissipation rates curves under the three conditions all begin to decrease, but the reduced amplitudes of dissipation rates curves in the cases of positive and zero stiffness are obviously higher than the case of negative stiffness. For the range of low impulsive inputs, we can also obviously observe the high amplitudes and frequency variation of the dissipation rates curves of the EI-BNES under the positive and zero stiffness, but the dissipation rates curve under the negative stiffness has slight change which shows the great robustness of the EI-BNES for different level impulsive inputs. Comparing the Fig. 2(a) with the Fig. 2(b), when the stiffness coefficients k 2 are chosen as positive or zero and the impulsive excitations X = 0 ∼ 0.4, the energy dissipation rates E NES (t) for t = 200s are significantly higher than the dissipation rates E NES (t) for t = 100s, the maximal amplitude difference of the dissipation rates can even reach about 15% within the two time interval. It can be seen that when the time t is slightly larger than 100s, the value of E NES (t) is still getting bigger, that is to say, there is still dissipation phenomena in the EI-BNES when the coefficient k 2 is positive or zero. However, the dissipation rates of the EI-BNES only slightly change between t = 100s and t = 200s when the k 2 is negative which means both the EI-BNES and the LO have almost reached the stable states at t = 100s and the presented EI-BNES with negative stiffness can even complete the vibration reduction in a shorter time. In order to observe clearly the level of dissipation rates in the presented EI-BNES with negative stiffness, the energy dissipation rates for k 2 = −0.05 within 100s under different impulsive excitations X are displayed in Figure 3, which shows initial input energy at more than 96.5% is dissipated by the damped EI-BNES. How to realize the negative stiffness impacts from physical models will be discussed in the next subsection, in fact, it is easy to explain the reason of the higher energy dissipations within shorter times in physical sense. When the mass m 2 and the k 2 -spring impact instantaneously, the k 2 -spring with negative stiffness will give the mass m 2 a force in the same direction as its own motion and cause a larger velocity difference between the LO and the EI-BNES. Therefore, according to the equation (3), we can easily get that there indeed has a higher dissipation rate in the EI-BNES under this case.
Next, the changes of energy dissipation rates E NES (t) with time t will be analysed by respectively taking X = 0.2, X = 0.5, X = 0.8 as the representative in low, medium and high impulsive regions. The stiffness coefficient k 2 = −0.05 and the other parameters are the same as above. The changes of dissipation rates E NES (t) with time t are presented in Figure 4. For X = 0.2, X = 0.5 and X = 0.8, the energy dissipation rates of the EI-BNES respectively tend to be stable at about t = 90s, t = 50s and t = 60s in the Fig. 4. The energy dissipation rates E NES (t) of the three conditions are all nearly 100%. That is to say, most of the energy of the model has been dissipated by the damped EI-BNES, which directly shows one way TET from the LO to the EI-BNES. Moreover, when the energy impulsive excitation X is in the medium level, the energy dissipation process is faster to be completed in less than one minute.

Time(s)
Next, what k 2 -value can make the EI-BNES has better energy dissipation effect on the whole under the impulsive internal X = 0 ∼ 1 will be discussed. The variation range k 2 = −0.08 ∼ −0.01 will be set the other parameters are the same as above. The relationships among different impulsive excitations X, energy dissipation rates E NES (t) and different stiffness coefficients k 2 for t = 100s are shown in Figure 5(a) and Figure 5(b). It is found that the energy dissipation effect of the EI-BNES is better at k 2 = −0.05 which is shown in the red line for different impulsive excitations X applied to the LO in the Fig. 5(a).
Therefore, how to develop analytical methods to reveal the TET mechanism and optimize the tuned EI-BNES to achieve the highest dissipation rates and the fastest dissipation process is a challenging topic. These are directly related to the presented GMRA and model optimization which will be discussed in detail in part 4.2.

Design of the model with negative stiffness impacts.
According to the EI-BNES in section 2.1 and the analysis of energy dissipation rates in section 2.2, it is found the dissipation rates of the EI-BNES can be all up to 96.5% and the EI-BNES completes the vibration reduction process within 1 minute when the coefficient k 2 = −0.05N/m and the impulsive inputs X = 0 ∼ 1m/s. How to realize the negative stiffness impacts in the real physical model is a new problem. In this paper, the realization of the negative stiffness impacts is discussed in detail and it is inspired by the design of a magnetic vibration absorber with tunable stiffness presented in [36]. The specific realization model is shown in Figure 6. The mass m 2 is vertically connected with a pair of symmetric linear k 1 -springs and it slides on the thin rod. A vicious c 2 -damping is also considered in the process of sliding. At the left (or right) side of the mass m 2 , a pair of symmetrical linear k 3 -springs are connected with a wood block ignored size and mass. The wood block passes through the thin rod and it can also slide on the thin rod. When the mass m 2 is in the center of the rod, the k 1 -springs are compressed to l-lengths. Furthermore, when the elastic impacts between the m 2 and the wood block happen instantaneously, the k 1 -springs and the k 3 -springs are both in their original state without compressed or stretched, and their original lengths of the k 1 -springs and the k 3 -springs are respectively denoted as L and l. The four corners of the frame are embedded three permanent magnets with paramagnetic poles. At the leftmost and rightmost sides of the rod are both inlaid a permanent magnet. In addition, two sides of the mass m 2 are respectively inlaid with a permanent magnet. The magnetic poles of the permanent ring magnets at all positions are shown in the Fig. 6 (black and white show different magnetic poles).
The width of the permanent magnet is L * and the horizontal distance from the center of the three permanent magnets to the center of the mass m 2 is R. The distance between the center of the permanent magnet on the leftmost (rightmost) side of the bar and the center of the three permanent magnets embedded in the upper or lower corner of the leftmost (rightmost) end of the EI-BNES frame is d. The height of the permanent magnet is D 1 , the width of the rod is D 2 and the distance between the moving mass m 2 and the origin of the mass m 2 is x, where, r = R + L * . The mass m 2 will move from the left stable state under a certain impulsive excitation. If the moving distance x of the mass m 2 reaches a certain constant a or −a, the elastic impact happens between the mass m 2 and the wood block at this time. When the mass m 2 moves to the center of the rod, the mass m 2 is in an unstable state because of the k 1 -springs are in compression state. In additional, the bistability will also appear because of the influence of the springs with stiffness coefficients k 1 and k 3 and the magnetic force. When the mass m 2 compresses the wood block to the left (right), the permanent magnet embedded in the leftmost side (rightmost side) of the frame of the EI-BNES will act on the permanent magnet embedded on the left (right) side of the mass m 2 . In this process, it is assumed that the wood block will not impact with the permanent magnets at the leftmost and rightmost sides of the thin rod.
The coordinate origin is selected at the center of the frame of this model shown in the Fig. 6. According to the Newton's second laws and geometric analysis, the restoring force induced by the permanent magnets and the k 3 -springs is obtained as follows: where K 1 and K 3 are the stiffness coefficients of the permanent magnets forces presented in [36].

The derivation of the negative stiffness coefficients
Theorem 1. The aforementioned restoring force f (x) can be obviously written as if the stiffness coefficient of the k 3 -spring is taken as k 3 = −K 3 l 2 for K 3 < 0.

Proof:
The case x ≥ a will be only discussed such that the equation (4) becomes The Taylor expansion of the third items of the equation (6) is carried out at x = a and the following equation can be obtained If K 3 < 0 and the following equation is considered Then, f (x) can be approximately rewritten as follows: The case for x ≤ a can be discussed similarly, and the final restoring force is given by Remark 1. The Theorem 1 has shown that if the stiffness coefficients K 1 and K 3 are both negative, that is to say, K 1 < 0 and K 3 < 0, then the negative stiffness impacts can be realized.
Next, whether the negative stiffness coefficients K 1 and K 3 can be obtained by adjusting the parameters of the model will be discussed.  According to the above analysis, it is easy to find that the repulsive force on the left side permanent magnet of the mass m 2 can be adjusted by the length r. Similarly, the attractive force on it can be adjusted by the length R and d. Because of r = R + L * and L * is a fixed value, so the modifications of r and d are only needed to consider to adjust the coefficients K 1 and K 3 such that K 1 < 0 and K 3 < 0. Next, the numerical simulations of looking for the magnetic force coefficients K 1 and K 3 under different parameters r and d will be considered. Among them, the parameters L * , D 1 , D 2 , m r , m c , u 0 and m * are shown in Table 1. u 0 is the void permeability, m r is the magnetic moment of repulsive magnet, m c is the magnetic moment of attraction magnet, m * is the magnetic moment of the permanent magnet inlaid on the mass m 2 . The variable conditions of the magnetic coefficients K 1 and K 3 with r and d are shown in Figure 7(a) and Figure 7(b), respectively.  Fig. 7. The coefficients K 1 and K 3 can both be two conditions which are positive, negative. Furthermore, the parameters K 1 and K 3 can both be taken relatively large negative values. That is to say, k 3 as a relatively small positive number can be obtained because of the relation k 3 = −K 3 L 2 1 . In this case, the value k 3 is also in line with other stiffness coefficients of springs. Even more fortunate is that K 1 = −0.05N/m can be obtained which is agreement with the discussion of negative stiffness impacts in section 2.2. The concrete parameters K 1 and K 3 under the two conditions are shown in Table 2.
According to the above analysis, the cubic term generated by magnetic force and the cubic term produced by the k 3 -springs can be eliminated, so only the linear term generated by magnetic force is left. When r = 0.171m, d = 0.049m, K 1 = −0.05N/m, combining with the discussion of motion equations in section 2.3, the model can realize negative stiffness impacts. Therefore, the feasibility of the negative stiffness impacts of the EI-BNES is fully verified.

Dynamic analysis of the EI-BNES
Next, the time history of the LO and the EI-BNES, the phase portraits of the EI-BNES, and the corresponding frequency spectrums will be presented with the changes of the input impulsive loads to show the internal resonant capture of the nonlinear TET from the LO to the EI-BNES. The parameters L = 1, l = 0.9, m 1 = 1, m 2 = 0.05, k 2 = −0.05, k p = 1, c 1 = 0.001, c 2 = 0.01, k 1 = 0.2 are selected. The impulsive excitations X are chosen as X = 0.05, X = 0.2, X = 0.5, X = 0.8, respectively. Numerical simulations for 200s with the aforementioned four initial impulsive conditions are respectively shown from Figure 8 to Figure 10. In the Fig. 8(a), the LO and the EI-BNES almost stop vibration at about 60s and 100s, respectively, and they both tend to their equilibrium positions at the end. It is also found that the EI-BNES produces alternating in-well and cross-well oscillations in 0 ∼ 40s and vibrates around one of the two equilibrium positions in the following stage in the Fig. 8(a) and the Fig. 8(b). Moreover, the EI-BNES can generally complete the energy dissipate process for no more than 60s which sufficiently implies that the process of vibration reduction of the EI-BNES is fast. According to the Fig. 8(c) and the Fig. 8(d), the wavelet spectrums show that the 1: 1 resonance and the low-frequency resonance of the EI-BNES are excited which increases the dissipation rates of the EI-BNES. In the Fig. 9(a), the LO and the EI-BNES almost stop vibration at about 40s and 60s, respectively, and they both tend to their respective stable positions. In addition, it is seen that the EI-BNES generates large amplitude cross-well oscillation within about 0 ∼ 30s in the Fig. 9(a) and Fig. 9(b). According to the corresponding spectrums Fig. 9(c) and the Fig. 9(d), the time-frequency analysis indicates the LO and EI-BNES have the same resonance frequency (1:1 resonance capture occurs between the LO and EI-BNES). In the Fig. 10(a) and the Fig. 10(b), the LO and the EI-BNES almost stop vibration at about 60s and 90s, respectively, and they also both tend to their respective stable positions. In additional, it is seen that the EI-BNES generates large amplitude cross-well oscillation within about 0 ∼ 60s. The corresponding spectrums in the Fig. 10(c) and the Fig. 10(d) also reflect that the LO and the EI-BNES have the same resonance frequency (1:1 resonance capture occurs between the LO and EI-BNES).

Dimension reduction of the motion equations
The Melnikov method for non-smooth systems is an important analytical method [38][39][40][41][42][43][44] used to analyze the global chaos of piecewise-smooth systems. Hence, based on this mechanism, it is necessary to reduce the dimensions of the EI-BNES from four to two and analyse the sufficient conditions to induce transient chaos by means of the Melnikov function of non-smooth systems.
Firstly, let Furthermore, the following equation is considered Substituting the equations (11) and (12) into the equations (1) and (2), the following equations can be obtained |µ| ≤ a, |µ| > a, where the derivatives with respect to τ. Next, the terms with ε which are regarded the perturbation parts of the equations (13) and (14) are discarded, then, the following equations can be got |µ| ≤ a,ẍ |µ| > a,ẍ The analytic solution ofẍ 1 + x 1 = 0 is easy to be obtained as follows: Substituting the equation (17) into the second equation of the equations (15) and (16), the following results can be obtained |µ| ≤ a,μ |µ| > a,μ Let then, the equations (18) and (19) can become |µ| ≤ a,μ |µ| > a,μ Furthermore, because of then, the equations (21) and (22) become |µ| ≤ a, The corresponding Hamiltonian functions of the equations (24) and (25) can be obtained as follows: If we consider H(µ, ν) = 0, the following equations can be got Furthermore, if H(µ 0 , 0) = 0 is considered, the following equation can be obtained where µ 0 > 0 can be obtained by finding the root of the equation (28) and (µ 0 , 0) is the intersection point of the right half homoclinic orbit with the µ axis. Next, the potential energy function and homoclinic orbits are numerically simulated to observe the dynamics of the reduced system. In Figure. 11(a), the potential energy has the local maximum at u = 0, a global homoclinic bifurcation will be produced, which implies a transition from in-well oscillation to cross-well oscillation. Further, the potential energy about k 2 increases with the decrease of k 2 . In Figure.

Numerical solutions of the Melnikov function
The four dimensional system (1) and (2) have been respectively reduced to the two dimensional system (24) and (25). According to the analysis of Melnikov function [42], the following result can be got Considering the parity of the trigonometric functions, the equation (29) where The detailed derivation process of the equation (30) is shown in Appendix A.
According to the Melnikov's theorem, the following equations can be obtained Hence, the following equation can be obtained from the equation (32) Because of the equation (33) can become

Optimization of the EI-BNES dynamical model
The EI-BNES is mainly considered the characteristic of high dissipation rates. It is noticed that the dissipation rates E NES (t) relate to whether the EI-BNES happens transient chaotic cross-well motion occured near homoclinic orbit when the impulsive excitations X are in low impulsive excitations. Furthermore, it is found that there is a certain relationship about chaotic thresholds between I and A in the equation (35). It is easy to find that A is the impulsive excitations X of the LO and is given by the equation (17). Once the parameters of the EI-BNES are selected, I will be fixed according to the equation (35). At this time, if A is larger than I, then the chaos motions will occur. That is to say, the EI-BNES will have chaotic motions under any impulsive excitations X > I. It has been found that I 1 and I 2 mainly relate to k 1 , k 2 , and U mainly relates to k 1 in the equation (35). This means that I is mainly determined by k 1 and k 2 , hence, the relationship between k 1 and k 2 is first needed to find to better optimize the EI-BNES in order to improve its effect of vibration absorption.

Analysis of the relationship between k 1 and k 2
Next, we will consider the EI-BNES to be a single nonlinear oscillator to detect the relationship between k 1 and k 2 for bifurcations. The non-smooth conservative system describing the motion of the EI-BNES can be obtained as follows: where the piecewise-smooth restoring force f (x) is For the convenience of discussion, the values k 2 = −0.05, L = 1, l = 0.9 are fixed in this section. The restoring force f (x) for different k 1 is illustrated in Figure 12 in order to detect the bifurcation value k 1 .
It can be seen from the Fig. 12(a) and the Fig. 12(b) that k 1 = 0.025 (the yellow line) and k 1 ≈ 0.13 (the blue line) respectively are the bifurcation values such that the numbers of the zero points for the restoring force f (x) change from three to five and from five to three with the changes of k 1 . It is noted that the zero points x = 0, x = a and x = −a always exist for any k 1 and there exist two thresholds of k 1 such that the other two zero points will emerge or disappear, which will be detected for the accurate boundary of k 1 by the following analytical technique.
By observing the behavior of the critical yellow curve shown in the the Fig. 12(a), the critical vale of k 1 for the numbers of zero points changing from three to five can be obtained by studying the limit of the restoring forcing f (x) as t → ±∞. When the case x > a is considered (x < a can be discussed in a similar way), the restoring forcing f (x) can be written as and its limit lim x→+∞ f (x) =Ã exists if and only if k 1 = −0.5k 2 , which is the bifurcation value such that the numbers of zero points change from three to five. The bifurcation value is the same as the numerical simulation in Fig. 12(a) for k 1 = 0.025 assigned.
By observing the behavior of the critical blue curve shown in the the Fig. 12(b), the critical vale of k 1 for the numbers of zero points changing from five to three can be obtained by letting the right derivative f ′ + (a) = 0 or the left derivative f ′ − (−a) = 0 as follows: which can get the bifurcation value such that the numbers of zero points of the restoring forcing f (x) change from five to three. It is also noted that the critical value k * 1 ≈ 0.13 is also the same as the numerical simulation in Fig. 12(b). Next, the changes of stability of these equilibrium points and the global dynamics of the equation (36) will be analysed and depicted in Figure 13. Before further analysis, we give a new definition of saddle-center equilibrium points and the corresponding unique bifurcations of non-smooth systems. In the Fig. 13(a), the red solid lines represent the center equilibrium points, the black dotted lines represent the saddle-type equilibrium points, and the red and black intersecting lines represent the saddle-center equilibrium points. Furthermore, the EI-BNES has bifurcations not only about the number of equilibrium points, but also about the stability and the types of the equilibrium points. There are three equilibrium points when 0 ≤ k 1 ≤ −0.5k 2 (A area), five equilibrium points when −0.5k 2 ≤ k 1 ≤ k * 1 (B area) and three equilibrium points when k 1 > k * 1 (C area). It is easy to observe that saddle-center splitting bifurcation occurs at k 1 = −0.5k 2 and saddle-center and central convergence bifurcation occurs at k 1 = k * 1 . Combining the actual physical model, the bifurcations of these equilibrium points and their corresponding dynamics will be analysed in detail through phase portraits by selecting some specific values in the three areas, which are respectively displayed from Figure 13(b) to Figure 13(f). The concrete assignments for k 1 are k 1 = 0.02 in the area A, k 1 = 0.05 and k 1 = 0.08 in the area B, k 1 = 0.13 and k 1 = 0.16 in the area C. It is easy to find that (0, 0) always exists as the saddle-type equilibrium point and (±a, 0) are the saddle-center type equilibrium points for k 1 < k * 1 . Two other center equilibrium points are split at the bifurcation value k 1 = −0.5k 2 . They are close to the saddle-center equilibrium points (±a, 0) with the increase of the parameter k 1 and finally collide the (±a, 0) and merge into center equilibrium points at (±a, 0) for k 1 = k * 1 . The Fig. 13(b) displays a representative phase portrait for a given k 1 = 0.02 satisfying k 1 < −0.5k 2 . There exist one saddle-type and two saddle-center equilibria such that the EI-BNES shows strong instability and sensitive dependence on initial conditions. The main reason is the negative stiffness characteristics of the structure due to the double role of impacts with negative stiffness and bistability. Once the movement displacement of the mass m 2 satisfies |x| > a, the speed and the displacement of the mass m 2 will diverge to infinite so that the mass m 2 has the risk of hitting the leftmost or rightmost fixed boundary of the frame. However, this case should be avoided in our design of EI-BNES. Hence, we must assure k 1 > −0.5k 2 in the following discuss.
The Fig. 13(c) and Fig. 13(d) respectively display two representative complicated phase portraits for two given k 1 = 0.05 and k 1 = 0.08 satisfying −0.5k 2 < k 1 < k * 1 . There exist one saddle at (0, 0), two semi-stable saddlecenter equilibria at (±a, 0) and two stable states as center-type equilibria located in the region |x| > a. A pair of homoclinic orbits noted as type I connecting the origin (0, 0) to itself transversally cross two switching manifolds Σ ± = {(±a, y)|∀y ∈ R} as an isoenergetic separatrix to display large amplitude cross well motions outside it. There also exist a pair of symmetry connecting orbits noted as type II in the region |x| > a respectively homoclinic to two semi-stable saddle-center equilibria (±a, 0). The region between the type-I and type-II homoclinic orbits are full of periodic orbits traversing one of the switching manifolds Σ ± with any period and the area is getting larger and larger as k 1 increasing. There are also full of in-well periodic motions with arbitrary frequency located inside of the type-II homoclinic orbit. However, the area occupied by these periodic orbits satisfying |x| > a gradually shrinks to (±a, 0) as the increasing of k 1 to k * 1 . These periodic orbits with arbitrary frequency provide the geometrical structure of multifrequency internal resonance with the main structure, which has huge potential to increase the energy dissipation rates by internal resonance capture and obtain better energy dissipation effect as a passive EI-BNES. Furthermore, the instability of the saddle-center type equilibrium points at (±a, 0) will cause rapid increase in the velocity of the EI-BNES, which can be employed to acquire energy dissipation rates of the EI-BNES.
The Fig. 13(e) and Fig. 13(f) respectively display two representative piecewise bi-stable phase portraits for two given k 1 = 0.13 and k 1 = 0.16 satisfying k 1 > k * 1 . There exists a saddle (0, 0), two center equilibria (±a, 0) such that the EI-BNES shows stability at x = ±a. There are two piecewise-smooth homoclinic orbits connecting (0, 0) to itself that are symmetric about the y-axis. Moreover, the interiors of the homoclinic orbits are filled with periodic orbits with arbitrary frequency circling around the center (±a, 0). In the same way, these periodic orbits can also increase the energy dissipation rates of the EI-BNES because of the potential multi-frequency internal resonate capture with the main structure LO. Considering the actual physical model, the stable states of the mass m 2 occur at x = a and x = −a, which further illustrates that the effect of k 1 plays a major role when the stiffness coefficient satisfies k 1 > k * 1 even if the negative stiffness k 2 works.

Numerical simulations of the optimized EI-BNES
Next, according to the relationship k 1 > −0.5k 2 obtained from the above analysis, the values k 1 and k 2 will be modified to detect what conditions can make the EI-BNES have the better effects of energy dissipation. Since the energy dissipation in the section 2 is already shown to be better at k 1 = 0.2 when k 2 = −0.05, this set of parameters will be used to verify the effectiveness of the above proposed GMRA. Furthermore, k 2 = −0.15 and k 2 = −0.23 are respectively selected to obtain the optimized values of k 1 to match with k 2 such that the EI-BNES has better effects of energy dissipation. The thresholds of chaotic oscillations between A and k 1 for the EI-BNES are shown in Figure 14 when k 2 are selected for k 2 = −0.05, k 2 = −0.15, k 2 = −0.23, respectively. The Fig. 14 shows the threshold function A of k 1 can achieve a minimum for some assigned k 2 , implying that the Melnikov function itself will be more sensitive to certain spring coefficient ranges than to others. In particular, A respectively reach the minimum at k 1 = 0.18, k 1 = 0.31, k 1 = 0.41 when k 2 = −0.05, k 2 = −0.15, k 2 = −0.23 are assigned, implying that the Melnikov function will be most influenced by the specific spring coefficients k 1 such that the EI-BNES will be most susceptible to occur homoclinic bifurcations and therefore more likely to have crosswell oscillations, which is beneficial to the effects of energy dissipation. For k 2 = −0.05 the minimum A is obtained at k 1 = 0.18 which happens to be almost equal to the value k 1 = 0.2 known to achieve the best effect of energy dissipation. Therefore, the proposed GMRA and Melnikov function may become an effective analytical method for structural optimization of the EI-BNES in order to obtain higher energy dissipation effects. In additional, A respectively obtain the minimum at k 1 = 0.31 and k 1 = 0.41 when k 2 = −0.15 and k 2 = −0.23. The increasing of the k 2 -stiffness value will result in a decreased amplitude threshold for homoclinic bifurcation.
Next, the effects of energy dissipation for k 2 = −0.15, k 1 = 0.31 and k 2 = −0.23, k 1 = 0.41 is compared with the value k 2 = −0.05, k 1 = 0.18 and is simulated in Figure 15. It can be observed that the energy dissipation rates are high in the three cases. In the case of medium impulse excitations X, the effects of energy dissipation are better when k 1 = 0.41, k 2 = −0.23. In the cases of low and high impulse excitations X, the effects of energy dissipation are better when k 1 = 0.18, k 2 = −0.05. Hence, once the value k 2 is selected, we can use GMRA to detect the optimal range of k 1 . The achieve of the better effects of vibration reduction is also needed to consider to choose appropriate parameters of k 1 and k 2 according to the actual impulsive excitation X of the environment.

Conclusions
In this paper, a new highly efficient EI-BNES based on magnetic-elastic impacts with negative stiffness and bistability is proposed and realized by reasonable layout of permanent ring magnets and springs. Global dynamics and optimal design of the EI-BNES are studied emphatically in detail. A new GMRA is proposed to study global dynamics and homoclinic bifurcations of the reduced two-dimensional subsystem to explain the mechanism of TET and detect the chaos thresholds of the EI-BNES. A special type of saddle-center equilibrium points is also found in the non-smooth system of the EI-BNES and can be used to effectively increase the energy dissipation rates.Finally, the optimal design criterion of the EI-BNES for better dissipation performance is also first discussed and simulated based on the GMRA and the Melnikov function of the non-smooth systems. The effectiveness of the analytical GMRA is also verified by numerical simulations.