Analysis of Nonlinear Dynamic Characteristics of a Mechanical-Electromagnetic Vibration System with Rubbing

: In this study, a rotor-bearing-runner system (RBRS) considering multiple nonlinear factors is established, and the complex nonlinear dynamic behavior of the coupling system is studied. The e ﬀ ects of excitation current, radial sti ﬀ ness, and friction coe ﬃ cient on dynamic characteristics are analyzed by numerical simulation. The research results show that the dynamic properties of the coupling system caused by di ﬀ erent nonlinear factors are interactional. With the changes of di ﬀ erent parameters, the RBRS presents multiple motion states, including periodic-n, quasi-periodic, and chaotic motion. The increase of the excitation current I j has a certain inhibitory e ﬀ ect on the response amplitude of the system and makes the motion state of the system more complex, the chaotic motion wider, and the jump discontinuity enhanced. With the increase of radial sti ﬀ ness k r , the motion complexity of the coupling system increases, the chaotic region increases, the response amplitude increases, and the vibration intensity increases. With the increase of the friction coe ﬃ cient µ , the chaotic region increases ﬁrst and decreases, the di ﬀ erent motions alternate frequently, and the response amplitude gradually increases. This study can not only help to understand the dynamic characteristics of RBRS, but also help the stable operation of the generator set.


Introduction
A large-scale rotating generator set is a complex system involving many factors, such as mechanical and electromagnetic. The vibration failure during its operation seriously threatens the safe and stable operation of large rotating machinery. The unbalanced magnetic pull force of the rotor-stator and the rubbing of the bearing-shaft system are two important vibration fault sources of the generator set. In order to improve the vibration characteristics of the system, it is very important to establish a precise dynamic model and analyze the dynamic characteristics of the generator set. Researchers conducted a large number of analytical studies on this. Ma et al. [1] and Wen [2] gave a general introduction to the nonlinear vibration system of rotating machinery such as hydraulic turbines. Chu et al. [3] designed a rotor system test bench that can simulate rotor-stator friction and analyzed the harmonic components under different conditions. Patel et al. [4] discussed the torsional and bending vibrations of the rotor system due to friction and cracking. In Cao et al. [5]'s paper, nonlinear dynamics of the fractionally damped rub-impact rotor system are investigated. The rubbing forces consist of a radial elastic force and a tangential Coulomb friction force. The influence of speed, damping coefficient, and eccentricity is discussed. Focusing on the local rubbing problem of the generator rotor, Huang et al. [6] established a dynamic model of rotor misalignment and frictional rubbing coupling, and analyzed the dynamic In order to study the dynamic behavior of RBRS effectively, the finite element model of double rotor-bearing system of hydraulic turbine is modeled and simplified as four lumped mass points according to the structural characteristics of the system. The rotor and the runner are located at node 1 and node 4, respectively, and the combination bearing and turbine guide bearing are located at node 2 and node 3, respectively. The schematic diagram of RBRS is shown in Figure 1. The rotation center coordinates of the turbine rotor and the runner are O1(x1, y1) and O4(x4, y4), and the centroid coordinates of the rotor and the runner are G1(xg1, yg1) and G4(xg4, yg4).
The dynamic equation of RBRS with 16 degrees of freedom can be deduced as follows: ( ) e w r ρ + + +   Mq+ C+G q+Kq= F F F F (1) where M, C, G, K, q are mass matrix, damping matrix, gyroscopic matrix, stiffness matrix, and displacement vector of RBRS, respectively .   T   1  1  2  2  3  3  4  4  1  1  2  2  3  3  4  4 , , ,  (2) where xi, yi (i = 1,2,3,4) and θxi, θyi (i = 1,2,3,4) are the displacements of the rotor, combination bearing, turbine guide bearing, and runner in x and y directions and angles of orientation associated with the x and y axes, respectively. Suppose that RBRS is composed of an elastic shaft element and the generalized coordinates of the shaft element are the displacements of nodes at both ends. The rotor system in this study only considers bending deformation and torsional deformation in horizontal and vertical directions, ignoring axial deformation.  q + (C + G) where M, C, G, K, q are mass matrix, damping matrix, gyroscopic matrix, stiffness matrix, and displacement vector of RBRS, respectively. q = x 1 , θ y1 , x 2 , θ y2 , x 3 , θ y3 , x 4 , θ y4 , y 1 , θ x1 , y 2 , θ x2 , y 3 , θ x3 , y 4 , θ x4 T (2) where x i , y i (i = 1,2,3,4) and θ xi , θ yi (i = 1,2,3,4) are the displacements of the rotor, combination bearing, turbine guide bearing, and runner in x and y directions and angles of orientation associated with the x and y axes, respectively. Suppose that RBRS is composed of an elastic shaft element and the generalized coordinates of the shaft element are the displacements of nodes at both ends. The rotor system in this study only considers bending deformation and torsional deformation in horizontal and vertical directions, ignoring axial deformation.
Appl. Sci. 2019, 9,4612 4 of 23 where M sT , M sR ,G s (G s = ωJ s ), and K s represent moving inertia matrix of the unit, rotational inertia matrix of the unit, rotation matrix of the unit, and stiffness matrix of the unit, respectively. µ, l, and EI are the quality of the unit length, the length of the unit, and bending stiffness. 0 diag 0, 0, , 0, , 0, 0, 0, 0, 0, , 0, , 0, 0, 0 where kz2x, kz3x, kz2y, and kz3y are the stiffness of the combination bearing and turbine guide bearing in x and y directions, respectively.
The total damping matrix (Rayleigh damping form is applied to determine the viscous part) can be obtained by the following formula [30]: Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 24   2  2  2   2  2   36  3  36 3  3  4  3  36  3  36  3  120  3  3 4   sR   l  l  l  l  l  l  r  l  l  l  l 2   36  3  36 3  3  4  3  36  3  36  3  60  3  3 where MsT, MsR, Gs (Gs = ωJs), and Ks represent moving inertia matrix of the unit, rotational inertia matrix of the unit, rotation matrix of the unit, and stiffness matrix of the unit, respectively. μ, l, and EI are the quality of the unit length, the length of the unit, and bending stiffness.
The total damping matrix (Rayleigh damping form is applied to determine the viscous part) can be obtained by the following formula [30]: where k z2x , k z3x , k z2y , and k z3y are the stiffness of the combination bearing and turbine guide bearing in x and y directions, respectively.
The total damping matrix (Rayleigh damping form is applied to determine the viscous part) can be obtained by the following formula [30]: where ω 1 and ω 2 are the first and second natural frequencies of the system, respectively. ξ 1 and ξ 2 are the first and second modal damping ratios, respectively.
where c z2x , c z3x , c z2y , and c z3y are the damping of the combination bearing and turbine guide bearing in x and y directions, respectively.
where m 1 ρ 1 and m 4 ρ 4 denote eccentricity of unbalance mass of the generator rotor and the turbine runner. F e = F ex , 0, 0, 0, 0, 0, 0, 0, F ey , 0, 0, 0, 0, 0, 0, 0 T where F ex and F ey are nonlinear electromagnetic force of the rotor in x and y directions. Nonlinear electromagnetic force can be expressed as: where k n , ε and θ r represent electromagnetic stiffness, relative eccentric, and rotation angle of the rotor. α and δ 0 are variation coefficient of the electromagnetic stiffness and average air gap length. The functions k n , ε, and θ r are respectively given in Equations (20)- (22): where µ 0 , r, L, K j , and I j are permeability of vacuum, rotor radius, rotor height, coefficient of air gap fundamental wave magnetomotive force, and excitation current of the generator, respectively.
where F wx and F wy are the hydraulic force of the runner in x and y directions. Hydraulic force can be expressed as: where K w is coefficient of the hydraulic force.
F r = 0, 0, F rx , 0, 0, 0, 0, 0, 0, 0, F ry , 0, 0, 0, 0, 0 where F rx and F ry are nonlinear rubbing force of the combination bearing in x and y directions. Nonlinear rubbing force can be expressed as: where r represents vibration amplitude of the rotor. k r , δ, and f are radial stiffness of the stator, gap, and friction coefficient of rubbing. Where the functions r is given in Equation (27): where ξ is the relative position of rotor centroid and stator centroid.

