Vibration Analysis of a 1-DOF System Coupled with a Nonlinear Energy Sink with a Fractional Order Inerter

The fluid inerter described by the fractional derivative model is introduced into the traditional nonlinear energy sink (NES), which is called fractional-order NES in this paper. The slowly varying dynamic equation (SVDE) of the system coupled with fractional-order NES is obtained by the complex averaging method, in which the fractional derivative term is treated using the fractional Leibniz theorem. Then, the discriminants (Δ, Δ1, and Δ2) of the number of equilibrium points are derived. By using the variable substitution method, the characteristic equation for judging the stability is established. The results show: (1) the approximate SVDE is sufficient to reflect the slowly varying characteristics of the system, which shows that the mathematical treatment of the fractional derivative term is reliable; (2) the discriminant conditions (Δ1, Δ2) can accurately reflect the number of equilibrium points, and the corresponding range of nonlinear parameter κ can be calculated when the system has three equilibrium points. The expressions of Δ1, Δ2 are simpler than Δ, which is suitable for analysis and design parameters; (3) the stability discrimination methods of schemes 1 and 2 are accurate. Compared with scheme 2, scheme 1 is more prone to various responses, especially various strongly and weakly modulated responses. In scheme 2, the inertia effect of mass can be completely replaced by integer order inerter. Compared with integer order inerter, the introduction of fractional order inerter, whether in series or in parallel, means that the amplitude of the equilibrium point on the NES vibrator is smaller, but it is also for this reason that it is not easy to produce a modulated response with scheme 2, and the vibration suppression effect of the main structure is not good.


Introduction
When a fluid inerter works [1][2][3][4][5], it produces inertia force, usually accompanied by large parasitic damping force, which cannot be ignored. Therefore, the traditional mathematical model of an ideal inerter, such as a mechanical inerter [6][7][8], cannot describe the characteristics of a fluid inerter. The approximate model of its output force is mainly obtained by referring to the empirical formula. For example, the output force model of the fluid inerter proposed in [9][10][11] is obtained using the laminar flow empirical formula, and the output force model is established using the turbulence empirical formula in [4,5]. Some of these models consider laminar flow, and others consider the form of turbulence, but the actual situation may be both. There is still room for discussion on the model of fluid inerters. In reference [12], the author noted that describing the multiphase mechanical properties by simply adding the integer derivatives of different derivative orders is not satisfactory. From the perspective of a modeling fractional differential equation, Westerlund [13] proposed the use of a unified fractional model containing multiphase mechanical properties. In Obviously, when b n = 0, it is the vibration dynamic model coupled with the traditional NES. Therefore, the dynamic models of the two schemes can be obtained from Newton's second law, as shown in Equations (2) and (3).
S.1:  According to reference [32], when the generalized harmonic analysis is performed, the fractional derivative terms in vibration Equations (2) and (3) can be approximately equivalent in the following integer order as follows: In this paper, the nonlinear term is not limited to a small quantity, so system (5) can be a strongly nonlinear system. The complex averaging method (CA-X) [25,26] is used to derive the SVDE of system (5). Internal resonance of the system is one of the preconditions to ensure TET. Therefore, this paper studies the characteristics of the equilibrium point at 1:1 internal resonance i.e., Ω =1, and introduces complex variables as follows: ϕ j e iτ = z j + iz j , ϕ * j e −iτ = z j − iz j (j = 1, 2).
The following relations can be obtained by further derivation of the (7): , z j = ϕ j e iτ + ϕ * j e −iτ 2 , z j = .
To compare the slow variable model (12) with the z j obtained from the original system (5), the slow variable obtained from (12) is ϕ j , to be converted into a response in system (5). According to (7), z j can be determined by ϕ j , the formula is z j = imag(ϕ j e iτ ), z j = real(ϕ j e iτ ), j = 1,2. The fractional derivative term in the original system (5) is solved by using the improved Oustaloup's method [33]. Figure 2 (ε 1 = 0) and Figure 3 (ε 1 = 0.05) are displacement responses: (a) main structure (z 1 ), (b) NES (z 2 ), respectively (other parameters: ε = 0.05, µ = 1.9, ω 0 = 2.5, ξ s = 0.01, ξ n = 0.005, f e = 0.01, κ = 1). It can be seen from Figures 2 and 3 that the responses of the slow varying model are in good agreement with the numerical solution in an abbreviated period at the beginning, and has a small deviation from the numerical solution after a long period of time. In general, the SVDE (12) can reflect the change trend of the original system response, and can be used for the comparison of the later stable responses. varying model are in good agreement with the numerical solution in an abbreviated period at the beginning, and has a small deviation from the numerical solution after a long period of time. In general, the SVDE (12) can reflect the change trend of the original system response, and can be used for the comparison of the later stable responses.