Nonlinear Dynamic Analysis of the RBRS
In the previous section, the dynamics model of the rotor-bearing-runner system (RBRS) is introduced and the system dynamics equation is obtained. In this section, the influence of different system parameters on nonlinear systems under the coupling of mechanical unbalance force, nonlinear magnetic pull force, hydraulic unbalance force, and nonlinear rubbing force is studied, through Bifurcation diagram, 3-D frequency spectrum, Waveform, Frequency, Axis orbit, and Poincaré map. The parameters for the coupling system are presented in Table 1. In addition, the detailed simulation condition schematic is displayed in Figure 2.

Model Verification
In order to verify the correctness of the non-linear dynamic model built in this paper, the accuracy of the model was verified by comparing it with the models built in the existing literature. In simulation, the parameters given in Table 1 were used in both models. As shown in Figure 3, the three-dimensional spectrum of the model in the Y direction of the rotor is given. Figure 3a is the threedimensional spectrum of the model in reference [31]. Figure 3b is the calculation result of the model established in this paper. By comparison, it can be found that when the rotating speed is small, the frequency components of the two models are mainly one times frequency conversion fn1, and no rubbing occurs in the shafting. With the increase of rotational speed, the main components of Figure  3a are the frequency multiplication components fn1, 2fn1, 3fn1, and the frequency demultiplication component is not obvious. In addition, the amplitude of one times frequency conversion fn1 shows a fluctuating growth, while the amplitude of other frequency multiplication components increases first and then decreases. In Figure 3b, besides the frequency conversion and frequency multiplication components, there are also obvious frequency demultiplication components, such as fn1/2, 3fn1/2, 5fn1/2, 7fn1/2, etc. The amplitude of each frequency component varies nonlinearly, and at a certain

Model Verification
In order to verify the correctness of the non-linear dynamic model built in this paper, the accuracy of the model was verified by comparing it with the models built in the existing literature. In simulation, the parameters given in Table 1 were used in both models. As shown in Figure 3, the three-dimensional spectrum of the model in the Y direction of the rotor is given. Figure 3a is the three-dimensional spectrum of the model in reference [31]. Figure 3b is the calculation result of the model established in this paper. By comparison, it can be found that when the rotating speed is small, the frequency components of the two models are mainly one times frequency conversion f n1 , and no rubbing occurs in the shafting. With the increase of rotational speed, the main components of Figure 3a are the frequency multiplication components f n1 , 2f n1 , 3f n1 , and the frequency demultiplication component is not obvious. In addition, the amplitude of one times frequency conversion f n1 shows a fluctuating growth, while the amplitude of other frequency multiplication components increases first and then decreases. In Figure 3b, besides the frequency conversion and frequency multiplication components, there are also obvious frequency demultiplication components, such as f n1 /2, 3f n1 /2, 5f n1 /2, 7f n1 /2, etc. The amplitude of each frequency component varies nonlinearly, and at a certain speed, there is an obvious continuous spectrum component in the three-dimensional spectrum. At this time, the system is in a chaotic state, but not in Figure 3a. When the speed reaches 71 rad/s, there are only obvious f n1 and 2f n1 components in Figure 3a. Besides frequency multiplication components, Figure 3b has more obvious frequency demultiplication components.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 24 speed, there is an obvious continuous spectrum component in the three-dimensional spectrum. At this time, the system is in a chaotic state, but not in Figure 3a. When the speed reaches 71 rad/s, there are only obvious fn1 and 2fn1 components in Figure 3a. Besides frequency multiplication components, Figure 3b has more obvious frequency demultiplication components. In order to compare the frequency components of the two models more clearly, Table 2 gives the comparison of frequency characteristics in three-dimensional spectrum. By comparing Figure 3 with Table 2, it can be seen that the model established in this paper takes into account the electromagneticrubbing coupling, non-linear electromagnetic force, hydraulic unbalance force, and other factors, which can better reflect the non-linear vibration characteristics of multi-physical field coupling. Note: "0" represents no rubbing; "1" represents rubbing.