The Discriminant Derivation of the Number of Equilibrium Points of S.1
To discuss the steady-state characteristics of the SVDE (S.1), it is difficult to calculate the equilibrium point of Equation (12) by directly making the derivative term 0. Therefore, the following variables are substituted for Equation (12): Substituting Equation (13) into Equation (12), the following relation can be obtained: If the derivative term in the above equation is zero, we obtain the equilibrium equation of the system as follows: ( ) 6  4  2  2  2  2  2  2  3  2  2  3  2  1   9 3 and

The Discriminant Derivation of the Number of Equilibrium Points of S.1
To discuss the steady-state characteristics of the SVDE (S.1), it is difficult to calculate the equilibrium point of Equation (12) by directly making the derivative term 0. Therefore, the following variables are substituted for Equation (12): Substituting Equation (13) into Equation (12), the following relation can be obtained: If the derivative term in the above equation is zero, we obtain the equilibrium equation of the system as follows: and where are a set of uncoupled nonlinear equations, and the equilibrium points of the system can be calculated from top to bottom. It can be seen from Equations (15a)-(15c) that (15b) and (15c) have only a single real root, so the number of equilibrium points of the system is determined by the number of real roots of (15a). The number of real roots of the polynomial can be determined by the Cardano discriminant. For equation a 2 x 3 + a 1 x + a 0 = 0, the Cardano discriminant is: When ∆ > 0, (15a) has a real root. When ∆ = 0, (15a) has three real roots, of which at least two are equal. When ∆ < 0, (15a) has three unequal real roots. To determine the number of real roots of (15a), it is necessary to replace them with variables, as shown below: Then, (15a) can be rewritten as: Therefore: Particularly, when ∆ = 0, the following relation can be obtained: Obviously, the number of κ satisfying ∆ = 0 is determined by the positive and negative of (21) as follows: ∆ 1 is the judgment condition for the quantity whose ∆ is equal to 0. When ∆ 1 = 0, ∆ = 0, the system has three equilibrium points, and in other positions when ∆ > 0, the system has only one equilibrium point. When ∆ 1 > 0, ∆ has two zero points, and there are three equilibrium points in the system between the two zero points, and there is only one equilibrium point at other positions. When ∆ 1 < 0, ∆ has no zero points, and there is only one equilibrium point at all positions.

The Discriminant Derivation of the Number of Equilibrium Points of S.2
The derivation method of S.2 is the same as S.1, so the derivation process is simplified, and the complex variables are substituted into (6) to obtain the SVDE of S.2 at 1:1 internal resonance, as shown below: Selecting the same variable as in the previous section to replace (Ψ 1 = ϕ 1 , Ψ 2 = ϕ 1 − ϕ 2 ), and the SVDE of S.2 is rewritten as: If the derivative term in (23) is 0, the equilibrium point equations of S.2 can be obtained as follows: 9 16 To judge the number of real roots of (24a), comparing (24a) with (15a) to find that the same variable can be used to replace. Substituting (17) into (24a), (24a) can be rewritten as follows: 9 16 The discriminant is obtained by substituting the coefficient of (25) into (16): Similarly, assuming ∆ = 0, the expression of κ can be simplified as follows: Therefore, the discriminant condition for the number of κ is: To understand the influence of various parameters on the number of equilibrium points when there is no inerter (i.e., b n = 0), according to Equations (2) and (3), the dynamic equations of the two schemes are the same. Therefore, taking S.2 as an example, this paper gives the discriminant ∆ 2 when b n = 0 by replacing parameters, as shown below: Observing the above formula, the value of ∆ 2 is determined by ξ n , ξ s , and ε, and the clear relationship is that the smaller ξ n is, the larger ∆ 2 is.