Influence of Excitation Current on Dynamic Characteristics of RBRS at Variable Speed
Nonlinear electromagnetic force has an important influence on the inherent characteristics and dynamic response of RBRS, which is one of the important reasons for the nonlinear characteristics of RBRS. In order to study the influence of different excitation currents of the generator Ij on the dynamic characteristics of RBRS, Figures 4 and 5, respectively, give the bifurcation diagram and 3-D frequency spectrum of RBRS with rotational speed ω (ω = 7.0~71 rad/s) as control parameter at Ij = 1000 A.
When the rotational speed is small, the motion state of RBRS is periodic motion, and the response amplitude is relatively stable. At the same time, the frequency component is not very obvious, mainly based on one times the frequency conversion fn. Figure 6 shows the characteristics of periodic motion at ω = 19.4 rad/s, through waveform, frequency, axis orbit, and Poincaré map. With the increase of the rotational speed, the system passes through saddle-node bifurcation at ω = 20 rad/s and enters chaotic motion. In the range of ω ∈ [20.0,59.8] rad/s, the system mainly exhibits chaotic motion, except for a few points where the vibration system presents periodic-n motion at ω = 24.4 rad/s, ω = In order to compare the frequency components of the two models more clearly, Table 2 gives the comparison of frequency characteristics in three-dimensional spectrum. By comparing Figure 3 with Table 2, it can be seen that the model established in this paper takes into account the electromagnetic-rubbing coupling, non-linear electromagnetic force, hydraulic unbalance force, and other factors, which can better reflect the non-linear vibration characteristics of multi-physical field coupling.

Influence of Excitation Current on Dynamic Characteristics of RBRS at Variable Speed
Nonlinear electromagnetic force has an important influence on the inherent characteristics and dynamic response of RBRS, which is one of the important reasons for the nonlinear characteristics of RBRS. In order to study the influence of different excitation currents of the generator I j on the dynamic characteristics of RBRS, Figures 4 and 5, respectively, give the bifurcation diagram and 3-D frequency spectrum of RBRS with rotational speed ω (ω = 7.0~71 rad/s) as control parameter at I j = 1000 A. fn/3 is the main frequency component, and the amplitude gradually increases and is significantly larger than other frequency components. The vibration responses of the RBRS at ω = 59.8 rad/s are displayed in Figure 8. As can be seen from the figure, the coupling system is quasi-periodic motion. The main frequency components are fn, 2fn, and 3fn. Meanwhile, the attractors form a closed circle in the Poincaré map.   larger than other frequency components. The vibration responses of the RBRS at ω = 59.8 rad/s are displayed in Figure 8. As can be seen from the figure, the coupling system is quasi-periodic motion. The main frequency components are fn, 2fn, and 3fn. Meanwhile, the attractors form a closed circle in the Poincaré map.   When the rotational speed is small, the motion state of RBRS is periodic motion, and the response amplitude is relatively stable. At the same time, the frequency component is not very obvious, mainly based on one times the frequency conversion f n . Figure 6 shows the characteristics of periodic motion at ω = 19.4 rad/s, through waveform, frequency, axis orbit, and Poincaré map. With the increase of the rotational speed, the system passes through saddle-node bifurcation at ω = 20 rad/s and enters chaotic motion. In the range of ω ∈ [20.0,59.8] rad/s, the system mainly exhibits chaotic motion, except for a few points where the vibration system presents periodic-n motion at ω = 24.4 rad/s, ω = 24.6 rad/s, ω = 42.6 rad/s, and ω = 42.8 rad/s, and has different bifurcations such as saddle-node bifurcation, period-doubling bifurcation, inverse period-doubling bifurcation; in the three-dimensional spectrum, the continuous and discrete spectral components alternate and the amplitude changes nonlinearly. The vibration responses of the RBRS at ω = 33.6 rad/s are given in Figure 7. It can be seen from the figure that the amplitude of the system vibration response increases and changes aperiodically; in the frequency domain, besides the frequency conversion and frequency multiplication components, the continuous frequency components also appear, and the axis orbit of the system is disorderedly.
On the corresponding Poincaré map, the attractors are irregular discrete points, and the system exhibits chaotic motion. With the speed increasing further, when ω = 59.6 rad/s, the system moves from a short quasi-periodic motion to a stable periodic motion via the Hopf bifurcation; except for the frequency multiplication components (f n , 2f n , 3f n ,...) in the 3-D frequency spectrum, there are also frequency demultiplication components (f n /3, 2f n /3, 4f n /3, 5f n /3,...), in which f n /3 is the main frequency component, and the amplitude gradually increases and is significantly larger than other frequency components. The vibration responses of the RBRS at ω = 59.8 rad/s are displayed in Figure 8. As can be seen from the figure, the coupling system is quasi-periodic motion. The main frequency components are f n , 2f n , and 3f n . Meanwhile, the attractors form a closed circle in the Poincaré map.   In order to further comprehensively study the effects of different excitation current of the generator Ij on RBRS, Figures 9-12, respectively, give the bifurcation diagram and 3-D frequency spectrum of RBRS at Ij = 700 A and Ij = 1200 A. By comparing the bifurcation diagrams and the 3-D frequency spectrum, it can be seen that the system shows different vibration forms and different bifurcation paths under various excitation currents. When the excitation current of the generator Ij is small, the dynamic behavior of the system is relatively simple. When ω < 20.2 rad/s, the vibration system behaves as a periodic motion, and the frequency components is not only the frequency conversion fn, but also the unobvious frequency multiplication component. As the rotational speed increases, the system transits from short-term periodic motion, chaotic motion, and quasi-periodic motion to periodic motion, successively undergoing saddle-node bifurcation, period-doubling bifurcation, inverse period-doubling bifurcation, and Hopf bifurcation. Frequency multiplication is the main frequency component in the 3-D frequency spectrum, and the peak value appears at 7fn when ω = 21.0 rad/s. When the rotational speed increases to ω = 54.8 rad/s, the system enters the chaotic motion and transits to stable periodic motion through the addle-node bifurcation, and the amplitude of response at ω = 68.0 rad/s appears obvious jumping discontinuity. Comparing Figures  4 and 5, Figures 9 and 10, and Figures 11 and 12 shows that the dynamic behavior of the system is   In order to further comprehensively study the effects of different excitation current of the generator Ij on RBRS, Figures 9-12, respectively, give the bifurcation diagram and 3-D frequency spectrum of RBRS at Ij = 700 A and Ij = 1200 A. By comparing the bifurcation diagrams and the 3-D frequency spectrum, it can be seen that the system shows different vibration forms and different bifurcation paths under various excitation currents. When the excitation current of the generator Ij is small, the dynamic behavior of the system is relatively simple. When ω < 20.2 rad/s, the vibration system behaves as a periodic motion, and the frequency components is not only the frequency conversion fn, but also the unobvious frequency multiplication component. As the rotational speed increases, the system transits from short-term periodic motion, chaotic motion, and quasi-periodic motion to periodic motion, successively undergoing saddle-node bifurcation, period-doubling bifurcation, inverse period-doubling bifurcation, and Hopf bifurcation. Frequency multiplication is the main frequency component in the 3-D frequency spectrum, and the peak value appears at 7fn when ω = 21.0 rad/s. When the rotational speed increases to ω = 54.8 rad/s, the system enters the chaotic motion and transits to stable periodic motion through the addle-node bifurcation, and the amplitude of response at ω = 68.0 rad/s appears obvious jumping discontinuity. Comparing Figures  4 and 5, Figures 9 and 10, and Figures 11 and 12 shows that the dynamic behavior of the system is   In order to further comprehensively study the effects of different excitation current of the generator Ij on RBRS, Figures 9-12, respectively, give the bifurcation diagram and 3-D frequency spectrum of RBRS at Ij = 700 A and Ij = 1200 A. By comparing the bifurcation diagrams and the 3-D frequency spectrum, it can be seen that the system shows different vibration forms and different bifurcation paths under various excitation currents. When the excitation current of the generator Ij is small, the dynamic behavior of the system is relatively simple. When ω < 20.2 rad/s, the vibration system behaves as a periodic motion, and the frequency components is not only the frequency conversion fn, but also the unobvious frequency multiplication component. As the rotational speed increases, the system transits from short-term periodic motion, chaotic motion, and quasi-periodic motion to periodic motion, successively undergoing saddle-node bifurcation, period-doubling bifurcation, inverse period-doubling bifurcation, and Hopf bifurcation. Frequency multiplication is the main frequency component in the 3-D frequency spectrum, and the peak value appears at 7fn when ω = 21.0 rad/s. When the rotational speed increases to ω = 54.8 rad/s, the system enters the chaotic motion and transits to stable periodic motion through the addle-node bifurcation, and the amplitude of response at ω = 68.0 rad/s appears obvious jumping discontinuity. Comparing Figures  4 and 5, Figures 9 and 10, and Figures 11 and 12 shows that the dynamic behavior of the system is more complicated, and the discontinuous jump is obviously increased. As the rotational speed In order to further comprehensively study the effects of different excitation current of the generator I j on RBRS, Figures 9-12, respectively, give the bifurcation diagram and 3-D frequency spectrum of RBRS at I j = 700 A and I j = 1200 A. By comparing the bifurcation diagrams and the 3-D frequency spectrum, it can be seen that the system shows different vibration forms and different bifurcation paths under various excitation currents. When the excitation current of the generator I j is small, the dynamic behavior of the system is relatively simple. When ω < 20.2 rad/s, the vibration system behaves as a periodic motion, and the frequency components is not only the frequency conversion f n , but also the unobvious frequency multiplication component. As the rotational speed increases, the system transits from short-term periodic motion, chaotic motion, and quasi-periodic motion to periodic motion, successively undergoing saddle-node bifurcation, period-doubling bifurcation, inverse period-doubling bifurcation, and Hopf bifurcation. Frequency multiplication is the main frequency component in the 3-D frequency spectrum, and the peak value appears at 7f n when ω = 21.0 rad/s. When the rotational speed increases to ω = 54.8 rad/s, the system enters the chaotic motion and transits to stable periodic motion through the addle-node bifurcation, and the amplitude of response at ω = 68.0 rad/s appears obvious jumping discontinuity. Comparing Figures 4 and 5, Figures 9 and 10, and Figures 11 and 12 shows that the dynamic behavior of the system is more complicated, and the discontinuous jump is obviously increased. As the rotational speed increases, the chaotic region becomes wider. In the range of ω ∈ [20.0,22.4] rad/s, except for the stable period-doubling motion in a small region, the system is of continuous chaotic motion or intermittency chaos motion. When ω = 29.0 rad/s, the system enters the chaotic motion through the saddle-node bifurcation, and in the range of ω > 28.8 rad/s, the vibration system exhibits strong nonlinear characteristics of periodic motion and chaotic motion frequently alternating. With the increase of the rotational speed, different bifurcation paths, such as saddle-node bifurcation, period-doubling bifurcation, and inverse period-doubling bifurcation appear in the system. The corresponding vibration amplitude increases first, then decreases and then increases nonlinearly. In order to specify the main frequency components of the system under different excitation currents of the generator I j , the main frequency components (where f n is the frequency conversion) in different speed ranges are listed in Table 3.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 24 increases, the chaotic region becomes wider. In the range of ω ∈ [20.0,22.4] rad/s, except for the stable period-doubling motion in a small region, the system is of continuous chaotic motion or intermittency chaos motion. When ω = 29.0 rad/s, the system enters the chaotic motion through the saddle-node bifurcation, and in the range of ω > 28.8 rad/s, the vibration system exhibits strong nonlinear characteristics of periodic motion and chaotic motion frequently alternating. With the increase of the rotational speed, different bifurcation paths, such as saddle-node bifurcation, perioddoubling bifurcation, and inverse period-doubling bifurcation appear in the system. The corresponding vibration amplitude increases first, then decreases and then increases nonlinearly. In order to specify the main frequency components of the system under different excitation currents of the generator Ij, the main frequency components (where fn is the frequency conversion) in different speed ranges are listed in Table 3.   increases, the chaotic region becomes wider. In the range of ω ∈ [20.0,22.4] rad/s, except for the stable period-doubling motion in a small region, the system is of continuous chaotic motion or intermittency chaos motion. When ω = 29.0 rad/s, the system enters the chaotic motion through the saddle-node bifurcation, and in the range of ω > 28.8 rad/s, the vibration system exhibits strong nonlinear characteristics of periodic motion and chaotic motion frequently alternating. With the increase of the rotational speed, different bifurcation paths, such as saddle-node bifurcation, perioddoubling bifurcation, and inverse period-doubling bifurcation appear in the system. The corresponding vibration amplitude increases first, then decreases and then increases nonlinearly. In order to specify the main frequency components of the system under different excitation currents of the generator Ij, the main frequency components (where fn is the frequency conversion) in different speed ranges are listed in Table 3.         Through comparison, it can be seen that the increase of the excitation current has a certain inhibitory effect on the response amplitude of the system, but it has no inhibitory effect on the motion state, which makes the motion state of the system more complex, the chaotic motion widens relatively, and the jump discontinuity is enhanced. At higher rotational speeds, the number of high-order bifurcation increases with the increase of excitation current. Therefore, it is necessary to control the excitation current within a certain range in order to enable the vibration system to operate stably.  f n , 2f n , 3f n , 4f n , . . . , mf n /2 19.8 < ω ≤ 23.6 f n , 2f n , 3f n , 4f n , . . . , continuous frequency components 23.6 < ω ≤ 24.6 f n , 2f n , 3f n , 4f n , . . . , mf n /2 24.6 < ω ≤ 42.4 f n , 2f n , 3f n , 4f n , . . . , continuous frequency components 42.4 < ω ≤ 42.8 f n , 2f n , 3f n , 4f n , . . . . . . 42.8 < ω ≤ 54. 6 f n /2, f n , 3f n /2, 2f n , 5f n /2, 3f n , 7f n /2, 4f n , continuous frequency components 54.6 < ω ≤ 59. 6 f n , 2f n , 3f n , continuous frequency components 59.6 < ω ≤ 71 f n , 2f n , 3f n , mf n /3