Stability Analysis Based on Lyapunov Theory
The equilibrium obtained from the equations derived in the previous section may be unstable. If the equilibrium point is unstable, the system may have a weakly modulated response (WMR) or SMR. As the response type of the system has a very important impact on the vibration suppression effect, it is necessary to analyze the stability of the equilibrium point of the system. For S.1, the Equation (14) is expanded near the equilibrium point, i.e., where Ψ 10 and Ψ 20 are the equilibrium points of the system, and δ 1 and δ 2 are the small increments near the corresponding equilibrium points. After algebraic operation, the linear differential equations about the increment can be obtained as follows: By substituting the following formula into (31): and separating the real part from the imaginary part, the following relationship can be obtained: where 20 in the above equation can be obtained from the following equation: Obviously, Equation (33) is a linear constant system, the stability can be directly determined by Lyapunov theory, and the characteristic factor is defined: λ i are the eigenvalues of the coefficient matrix on the right side of (33). According to Lyapunov theory, if Ξ < 0, the equilibrium point of the system is asymptotically stable; if Ξ = 0, the equilibrium point is critically stable; if Ξ > 0, the equilibrium point is unstable. If a pair of eigenvalues cross the imaginary axis when they reach a certain point with a certain parameter, then this point is also called the Hopf bifurcation point of the system.
The characteristic matrix of the system is as follows: Therefore, according to the following formula: its characteristic polynomial can be obtained as follows: Since the forms of a i (i = 1, 2, . . . 5) are too complex, and the eigenvalues in (38) can be easily calculated in MATLAB, the specific forms are not given here. For S.2, using the same method above, the system matrix A 2 of its incremental state equation can be derived from (23) as follows: . The stability of the equilibrium point of S.2 is determined by the eigenvalue of the characteristic equation corresponding to the system matrix A 2 .

Numerical Simulation and Discussion
Since the cubic stiffness of the NES is not much different from the linear stiffness of the actual structure, this section only studies the characteristics of the system equilibrium point when κ changes within a certain range, but the method used is not limited by parameter changes. The system parameters discussed in this section are normalized parameters, in which the constant parameters m s = 10, k s = 0.6, and f e = 0.01.

Influence of Parameters on the Number of Equilibrium Point
To find the parameter combination of three equilibrium points in the system without inerter (b n = 0), Figure 4a-c are drawn according to ∆ 2 (28), where Figure 4a shows that the larger c n , the smaller ∆ 2 , which is consistent with the previous prediction. Figure 4b shows that when m 2 increases, ∆ 2 first increases and then decreases, while the Figure 4c shows that the influence of c s is just opposite to m 2 . The greater the c s , ∆ 2 first decreases and then increases. To make ∆ 2 > 0, m 2 should choose a small value, and c s should choose a relatively large value.  Then, selecting the appropriate parameters, and considering the NES of paralle series inerter, the effects of parameters μ and bn on Δ1 and Δ2, respectively, are discu Considering the influence of inerter (bn ≠ 0), for S.1, Figure 5a-c are drawn according (21), as shown below. Figure 5a shows that the two surfaces almost coincide, and m2 has an insignif effect on Δ1, which is in line with the prediction of the curve when cs = 0.003 in Figu Then, selecting the appropriate parameters, and considering the NES of parallel and series inerter, the effects of parameters µ and b n on ∆ 1 and ∆ 2, respectively, are discussed. Considering the influence of inerter (b n = 0), for S.1, Figure 5a-c are drawn according to ∆ 1 (21), as shown below. Figure 5a shows that the two surfaces almost coincide, and m 2 has an insignificant effect on ∆ 1 , which is in line with the prediction of the curve when c s = 0.003 in Figure 4b. In Figure 5b, the change of µ directly changes the slope of ∆ 1 with b n , especially when µ = 5/3, the slope is 0. As shown in Figure 5c, b n has no effect on the µ range of ∆ 1 > 0, which is always between µ > 5/3 and 2. However, the larger b n is, the larger ∆ 1 is when µ > 5/3. In addition, when there is no inerter (b n = 0), the slope is 0, and ∆ 1 remains unchanged and less than 0, the system does not have multiple equilibrium points. Therefore, in the following discussion, µ should be selected between (5/3,2). Then, selecting the appropriate parameters, and considering the NES of paralle series inerter, the effects of parameters μ and bn on Δ1 and Δ2, respectively, are discu Considering the influence of inerter (bn ≠ 0), for S.1, Figure 5a-c are drawn according (21), as shown below. Figure 5a shows that the two surfaces almost coincide, and m2 has an insignif effect on Δ1, which is in line with the prediction of the curve when cs = 0.003 in Figu In Figure 5b, the change of μ directly changes the slope of Δ1 with bn, especially whe 5/3, the slope is 0. As shown in Figure 5c, bn has no effect on the μ range of Δ1 > 0, w is always between μ > 5/3 and 2. However, the larger bn is, the larger Δ1 is when μ > 5 addition, when there is no inerter (bn = 0), the slope is 0, and Δ1 remains unchanged less than 0, the system does not have multiple equilibrium points. Therefore, in following discussion, μ should be selected between (5/3,2). When S.2 adopt the same parameters as S.1, Δ2 is always negative, indicating t is not easy to make Δ2 > 0 in S.2. To highlight the role of the inerter, for S.2, m2 = 0, an influence of the inerter parameters on the value of Δ2 is analyzed. Under the prem m2 = 0, the surfaces of Δ2 versus μ and bn is drawn for different cs (cs = 0.1, cs = 0.3) in F As can be seen from Figure 6a,c, in the case of other parameter phase diagram larger the cs, the larger the range of Δ2 greater than zero. It can be seen from Figure 6 the change law of the curve is consistent with that of Figure 4b, which shows th inertance (bn) of the series inerter can nominally replace m2. However, for S.1 the ch law is different from that of the curve in Figure 4b, which shows that bn cannot si replace m2 when the inerter is connected in parallel.