Effect of Radical Stiffness on Dynamic Characteristics of RBRS at Variable Speed
The nonlinear radial stiffness (k r ) has an important influence on the stable operation and dynamic characteristics of RBRS. In order to study the dynamic response characteristics of the coupling system with nonlinear rubbing forces at different speeds, Figures 13-18 give the bifurcation diagram and 3-D frequency spectrum of RBRS with different radial stiffness (k r = 1 × 10 8 N/m, k r = 2.6 × 10 8 N/m, and k r = 4 × 10 8 N/m). All other parameters are fixed. the speed increasing, the chaotic motion region is widened, the intermittent chaotic phenomenon intensifies, and the continuum region increases in the three-dimensional spectrum. At the same time, the vibration system has obvious frequency demultiplication components; mfn/2, mfn/3, mfn/4, and mfn/6 (m: positive integer), and the frequency amplitude of fn/2, fn/3, fn/4 is significantly larger than other frequency amplitudes. In the high-speed region, the system finally transits to a stable period 4 motion through the saddle-node bifurcation, and the jumping and discontinuity of amplitude is enhanced. Table 4 shows the intervals in which the frequency components change (where fn is the frequency conversion).              It can be seen from Figures 13 and 14 that when the radial stiffness is small, the coupling system behaves as period 1 and a double-cycle motion, and only the bifurcation path of the period-doubling bifurcation and the inverse period-doubling bifurcation appears, and discontinuous jumps occur in the range of ω ∈ [39.6,52.2] rad/s. In the three-dimensional spectrum, when the rotational speed is small, the frequency components are relatively simple, mainly the frequency multiplication components (f n , 2f n , 3f n ,...). With the gradual increase of the rotational speed, besides the frequency multiplication component, the coupling system has unobvious frequency demultiplication components (f n /2, f n /3, f n /4,...) and the frequency demultiplication components alternately appear. The amplitude of the main frequency components increases first and then decreases nonlinearly (when ω = 58.0 rad/s, the c amplitude peak appears at f n ). When the radial stiffness increases to k r = 2.6 × 10 8 N/m, it can be found that the period doubling region widens, the period-doubling bifurcation and the inverse period-doubling bifurcation alternate frequently, and obvious frequency demultiplication components appear, and the transient quasi-periodic motion in high speed region. There is no chaos in the coupling system in the whole speed changing region. As the radial stiffness continues to increase to k r = 4 × 10 8 N/m, the coupling system exhibits complex nonlinear behavior and more unstable regions; different bifurcation paths, such as saddle-node bifurcation, period-doubling bifurcation, and inverse period-doubling bifurcation, and strong nonlinear phenomena, such as discontinuous amplitude jumps, appear. The chaotic region of the system is mainly concentrated in the three regions of ω ∈ [18.8,52.2] rad/s, ω ∈ [25.6,43.6] rad/s, and ω ∈ [54.6,67.4] rad/s, and with the speed increasing, the chaotic motion region is widened, the intermittent chaotic phenomenon intensifies, and the continuum region increases in the three-dimensional spectrum. At the same time, the vibration system has obvious frequency demultiplication components; mf n /2, mf n /3, mf n /4, and mf n /6 (m: positive integer), and the frequency amplitude of f n /2, f n /3, f n /4 is significantly larger than other frequency amplitudes. In the high-speed region, the system finally transits to a stable period 4 motion through the saddle-node bifurcation, and the jumping and discontinuity of amplitude is enhanced. Table 4 shows the intervals in which the frequency components change (where f n is the frequency conversion). f n , 2f n , 3f n , 4f n , 5f n , . . . . . . 39.6 < ω ≤ 52.2 f n , 2f n , 3f n ; mf n /2, mf n /3, and mf n /4 alternate 52.2 < ω ≤ 59.8 f n , 2f n , 3f n , 4f n , . . . . . . 59.8 < ω ≤ 71 f n , 2f n , 3f n ; mf n /2 and mf n /4 alternate k r = 2.6 × 10 8 f n , 2f n , 3f n , 4f n , . . . . . . 29.6 < ω ≤ 32. 6 f n , 2f n , 3f n , 4f n , . . . , mf n /2 32.6 < ω ≤ 35.8 f n , 2f n , 3f n , 4f n , . . . . . . 35.8 < ω ≤ 40.2 f n , 2f n , 3f n , 4f n , . . . ; mf n /18, mf n /23, and mf n /27 alternate 40.2 < ω ≤ 47.8 f n , 2f n , 3f n , 4f n , . . . , mf n /7 47.8 < ω ≤ 51. 6 f n , 2f n , 3f n , 4f n , . . . , mf n /2 51.6 < ω ≤ 60.4 f n , 2f n , 3f n , 4f n , . . . . . . 60.4 < ω ≤ 66.8 f n , 2f n , 3f n , 4f n , . . . , first mf n /2, then mf n /3 66.8 < ω ≤ 71 f n , 2f n , 3f n , mf n /7 k r = 4 × 10 8 f n , 2f n , 3f n , 4f n , . . . . . . 18.8 < ω ≤ 22.2 f n , 2f n , 3f n , 4f n , . . . , mf n /2, continuous frequency components 22.2 < ω ≤ 25.6 f n , 2f n , 3f n , 4f n , . . . . . .
In order to analyze more deeply, the related dynamic responses of the vibration system (k r = 4 × 10 8 N/m) at ω = 19.2 rad/s, ω = 40.8 rad/s and ω = 69.8 rad/s are displayed in Figures 19-21 through waveform, frequency, axis orbit, and Poincaré map. When the rotational speed is small, i.e., ω = 19.2 rad/s, the vibration amplitude of the system is small, and the corresponding main frequency domain components are frequency multiplication (f n , 2f n , 3f n , . . . ) and mf n /2 (m: positive integer) with small amplitude. The coupling system behaves as periodic-2 motion, and its attractors are two points on the Poincaré map. When ω = 59 rad/s, the amplitude of the system vibration response increases and presents a non-periodic change; in the frequency domain, the continuous frequency spectrum components appear, corresponding to the chaotic axis orbit; the attractors on the corresponding Poincaré map are irregular scattered points, and then the system motion is in chaotic state. When the rotational speed increases to ω = 69.8 rad/s, the rotor returns to the periodic motion again, and the vibration response has different frequency characteristics; f n /4 is the main frequency component and has obvious peak value. At the same time, four irregular closed curves are shown on the axis orbit map, and the attractors on the corresponding Poincare map are four isolated points. components appear, corresponding to the chaotic axis orbit; the attractors on the corresponding Poincaré map are irregular scattered points, and then the system motion is in chaotic state. When the rotational speed increases to ω = 69.8 rad/s, the rotor returns to the periodic motion again, and the vibration response has different frequency characteristics; fn/4 is the main frequency component and has obvious peak value. At the same time, four irregular closed curves are shown on the axis orbit map, and the attractors on the corresponding Poincare map are four isolated points.   Poincaré map are irregular scattered points, and then the system motion is in chaotic state. When the rotational speed increases to ω = 69.8 rad/s, the rotor returns to the periodic motion again, and the vibration response has different frequency characteristics; fn/4 is the main frequency component and has obvious peak value. At the same time, four irregular closed curves are shown on the axis orbit map, and the attractors on the corresponding Poincare map are four isolated points.   Poincaré map are irregular scattered points, and then the system motion is in chaotic state. When the rotational speed increases to ω = 69.8 rad/s, the rotor returns to the periodic motion again, and the vibration response has different frequency characteristics; fn/4 is the main frequency component and has obvious peak value. At the same time, four irregular closed curves are shown on the axis orbit map, and the attractors on the corresponding Poincare map are four isolated points.