Verification of Equilibrium Point Discriminant
To verify the positive and negative correspondence between the numb equilibrium points and Δ, and the prediction results of Δ1 and Δ2 in the previous se the curves of Δ and |Ψ2| with κ are drawn. As can be seen from Figure 6a,c, in the case of other parameter phase diagrams, the larger the c s, the larger the range of ∆ 2 greater than zero. It can be seen from Figure 6b that the change law of the curve is consistent with that of Figure 4b, which shows that the inertance (b n ) of the series inerter can nominally replace m 2 . However, for S.1 the change law is different from that of the curve in Figure 4b, which shows that b n cannot simply replace m 2 when the inerter is connected in parallel.

Verification of Equilibrium Point Discriminant
To verify the positive and negative correspondence between the number of equilibrium points and ∆, and the prediction results of ∆ 1 and ∆ 2 in the previous section, the curves of ∆ and |Ψ 2 | with κ are drawn.

Verification of Equilibrium Point Discriminant
To verify the positive and negative correspondence between the number of equilibrium points and Δ, and the prediction results of Δ1 and Δ2 in the previous section, the curves of Δ and |Ψ2| with κ are drawn.
For S.1, the values of m2 = 0.3 and m2 = 4 are selected, and the other parameters are consistent with those in Figure 5b   Whether in S.1 or S.2, the positive and negative of Δ correspond to the number of equilibrium points of the system one by one, and can also correspond to the positive and negative prediction of Δ1 or Δ2 in the previous section. Under the same parameters, the smaller μ, the smaller the amplitude of the equilibrium point, and the less likely it is to  Whether in S.1 or S.2, the positive and negative of Δ correspond to the number of equilibrium points of the system one by one, and can also correspond to the positive and negative prediction of Δ1 or Δ2 in the previous section. Under the same parameters, the smaller μ, the smaller the amplitude of the equilibrium point, and the less likely it is to switch the number of equilibrium points, indicating that the inertia and damping effects Whether in S.1 or S.2, the positive and negative of Δ correspond to the number of equilibrium points of the system one by one, and can also correspond to the positive and negative prediction of Δ1 or Δ2 in the previous section. Under the same parameters, the smaller μ, the smaller the amplitude of the equilibrium point, and the less likely it is to switch the number of equilibrium points, indicating that the inertia and damping effects Whether in S.1 or S.2, the positive and negative of ∆ correspond to the number of equilibrium points of the system one by one, and can also correspond to the positive and negative prediction of ∆ 1 or ∆ 2 in the previous section. Under the same parameters, the smaller µ, the smaller the amplitude of the equilibrium point, and the less likely it is to switch the number of equilibrium points, indicating that the inertia and damping effects of the inerter exist simultaneously when µ is between 1 and 2. Here, κ 1,2 = (0.22, 0.34) is calculated by substituting the parameters into (20) when µ = 1.8 in Figure 7b, and κ 1,2 = (0.77, 1.26) is calculated by substituting the parameters into (27) when µ = 2 in Figure 9b. Obviously, the calculated results of the two groups of parameters are consistent with the positions of the actual simulation curves. In addition, as predicted in the previous section, it is not easy to produce multiple equations with S.2. Particularly, when ∆ 1 = 0 or ∆ 2 = 0, the system is always only one equilibrium point. It can be considered that the system is just in a critical state.