Influence of Friction Coefficient on Dynamic Characteristics of RBRS under Variable Speed
Friction coefficient µ has an important influence on the dynamic behavior of RBRS, and its size directly affects the dynamic characteristics of the coupling system. Figures 22-27 show the dynamic response (bifurcation diagram and three-dimensional spectrum diagram) of the system under different working conditions when the rotational speed is in the range of ω ∈ [7.0,71.0] rad/s and the friction coefficient is µ = 0.05, µ = 0.3, and µ = 0.78, respectively. By comparison it can be seen that when the friction coefficient is small, that is, µ = 0.05, the RBRS in the interval of ω ∈ [7.0,60.2] rad/s mainly exhibits the period-1 motion, and the response amplitude is relatively stable (when ω = 51.2 rad/s, the coupling system behaves as chaotic motion; when ω = 52.6 rad/s, the coupling system behaves as a period-3 motion), and there are many obvious discontinuous amplitude jumps. With the increase of the rotational speed, the system transits from periodic motion to chaotic motion through saddle-node bifurcation at ω = 60.4 rad/s, and exhibits multiple motion states including periodic-1, periodic-2, periodic-n, Quasi-periodic, and chaotic motion. The coupling system appears as intermittent chaos, with the saddle-node bifurcation, the period-doubling bifurcation, and the inverse period-doubling bifurcation appearing alternately. In addition to the frequency multiplication component in the three-dimensional spectrum, there are also frequency demultiplication components (f n /2, 3 f n /25, f n /7,...) and continuous spectrum components; when µ = 0.3, the dynamic characteristics of the coupling system become more complex and the chaotic movement is strengthened. At lower rotational speeds, the jump of the coupling system is weakened, and the motion state returns to the stable periodic motion at 21.0 < ω < 27.8 rad/s again from the periodic motion in the range of ω ∈ [7.0,20.0] rad/s through the chaotic motion in the narrow interval ω ∈ [20.2,21.0] rad/s. When ω > 27.8 rad/s, the system enters the chaotic region through the saddle-node bifurcation, and the periodic motion, the quasi-periodic motion, and the chaotic motion appear alternately until ω = 61.4 rad/s. When 61.4 rad/s <ω < 63.2 rad/s, RBRS is in stable period-4 motion and period-8 motion and is accompanied by the occurrence of inverse period-doubling bifurcation; the three-dimensional spectrum shows distinct frequency demultiplication components of f n /4 and f n /8, and RBRS re-enters into the chaotic region at ω = 63.2 rad/s through the saddle-node bifurcation; when the friction coefficient is increased to µ = 0.78, while ω < 28.0 rad/s, the vibration system performs chaotic motion and periodic multiple motion alternatively only in the narrow interval ω ∈ [21.8,24.8] rad/s, and in other cases perform stable periodic motion. With the increase of the rotational speed, the system enters the unsteady region through the saddle-node bifurcation at ω = 28.2 rad/s. The system alternates frequently between the period-doubling motion, the counter-periodic motion, and the chaotic motion, and the jumping discontinuity of the response amplitude is intensified. Table 5 shows the intervals in which the frequency components change (where f n is the frequency conversion).
behaves as a period-3 motion), and there are many obvious discontinuous amplitude jumps. With the increase of the rotational speed, the system transits from periodic motion to chaotic motion through saddle-node bifurcation at ω = 60.4 rad/s, and exhibits multiple motion states including periodic-1, periodic-2, periodic-n, Quasi-periodic, and chaotic motion. The coupling system appears as intermittent chaos, with the saddle-node bifurcation, the period-doubling bifurcation, and the inverse period-doubling bifurcation appearing alternately. In addition to the frequency multiplication component in the three-dimensional spectrum, there are also frequency demultiplication components (fn/2, 3 fn/25, fn/7,...) and continuous spectrum components; when μ = 0.3, the dynamic characteristics of the coupling system become more complex and the chaotic movement is strengthened. At lower rotational speeds, the jump of the coupling system is weakened, and the motion state returns to the stable periodic motion at 21.0 < ω < 27.8 rad/s again from the periodic motion in the range of ω ∈ [7.0,20.0] rad/s through the chaotic motion in the narrow interval ω ∈ [20.2,21.0] rad/s. When ω > 27.8 rad/s, the system enters the chaotic region through the saddle-node bifurcation, and the periodic motion, the quasi-periodic motion, and the chaotic motion appear alternately until ω = 61.4 rad/s. When 61.4 rad/s <ω < 63.2 rad/s, RBRS is in stable period-4 motion and period-8 motion and is accompanied by the occurrence of inverse period-doubling bifurcation; the three-dimensional spectrum shows distinct frequency demultiplication components of fn/4 and fn/8, and RBRS re-enters into the chaotic region at ω = 63.2 rad/s through the saddle-node bifurcation; when the friction coefficient is increased to μ = 0.78, while ω < 28.0 rad/s, the vibration system performs chaotic motion and periodic multiple motion alternatively only in the narrow interval ω ∈ [21.8,24.8] rad/s, and in other cases perform stable periodic motion. With the increase of the rotational speed, the system enters the unsteady region through the saddle-node bifurcation at ω = 28.2 rad/s. The system alternates frequently between the period-doubling motion, the counter-periodic motion, and the chaotic motion, and the jumping discontinuity of the response amplitude is intensified. Table 5 shows the intervals in which the frequency components change (where fn is the frequency conversion).                In order to better display the influence of friction coefficient µ on the dynamic characteristics of the system, Figures 28-30 show that the time domain map, frequency domain map, axis orbit map, and Poincaré map of the coupling system with µ = 0.78, when ω = 29.8 rad/s, ω = 43.2 rad/s and ω = 52.2 rad/s, respectively. It can be seen from the Figures 28-30 that the vibration response of the three cases are nonlinear; in the frequency domain, f n , 2f n , and 3f n are the main frequency components, and there are frequency multiplication and frequency demultiplication components whose amplitude is not obvious. In addition, the axis orbit map changes from chaotic closed curves to multiple non-coincident closed curves, and then to five closed ones; and the attractors on the corresponding Poincaré map show irregular scatters, isolated multiple points, and isolated five points, respectively. Then the system completes the alternative evolution of chaotic motion, periodic-n motion, and period-5 motion.
and Poincaré map of the coupling system with μ = 0.78, when ω = 29.8 rad/s, ω = 43.2 rad/s and ω = 52.2 rad/s, respectively. It can be seen from the Figure 28-30that the vibration response of the three cases are nonlinear; in the frequency domain, fn, 2fn, and 3fn are the main frequency components, and there are frequency multiplication and frequency demultiplication components whose amplitude is not obvious. In addition, the axis orbit map changes from chaotic closed curves to multiple noncoincident closed curves, and then to five closed ones; and the attractors on the corresponding Poincaré map show irregular scatters, isolated multiple points, and isolated five points, respectively. Then the system completes the alternative evolution of chaotic motion, periodic-n motion, and period-5 motion.   52.2 rad/s, respectively. It can be seen from the Figure 28-30that the vibration response of the three cases are nonlinear; in the frequency domain, fn, 2fn, and 3fn are the main frequency components, and there are frequency multiplication and frequency demultiplication components whose amplitude is not obvious. In addition, the axis orbit map changes from chaotic closed curves to multiple noncoincident closed curves, and then to five closed ones; and the attractors on the corresponding Poincaré map show irregular scatters, isolated multiple points, and isolated five points, respectively. Then the system completes the alternative evolution of chaotic motion, periodic-n motion, and period-5 motion.   Based on the above numerical simulation, it can be seen that the friction coefficient μ is an important parameter affecting the dynamic characteristics of the coupling system (RBRS). With the increase of the friction coefficient, the motion state of the system becomes more complex, the chaotic region increases first and then decreases, the complex motion and periodic motion alternate frequently, and the response amplitude increases gradually.

Conclusions
In this paper, a model considering multiple nonlinear factors is proposed. The complex nonlinear dynamic behavior of the coupling system is studied. The influence of the excitation current, radial stiffness, and friction coefficient on the dynamic characteristics of RBRS is analyzed. The following conclusions can be drawn from the study.
(1) For the RBRS coupling system, with the change of parameters, there are many kinds of motion states, including periodic-n, quasi-periodic, and chaotic motion. At the same time, the coupling system has intermittent chaos many times, and Hopf, saddle-node, period-doubling, and inverse period-doubling bifurcation appear alternately. The dynamic properties of the coupling system caused by different nonlinear factors are interactional. (2) When the excitation current of the generator Ij is small, the dynamic behavior of the system is relatively simple. As the rotational speed increases, the system transits from short-term periodic Based on the above numerical simulation, it can be seen that the friction coefficient µ is an important parameter affecting the dynamic characteristics of the coupling system (RBRS). With the increase of the friction coefficient, the motion state of the system becomes more complex, the chaotic region increases first and then decreases, the complex motion and periodic motion alternate frequently, and the response amplitude increases gradually.

Conclusions
In this paper, a model considering multiple nonlinear factors is proposed. The complex nonlinear dynamic behavior of the coupling system is studied. The influence of the excitation current, radial stiffness, and friction coefficient on the dynamic characteristics of RBRS is analyzed. The following conclusions can be drawn from the study.
(1) For the RBRS coupling system, with the change of parameters, there are many kinds of motion states, including periodic-n, quasi-periodic, and chaotic motion. At the same time, the coupling system has intermittent chaos many times, and Hopf, saddle-node, period-doubling, and inverse period-doubling bifurcation appear alternately. The dynamic properties of the coupling system caused by different nonlinear factors are interactional. (2) When the excitation current of the generator I j is small, the dynamic behavior of the system is relatively simple. As the rotational speed increases, the system transits from short-term periodic motion, chaotic motion, and quasi-periodic motion to periodic motion, successively undergoing saddle-node bifurcation, period-doubling bifurcation, inverse period-doubling bifurcation, and Hopf bifurcation. The increase of the excitation current I j can inhibit the response amplitude of the system to some extent, but it cannot inhibit the motion state of the system, which makes the motion state of the system more complex, the chaotic motion wider, and the jump discontinuity enhanced. At higher speeds, the number of high-order bifurcation increases as the excitation current I j increases. (3) Due to the influence of the radial stiffness k r and the friction coefficient, the coupled system shows more complex nonlinear behavior. With the increase of radial stiffness k r , the motion complexity of the coupling system increases, the chaotic region increases, the response amplitude increases, and the vibration intensity increases. At the same time, in the high-speed region, the periodic motion region becomes narrow and tends to disappear. With the increase of the friction coefficient, the chaotic region increases first and then decreases, the complex motion and periodic motion alternate frequently, and the response amplitude increases gradually.

Conflicts of Interest:
The authors declare no conflict of interest.