Stability Discrimination of Equilibrium Points
In this section, the Lyapunov method in Section 3.3 is used to judge the stability of the equilibrium points of S.1 and S.2, and to verify whether the calculation results of the equilibrium point are consistent with those calculated by the numerical algorithm. Of course, this is meaningful only when the equilibrium point is asymptotically stable.
Taking S.1 as an example, the verification method for slowly varying parameters is established. According to (7), after the equilibrium point of system (14) is calculated by (15), the state of the original system (5) corresponding to the equilibrium point when t = 0 can be expressed as: A small disturbance (δ 1 , δ 2 , δ 3 , δ 4 ) is added to (40) and brought into the original system (5) as an initial condition. If the equilibrium point of the system is stable, the time response is stable after a period.
If Ψ 2 is expanded in the complex plane with N 2 e iθ , then |Ψ 2 | = N 2 . The steady-state amplitude |Ψ 2 | is obtained from (15a). Since Ψ 2 = ϕ 1 − ϕ 2 , the response of the original system can be converted to Ψ 2 , that is, Ψ 2 = −z 1 − z 2 + (z 1 − z 2 )i. The original system (5) is calculated by the fourth order Runge-Kutta (R-K) method, and its steady-state value N 2 can be obtained through conversion, as shown below: where < > (τ1 , τ 2) represents the weighted average in the (τ 1 , τ 2 ) period. Obviously, S.2 can also be treated in the same way. Adding a disturbance δ 1 = 0.01 to z 1 (0), the stability of the equilibrium point of the original system is analyzed. For S.1, the different stabilities of the equilibrium points are drawn, as shown in Figure 11, and the parameters of Figures 7b and 8b are used in Figure 11a,b, respectively. For S.2, the different stabilities of the equilibrium points are drawn, as shown in Figure 12, and the parameters of Figures 9b and 10b are used in Figure 12a,b, respectively. In Figures 11 and 12, the red squares represents the stable equilibrium points obtained by the R-K method and the black dots represents the stable equilibrium points obtained by (15a). The other represent unstable equilibrium points are also represented.
Firstly, from the above four pictures, the stable equilibrium points obtained by the R-K method are the same as that obtained by the calculation in (15a) or (24a), which verifies the accuracy of the analytical method. No matter which scheme, when the system has three equilibrium points, the switching of the stability always occurs near the turning point of the equilibrium point curve, and ∆ 1 and ∆ 2 proposed in this paper can only predict the turning point position, and the parameters are few and easy to analyze. original system is analyzed. For S.1, the different stabilities of the equilibrium poin drawn, as shown in Figure 11, and the parameters of Figures 7b and 8b are used in 11a,b, respectively. For S.2, the different stabilities of the equilibrium points are dra shown in Figure 12, and the parameters of Figures 9b and 10b are used in Figure  respectively. In Figures 11 and 12, the red squares represents the stable equilibrium obtained by the R-K method and the black dots represents the stable equilibrium obtained by (15a). The other represent unstable equilibrium points are also represe (a) (b)  Firstly, from the above four pictures, the stable equilibrium points obtained b R-K method are the same as that obtained by the calculation in (15a) or (24a), w verifies the accuracy of the analytical method. No matter which scheme, when the sy has three equilibrium points, the switching of the stability always occurs near the tur point of the equilibrium point curve, and Δ1 and Δ2 proposed in this paper can only pr the turning point position, and the parameters are few and easy to analyze. Figure 11a,b show that both schemes have three equilibrium points, but Figur shows two unstable equilibrium points with large amplitude and one stable equilib point with small amplitude, and Figure 11b shows an unstable equilibrium sandwiched between the two stable equilibrium points, indicating that there are different bifurcation phenomena in the system. When bn is calculated according to Δ or Δ2 = 0 is taken as the inertance of the system, the equilibrium point of the syste always stable with the change of κ. It is more difficult to produce multiple equilib points in S.2 ( Figure 12) than S.1 (Figure 11), and the range of unstable paramete relatively small.

The Time-Domain Response
Next, the time-domain responses are drawn to verify the accuracy of the a equilibrium stability analysis results, and the influence of equilibrium stability on response type when the slow varying system has multiple solutions is expl Particularly, to maintain the consistency between the steady-state amplitude N2 o system and |Ψ2| in the previous section, "z1-z2" is drawn instead of "z2".
In Figure 13a-d, the corresponding initial amplitude|Ψ2|(0) is (0.11, 0.41, 0.56, just between the three equilibrium points (|Ψ2|1,2,3 = 0.3109, 0.49, 0.79). The results s  Figure 11a,b show that both schemes have three equilibrium points, but Figure 11a shows two unstable equilibrium points with large amplitude and one stable equilibrium point with small amplitude, and Figure 11b shows an unstable equilibrium point sandwiched between the two stable equilibrium points, indicating that there are two different bifurcation phenomena in the system. When b n is calculated according to ∆ 1 = 0, or ∆ 2 = 0 is taken as the inertance of the system, the equilibrium point of the system is always stable with the change of κ. It is more difficult to produce multiple equilibrium points in S.2 ( Figure 12) than S.1 (Figure 11), and the range of unstable parameters is relatively small.

The Time-Domain Response
Next, the time-domain responses are drawn to verify the accuracy of the above equilibrium stability analysis results, and the influence of equilibrium stability on the response type when the slow varying system has multiple solutions is explored. Particularly, to maintain the consistency between the steady-state amplitude N 2 of the system and |Ψ 2 | in the previous section, "z 1 − z 2 " is drawn instead of "z 2 ". For S.1, we selected the parameters when k = 0.11 in Figure 11a (µ = 2 and µ = 1.9), and plotted the time-domain response of the original system ( Figure 13 (µ = 2) and Figure 14 (µ = 1.9)) under different initial interference (δ 1 ).
In Figure 13a-d, the corresponding initial amplitude|Ψ 2 |(0) is (0.11, 0.41, 0.56, 0.91), just between the three equilibrium points (|Ψ 2 | 1,2,3 = 0.3109, 0.49, 0.79). The results show that the response can be stabilized at the small stable equilibrium point only when it is between the small equilibrium point and the intermediate equilibrium point, while the system response is in the SMR when other disturbances occur, which is consistent with the conclusion when µ = 2, κ = 0.11 in Figure 11a. The four δ 1 are shown in Figure 14a-d, and the corresponding initial amplitudes (|Ψ 2 |(0)) are (0.07, 0.43, 0.73, 0.93) are just between the three equilibrium points (|Ψ 2 | 1,2,3 = 0.24, 0.63, 0.81). If the initial amplitude is less than the second largest equilibrium point, the system response asymptotically tends to the small stable equilibrium point after a long time of WMR (as shown in Figure 14a,b), while when it is greater than the second largest equilibrium point, the system response is in the SMR (as shown in Figure 14c,d). Then, the parameters (m2 = 4, κ = 0.1, μ = 2) and (m2 = 4, κ = 0.1, μ = 1.9) in Figure 11b are selected to draw the time-domain response of the original system (5) (as shown in Figures  15 and 16). Figure 15a-d show that no matter what kind of initial interference, the system response tends to be stable after a long time of WMR. When the initial value is between two stable equilibrium points, it is always stable at the small equilibrium point, otherwise, it is stable at the large stable equilibrium point. The four δ1 are shown in Figure 16a-d, and the corresponding initial amplitudes (|Ψ2|(0)) are (0.06, 0.33, 0.73, 0.93) are just between the three equilibrium points (|Ψ2|1,2,3 = 0.24, 0.67, 0.84). Basically consistent with the situation of μ = 2 (Figure 15), the system response can be stabilized at two stable equilibrium points, but the time to reach the stable equilibrium point is much shorter, and there is no obvious modulated process. Then, the parameters (m2 = 4, κ = 0.1, μ = 2) and (m2 = 4, κ = 0.1, μ = 1.9) in Figure 11b are selected to draw the time-domain response of the original system (5) (as shown in Figures  15 and 16). Figure 15a-d show that no matter what kind of initial interference, the system response tends to be stable after a long time of WMR. When the initial value is between two stable equilibrium points, it is always stable at the small equilibrium point, otherwise, it is stable at the large stable equilibrium point. The four δ1 are shown in Figure 16a-d, and the corresponding initial amplitudes (|Ψ2|(0)) are (0.06, 0.33, 0.73, 0.93) are just between the three equilibrium points (|Ψ2|1,2,3 = 0.24, 0.67, 0.84). Basically consistent with the situation of μ = 2 (Figure 15), the system response can be stabilized at two stable equilibrium points, but the time to reach the stable equilibrium point is much shorter, and there is no obvious modulated process. Then, the parameters (m 2 = 4, κ = 0.1, µ = 2) and (m 2 = 4, κ = 0.1, µ = 1.9) in Figure 11b are selected to draw the time-domain response of the original system (5) (as shown in Figures 15 and 16). Figure 15a-d show that no matter what kind of initial interference, the system response tends to be stable after a long time of WMR. When the initial value is between two stable equilibrium points, it is always stable at the small equilibrium point, otherwise, it is stable at the large stable equilibrium point. The four δ 1 are shown in Figure 16a-d, and the corresponding initial amplitudes (|Ψ 2 |(0)) are (0.06, 0.33, 0.73, 0.93) are just between the three equilibrium points (|Ψ 2 | 1,2,3 = 0.24, 0.67, 0.84). Basically consistent with the situation of µ = 2 (Figure 15), the system response can be stabilized at two stable equilibrium points, but the time to reach the stable equilibrium point is much shorter, and there is no obvious modulated process. For S.2: from Figure 10a,b in Section 4.2 and Figure 12a,b in Section 4.3, we can see that it is not easy to maintain multiple solutions, and the stability of the three equilibrium points is the same as that of S.1 (Figure 11b). To verify the special case of S.2, that is, whether the inerter can replace the mass function of traditional NES, we selected three groups of parameter schemes for comparison: T1 (μ = 2, m2 = 0, bn = 0.4), T2 (bn = 0, m2 = 0.4), T3 (μ = 1.9, m2 = 0, bn = 0.4) and the other parameters are selected when κ = 1.1 in Figure 12a. Figure 17a,b, and c show the phase diagrams of the three parameter schemes of the original system (6). The initial interference selection method is the same as that of S.1, and the four initial interferences (δ1) are (−0.01, 0.0, 8, 0.14, 0.3). Figure 18a,b show the displacement responses (z1, z1-z2) of T1 and T2 when the interference is −0.01 and 0.14, respectively. Figure 18c,d show the displacement (z1, z1-z2) response of T3 and T2 when the interference is −0.01 and 0.14, respectively. For S.2: from Figure 10a,b in Section 4.2 and Figure 12a,b in Section 4.3, we can see that it is not easy to maintain multiple solutions, and the stability of the three equilibrium points is the same as that of S.1 (Figure 11b). To verify the special case of S.2, that is, whether the inerter can replace the mass function of traditional NES, we selected three groups of parameter schemes for comparison: T1 (μ = 2, m2 = 0, bn = 0.4), T2 (bn = 0, m2 = 0.4), T3 (μ = 1.9, m2 = 0, bn = 0.4) and the other parameters are selected when κ = 1.1 in Figure 12a. Figure 17a,b, and c show the phase diagrams of the three parameter schemes of the original system (6). The initial interference selection method is the same as that of S.1, and the four initial interferences (δ1) are (−0.01, 0.0, 8, 0.14, 0.3). Figure 18a,b show the displacement responses (z1, z1-z2) of T1 and T2 when the interference is −0.01 and 0.14, respectively. Figure 18c,d show the displacement (z1, z1-z2) response of T3 and T2 when the interference is −0.01 and 0.14, respectively.  Figure 12a,b in Section 4.3, we can see that it is not easy to maintain multiple solutions, and the stability of the three equilibrium points is the same as that of S.1 (Figure 11b). To verify the special case of S.2, that is, whether the inerter can replace the mass function of traditional NES, we selected three groups of parameter schemes for comparison: T 1 (µ = 2, m 2 = 0, b n = 0.4), T 2 (b n = 0, m 2 = 0.4), T 3 (µ = 1.9, m 2 = 0, b n = 0.4) and the other parameters are selected when κ = 1.1 in Figure 12a. Figure 17a,b, and c show the phase diagrams of the three parameter schemes of the original system (6). The initial interference selection method is the same as that of S.1, and the four initial interferences (δ 1 ) are (−0.01, 0.0, 8, 0.14, 0.3). Figure 18a,b show the displacement responses (z 1 , z 1 − z 2 ) of T 1 and T 2 when the interference is −0.01 and 0.14, respectively. Figure 18c,d show the displacement (z 1 , z 1 − z 2 ) response of T 3 and T 2 when the interference is −0.01 and 0.14, respectively.
It can be seen from Figure 18c,d that the time-domain responses of T 3 and T 2 are basically coincidental, which shows that the b n of S.2 can replace the role of m 2 equivalently when the integer order inerter is connected in series. From Figures 17a and 18a,b, T 1 and T 2 can be stabilized at small stable equilibrium points no matter how large the interference (δ 1 ) is. As can be seen from Figure 18b, when the initial amplitude is above the small equilibrium point, the amplitude of NES jumps from the large equilibrium point to the small equilibrium point due to the fractional order inerter in series in S.2, resulting in the inability to effectively suppress the displacement of the main structure. At this time, the effect of vibration suppression is not as good as that of series integer order inerter. It can be seen from Figure 18c,d that the time-domain responses of T3 and T2 are basically coincidental, which shows that the bn of S.2 can replace the role of m2 equivalently when the integer order inerter is connected in series. From Figures 17a and 18a,b, T1 and T2 can be stabilized at small stable equilibrium points no matter how large the interference (δ1) is. As can be seen from Figure 18b, when the initial amplitude is above the small equilibrium point, the amplitude of NES jumps from the large equilibrium point to the small equilibrium point due to the fractional order inerter in series in S.2, resulting in the inability to effectively suppress the displacement of the main structure. At this time, the effect of vibration suppression is not as good as that of series integer order inerter.
The analysis of different situations of the above two schemes shows that the stability discrimination in the previous section is accurate. No matter the scheme, μ has a great impact on response, it is sensitive to the initial value when only the inertia effect (i.e., μ = 2) exists, and it is easy for various modulated responses to appear in the time-domain response. Under the same conditions, S.1 is more sensitive to the nonlinear parameter κ than S.2, and it is easy to produce different system responses.

Conclusions
In this paper, the inerter described by fractional derivative is used in traditional NES for the first time. Combined with the CA-X method and fractional properties, especially fractional Leibniz theorem, the difficulty of deriving SVDE from fractional nonlinear dynamic equations is solved. Through the verifications of analytical and numerical methods, the accuracy of SVDE ensures that the analysis is reliable.
The discriminants Δ1 (S.1) and Δ2 (S.2) for the number of equilibrium points are proposed. Compared with Δ, their expressions are simpler, which can directly and accurately predict the number of equilibrium points of the equilibrium point equation, It can be seen from Figure 18c,d that the time-domain responses of T3 and T2 are basically coincidental, which shows that the bn of S.2 can replace the role of m2 equivalently when the integer order inerter is connected in series. From Figures 17a and 18a,b, T1 and T2 can be stabilized at small stable equilibrium points no matter how large the interference (δ1) is. As can be seen from Figure 18b, when the initial amplitude is above the small equilibrium point, the amplitude of NES jumps from the large equilibrium point to the small equilibrium point due to the fractional order inerter in series in S.2, resulting in the inability to effectively suppress the displacement of the main structure. At this time, the effect of vibration suppression is not as good as that of series integer order inerter.
The analysis of different situations of the above two schemes shows that the stability discrimination in the previous section is accurate. No matter the scheme, μ has a great impact on response, it is sensitive to the initial value when only the inertia effect (i.e., μ = 2) exists, and it is easy for various modulated responses to appear in the time-domain response. Under the same conditions, S.1 is more sensitive to the nonlinear parameter κ than S.2, and it is easy to produce different system responses.

Conclusions
In this paper, the inerter described by fractional derivative is used in traditional NES for the first time. Combined with the CA-X method and fractional properties, especially fractional Leibniz theorem, the difficulty of deriving SVDE from fractional nonlinear dynamic equations is solved. Through the verifications of analytical and numerical methods, the accuracy of SVDE ensures that the analysis is reliable.
The discriminants Δ1 (S.1) and Δ2 (S.2) for the number of equilibrium points are proposed. Compared with Δ, their expressions are simpler, which can directly and accurately predict the number of equilibrium points of the equilibrium point equation, The analysis of different situations of the above two schemes shows that the stability discrimination in the previous section is accurate. No matter the scheme, µ has a great impact on response, it is sensitive to the initial value when only the inertia effect (i.e., µ = 2) exists, and it is easy for various modulated responses to appear in the time-domain response. Under the same conditions, S.1 is more sensitive to the nonlinear parameter κ than S.2, and it is easy to produce different system responses.

Conclusions
In this paper, the inerter described by fractional derivative is used in traditional NES for the first time. Combined with the CA-X method and fractional properties, especially fractional Leibniz theorem, the difficulty of deriving SVDE from fractional nonlinear dynamic equations is solved. Through the verifications of analytical and numerical methods, the accuracy of SVDE ensures that the analysis is reliable.
The discriminants ∆ 1 (S.1) and ∆ 2 (S.2) for the number of equilibrium points are proposed. Compared with ∆, their expressions are simpler, which can directly and accurately predict the number of equilibrium points of the equilibrium point equation, and the range of κ corresponding to the system with three equilibrium points can be calculated by (20) and (27).
The verifications in Figures 11-18 show that the Lyapunov theory given in this paper is accurate in judging the stability of the equilibrium point. For S.2, when m 2 = 0, the introduction of the inerter does not increase the complexity of the structure, and the inerter can not only provide inertia, but also provide damping. S.1 has two unstable equilibrium points or two stable equilibrium points at the same time, so S.1 is more prone to various responses, especially various WMR and SMR. Compared with integer order inerter, the introduction of fractional order inerter, whether in series or in parallel, means that the amplitude of the equilibrium point on the NES vibrator is smaller, but it is also for this reason that it is not easy to produced modulated responses in S.2, and the vibration suppression effect of the main structure is not good.
Author Contributions: Writing-original draft preparation, Y.C.; writing-review and editing, Y.C., Y.T. and J.X.; supervision, X.X.; supervision, N.C. All authors have read and agreed to the published version of the manuscript.