Coupled Electromagnetic-Dynamic Modeling and Bearing Fault Characteristics of Induction Motors considering Unbalanced Magnetic Pull

Induction motors are complex energy conversion systems across the domains of dynamics, electricity, and magnetism. Most existing models mainly consider unidirectional coupling, such as the effect of dynamics on electromagnetic properties, or the effect of unbalanced magnetic pull on dynamics, while in practice it should be a bidirectional coupling effect. The bidirectionally coupled electromagnetic-dynamics model is beneficial to the analysis of induction motor fault mechanisms and characteristics. This paper proposes a coupled electromagnetic-dynamic modeling method that introduces unbalanced magnetic pull. By using the rotor velocity, air gap length, and unbalanced magnetic pull as the coupling parameters, the coupled simulation of the dynamic and electromagnetic models can be effectively realized. Simulation results for bearing faults show that the introduction of magnetic pull induces a more complex dynamic behavior of the rotor, which in turn leads to modulation in the vibration spectrum. The fault characteristics can be found in the frequency domain of the vibration and current signals. Through the comparison between simulation and experimental results, the effectiveness of the coupled modeling approach and the frequency domain characteristics caused by the unbalanced magnetic pull are verified. The proposed model can help to obtain a variety of information that is difficult to measure in reality and can also serve as a technical basis for further research on nonlinear characteristics and chaos in induction motors.


Introduction
Induction motors are the most widely used electromechanical energy conversion devices in the industrial production chain, which benefits from the significant advantages of their simple construction in terms of easy maintenance and reliable operation. However, while most induction motors are inexpensive and easy to replace, their importance in the task chain may not match the cost of maintenance. In mission situations where unplanned interruptions are not tolerated, the cost of induction motor faults can be incalculable. According to the statistics of IEEE-IAS and EPRI, for a large number of small and mediumsized motor faults, more than 40% of induction motor faults are related to bearings [1][2][3][4]. Therefore, the research on mechanism and prevention of induction motor bearing faults has been widely concerned.
Vibration from motor bearings can be transmitted to motor end shields and enclosures; thus, it seems that many of the vibration analysis methods applied to common bearings can be easily transferred to motors. Based on vibration signals, many research results on bearing faults have been published over the years, including, but not limited to, modeling and simulation, signal processing, feature extraction, life prediction, and intelligent diagnosis. Taking some of the recently published studies as examples, Zhao et al. [5] proposed a performance degradation prediction method based on high-order differential mathematical morphology gradient spectrum entropy, phase space reconstruction, and extreme learning machines to predict the performance degradation trend of rolling bearings. Zhu et al. [6] used a hidden Markov model to locate the fault occurrence time and solved the distribution discrepancy problem by a transfer learning method based on multiple layer perceptron. Cui et al. [7] proposed a coupled multistable stochastic resonance method with two first-order multistable stochastic resonance systems, which has higher output signal-to-noise ratio than multistable stochastic resonance. Dong et al. [8] used a bearing dynamics model to generate a large amount of simulated data for the small sample problem, and then implemented transfer learning between the diagnostic knowledge in the simulated data to the real scenario based on convolutional neural network and parameter transfer strategies. Considering the compensation balance excitation caused by the rotor mass eccentricity, Liu et al. [9] developed a dynamic model for the bearing rotor system of a high-speed train under variable speed conditions. The use of the angle iteration method in this case solves the problem that the space position of the rollers cannot be determined during bearing rotation. With the improvement of computer computing power and the rise of artificial intelligence algorithms, a large number of intelligent diagnostics works with excellent performance in public data sets have been published in recent years. However, how to explain the artificial intelligence model at the physical level and how to obtain enough fault sample data are still the key issues for the intelligent diagnosis method. By automatically generating data, the modeling approach that simulates machine faults makes it possible to solve the latter. Analytical models constructed according to the laws of physics can simulate faults of different sizes and locations and generate enough representative signals to train intelligent diagnostic algorithms. In addition, for complex cases occurring in the machine such as nonlinearity and coupling, an analytical model helps the understanding and analysis at the physical level [10].
Unlike general mechanical structures, induction motors are complex energy conversion systems with mechanism-electricity-magnetism coupling, which requires more detailed consideration in the research on it. However, as noted above, most studies focused on motor bearing faults have adopted analytical approaches similar to those of nonmotor bearing studies. Although they seem to work from the results, many simple transfer schemes ignore the electromagnetic force between the stator and rotor, i.e., the unbalanced magnetic pull (UMP) on the rotor. There have been studies and discussions about UMP since the beginning of the last century [11]; one of the reasons is to reduce the air gap in the design for improving the power factor [12]. The small air gap in rotating motors is susceptible to small variations in the dimensions of stator, rotor, bearings, etc. A nonuniform air gap will twist the flux density distribution, thus producing a UMP that tends to drag the rotor further away from its center position, which is difficult to ignore during motor operation. Therefore, UMP has been a long-standing concern in the electrical field, from early discussions of computational methods to current research on the dynamic behavior of motor rotors. Taking some of the research in the last three decades as an example, Dorrell and Smith [13,14] proposed a general analytical model for calculating the UMP in threephase squirrel-cage induction motors with static eccentricity and observed a favorable correlation between the measured and calculated results. In another study, Dorrell [15] considered the axial variation of the eccentricity in the calculation for the UMP. Guo [16] et al. derived an analytical expression of the UMP caused by eccentricity in a three-phase motor at no load and investigated the nonlinear vibration response excited by the UMP in the Jeffcott rotor model. Using the same Jeffcott rotor model, Chen et al. [17] derived an analytical expression for the UMP in a permanent magnet synchronous motor and introduced it into the nonlinear vibration equation to investigate the variation of the system response under different design parameters. To investigate the rotor dynamics stability, Han et al. [18] proposed a method to calculate the radial and axial motor eccentric forces by magnetic equivalent circuit modeling. Zhang et al. [19] developed a rotor bearing system dynamics model considering gyroscopic effect, nonlinear bearing force, and UMP force, On the other hand, studies on the dynamics of rotor-bearing systems have also led to many results. Cheng [20] et al. investigated the nonlinear dynamic behaviors of a rotor-bearing-seal coupled system by using the nonlinear seal fluid dynamic force model and nonlinear oil film force. Li et al. [21] proposed a generalized dynamics model for a rolling bearing-rotor system by coupling an explicit finite element model with a dynamic bearing model. Li et al. proposed a generalized dynamics model for a rolling bearingrotor system by coupling an explicit finite element model with a dynamic bearing model. Yan et al. [22] developed a fractional order mathematical model of the rotor-bearing-seal system and investigated the dynamic characteristics with changes in rotational speed, mass eccentricity of rotor, sealing clearance, and sealing pressure drop at a specific seal fractional order. Cheng et al. [23] investigated the dynamic properties of a rotor-bearing system with different compliances, radial clearances, rotor-stator rubbing, raceway defects, and surface waviness in an analytical model. Dynamic studies focusing on rotor-bearing systems have shown that faults occurring in the bearings will cause significant changes in rotor dynamic characteristics, which will affect the distribution of the air gap between the stator and rotor in the motor. A nonuniform and varying air gap will not only introduce the effect of UMP but will also affect the inductance between the stator and rotor, which will also affect the UMP through the current. Therefore, we believe that the modeling and fault characteristics study for motor bearings should be based on a bidirectional coupling of dynamic and electromagnetic models. Currently, only a few researchers have conducted studies on electromagnetic-dynamic coupled modeling for induction motors. Fourati et al. [24] proposed an electromagnetic-dynamic model to describe the dynamical behavior of an elastically supported rotating shaft coupled to a squirrel-cage induction motor. Han et al. [25] proposed an electromagnetic-dynamic coupling model based on magnetic equivalent circuits, which can consider nonlinear air gap permeance, nonlinear iron material, and magnetic saturation. The main difference between the models developed by the above studies is the choice of the induction motor electromagnetic model, which may produce differences in computational speed and accuracy. However, the existing coupled model only considers the unidirectional coupling from the dynamic model to the electromagnetic model, which may bring some unexplained simulation biases. Therefore, in order to investigate the effect of whether to consider the UMP on the simulation results in the coupled model, we propose an electromagnetic-dynamic coupled modeling method considering the UMP. UMP is introduced into the coupled model, which is based on rotorbearing dynamic model and multiple coupled circuit (MCC) theory for the first time. The dynamic response of the coupled model is solved by the numerical iterative algorithm. This paper is organized as follows. Section 2 briefly describes the modeling process of the coupled model, containing the theoretical formulations of the rotor-bearing dynamic model and the MCC model, the calculation methods of inductance and UMP, and the overall logic of the model. Section 3 presents the simulation, experimental results, and analysis of the coupled model. Finally, conclusions are given in Section 4.

Rotor-Bearing Dynamic Model
In contrast to the signal model, the simulation results provided by the dynamic model are more objective in reflecting the real physical situation. In this work, the nonlinearities of the rotor-bearing dynamics model are considered by introducing the bearing clearance and the additional radial clearance due to spalling. As shown in Figure 1, the radial force supporting the rotor rotation is provided by the Hertzian contact force between the ball and the inner ring in the bearing. The UMP acts directly on the rotor.
According to the Hertzian contact theory, the total restoring force of the bearing in the X and Y directions is Finally, considering the effect of UMP, the differential equation for the vibration of the rotor-bearing system is formulated as (1) Han et al. [25] simulated the additional radial clearance generated when the ball passes through the defective position by the half-wave sine function: According to the Hertzian contact theory, the total restoring force of the bearing in the X and Y directions is λ j = 1 when δ j > 0, otherwise, λ j = 0. Finally, considering the effect of UMP, the differential equation for the vibration of the rotor-bearing system is formulated as m r ..

Electromagnetic Model Based on Multiple Coupled Circuit Theory
In industrial applications, motors are inevitably operated under asymmetric supply conditions, accompanied by saturation effect, eddy current loss, and friction loss. If the influence of the above factors can be effectively introduced into the process of motor model construction, a simulation model closer to the real motor could be obtained. However, the consideration of detailed factors always implies the growth of model parameters and complexity of the structure, which leads to small improvements in results and huge increases in costs. It is not compatible with the cost-efficiency ratio requirements of industrial applications. Furthermore, the objective in this study is to construct coupled models that can effectively reflect the fault transmission paths and generate usable fault signals, which can help to qualitatively analyze and explain the fault mechanisms. Therefore, several popular simplifying assumptions were used in the construction of the MCC model to help improve the computational efficiency. These assumptions are as follows.

1.
The motor is powered by a balanced three-phase voltage source; 2.
Hysteresis loss, eddy current loss, and friction loss are ignored; 3.
Rotor bars are insulated from each other.
As shown in Figure 2, the MCC theory treats each stator winding phase as a series circuit of resistance and inductance and treats the squirrel-cage rotor as a loop evenly distributed in space. The voltage equation for a three-phase squirrel-cage induction motor when star-connected is where The motor dynamic equation:

Calculation of Inductance and Unbalanced Magnetic Pull
The rotor radial displacement determines the air gap distribution between the stator and rotor, as shown in Figure 3. Strictly speaking, the air gap in real conditions is always nonuniform due to manufacturing errors and gravity, etc., which is aggravated by the additional radial displacements caused by bearing faults. The bearing faults cause the dynamic eccentricity of the rotor due to the variation of the restoring force, so as to obtain the time-varying air gap. The distribution of the air gap can be calculated in successive numerical iterations. According to the displacement of the rotor in the x and y directions at moment t, the rotor eccentricity, eccentricity angle, and relative eccentricity introduced by the bearing faults can be obtained: where The electromagnetic torque equation in matrix form: The motor dynamic equation:

Calculation of Inductance and Unbalanced Magnetic Pull
The rotor radial displacement determines the air gap distribution between the stator and rotor, as shown in Figure 3. Strictly speaking, the air gap in real conditions is always nonuniform due to manufacturing errors and gravity, etc., which is aggravated by the additional radial displacements caused by bearing faults. The bearing faults cause the dynamic eccentricity of the rotor due to the variation of the restoring force, so as to obtain the time-varying air gap. The distribution of the air gap can be calculated in successive numerical iterations. According to the displacement of the rotor in the x and y directions at moment t, the rotor eccentricity, eccentricity angle, and relative eccentricity introduced by the bearing faults can be obtained: Entropy 2022, 24, x FOR PEER REVIEW 7 of 24 Hence, the distribution of the air gap length at moment t is The air gap length distribution is associated with the inductance through the winding function approach (WFA), which in turn affects the stator current in the MCC model. The WFA allows one to calculate the time-varying inductance between the stator and rotor by the motor winding structure and the air gap length distribution. The assumptions when applying the winding function method are as follows.
1. The flux linkage passes through the air gap radially, that is, the axial flux linkage is ignored; 2. The magnetic permeability of magnetic materials is infinite; 3. There is a negligible slot effect.
By derivation based on Gauss's law and Ampere's law, the equation for the WFA to calculate the mutual inductance between any coil 1 c and coil 2 c can be expressed as where the winding function    Hence, the distribution of the air gap length at moment t is The air gap length distribution is associated with the inductance through the winding function approach (WFA), which in turn affects the stator current in the MCC model. The WFA allows one to calculate the time-varying inductance between the stator and rotor by the motor winding structure and the air gap length distribution. The assumptions when applying the winding function method are as follows.

1.
The flux linkage passes through the air gap radially, that is, the axial flux linkage is ignored; 2.
The magnetic permeability of magnetic materials is infinite; 3.
There is a negligible slot effect.
By derivation based on Gauss's law and Ampere's law, the equation for the WFA to calculate the mutual inductance between any coil c 1 and coil c 2 can be expressed as where the winding function N c 2 (θ s ) is a function of the turn function n c 1 (θ s ) and the air gap length g(θ s ): Taking into account the relative motion between the stator and rotor, the WFA equation can also be rewritten as a function of the rotor's mechanical angular position or time t: The reason of UMP is that the nonuniform air gap generated by the eccentricity makes the flux density unevenly distributed. The presence of UMP tends to pull the rotor in an eccentric direction, which can lead to more complex rotor dynamics. Therefore, many studies have focused on the analytical calculation of the UMP. Based on the expression for the UMP proposed by Guo et al. [16], Zhou et al. [26] derived the general expression for the UMP at different number of poles: where: , n > 0 . Figure 4 shows the synthesis of the stator magnetomotive force (MMF) and rotor magnetomotive force when the motor is running with a load. Due to the assumption that hysteresis loss and eddy loss are not considered in the MCC model, the core loss angle α loss = 0 • . The rotor leakage reactance causes the rotor MMF to lag behind the induced electromotive force by an impedance angle ψ r . When the rotor leakage reactance is small, ψ r ≈ 0, and thus, ϕ sr = π − arccos F r F s .  Figure 4 shows the synthesis of the stator magnetomotive force (MMF) and rotor magnetomotive force when the motor is running with a load. Due to the assumption that hysteresis loss and eddy loss are not considered in the MCC model, the core loss angle loss 0    . The rotor leakage reactance causes the rotor MMF to lag behind the induced electromotive force by an impedance angle r  . When the rotor leakage reactance is small, r 0   , and thus, The stator and rotor MMF amplitudes are calculated by: The object of this paper is a 3-phase 2-pole squirrel-cage induction motor, whose stator windings are distributed in a single-layer and full-pitch manner. For the stator, the number of phases phase 3 n  , the total number of series turns per phase for single-layer winding distribution , the total number of series turns per phase 1 2 N  , and the fundamental winding factor w 1 k  .

Electromagnetic-Dynamic Coupled Model Framework
As shown in Figure 5, the bidirectional coupling between the rotor-bearing system dynamic model and the induction motor MCC model is achieved by the variables in the center coloring area. The blue frame represents the dynamic model of the rotor-bearing The stator and rotor MMF amplitudes are calculated by: The object of this paper is a 3-phase 2-pole squirrel-cage induction motor, whose stator windings are distributed in a single-layer and full-pitch manner. For the stator, the number of phases n phase = 3, the total number of series turns per phase for single-layer winding distribution N = p 2 Q s Z s , and the fundamental winding factor k w = k p k d = k p sin(Q s α/2) Q s sin(α/2) , where k p = 1 and α = p 2 · 2π Q s . For the rotor, the number of phases n phase = Q r p/2 , the total number of series turns per phase N = 1 2 , and the fundamental winding factor k w = 1.

Electromagnetic-Dynamic Coupled Model Framework
As shown in Figure 5, the bidirectional coupling between the rotor-bearing system dynamic model and the induction motor MCC model is achieved by the variables in the center coloring area. The blue frame represents the dynamic model of the rotor-bearing system, and the red frame represents the MCC model of the induction motor. In the dynamic model, the contact deformation equation includes a preset bearing fault state. Based on the Hertzian contact theory, the restoring force of the bearing in the X and Y directions can be calculated from the contact deformation and the ball position obtained by integrating the angular velocity. Since the UMP and restoring forces act on the rotor simultaneously, both are introduced into the vibration differential equation. The air gap length distribution is related to the stator-rotor inductance by the WFA, and it can be obtained by solving the vibration differential equation of the rotor-bearing system. The inductance between the stator and rotor directly affects the state of the motor MCC model and is reflected in the system's state parameters such as current and rotor angular velocity. The rotor angular velocity is fed back to the calculation of the bearing restoring force through the Hertzian contact theory, which is one of the paths by which the electromagnetic model affects the dynamic model. In addition, regarding the UMP calculated by current, relative eccentricity will be fed back to the vibration differential equation, potentially causing more complex rotor dynamic behavior.
on the Hertzian contact theory, the restoring force of the bearing in the X and Y directions can be calculated from the contact deformation and the ball position obtained by integrating the angular velocity. Since the UMP and restoring forces act on the rotor simultaneously, both are introduced into the vibration differential equation. The air gap length distribution is related to the stator-rotor inductance by the WFA, and it can be obtained by solving the vibration differential equation of the rotor-bearing system. The inductance between the stator and rotor directly affects the state of the motor MCC model and is reflected in the system's state parameters such as current and rotor angular velocity. The rotor angular velocity is fed back to the calculation of the bearing restoring force through the Hertzian contact theory, which is one of the paths by which the electromagnetic model affects the dynamic model. In addition, regarding the UMP calculated by current, relative eccentricity will be fed back to the vibration differential equation, potentially causing more complex rotor dynamic behavior.

Model Parameters and Experimental Conditions
The parameters in the coupled model are set according to the 3-phase 2-pole induction motor used in the experiments, as shown in Tables 1 and 2. Due to the limitations of measurement means and methods, a few mechanical and electrical parameters that were difficult to determine were estimated with reference to parameters from other works of the literature. The faults in the inner and outer races of the bearings are spalling faults with the same width. In order to highlight the bearing fault characteristics, the spalling position of the outer race fault is set directly below the inside of the outer race. Because of the complex rotational state of the ball in the raceway, this paper will not discuss the ball spalling fault in the bearing for the time being. The experimental bench and bearings used to simulate the faults are shown in Figure 6. The vibration signal on the motor end shield is measured by an IMI 608A11 accelerometer. The inverter provides a three-phase voltage at a specific frequency to the induction motor. The photoelectric sensor collects the pulse signal reflecting the rotational velocity by detecting the reflective strip located on the rotating shaft. The brake is connected to the induction motor via a rotor shaft and coupling to simulate the load at work. The bearing outer race spalling fault is simulated by cutting grooves on the inside of the outer race with a spalling width of about 3 mm, as shown in the Figure 6. Due to the low number of motor poles, this experiment was conducted at 20 Hz supply frequency for safety reasons.

Vibration Characteristics Analysis
This is performed for the purpose of studying the effect of UMP on the rotor dynamic behavior. The simulation time is set to 8 s, and the coupled model is used to simulate the normal, outer race fault, and inner race fault states of the induction motor in two cases: without magnetic pull and with magnetic pull. The rotor vibration time domain in Figure  7 shows the effect of whether or not UMP is introduced. The simulation results show that the healthy motor rotor has stationary vibration without UMP. In the presence of UMP, the vibration amplitude of the healthy motor rotor increases, and modulation occurs. The spalling located on the outer race of the bearing makes the rotor vibration waveform have regular shocks, while the presence of UMP seems to make the shock amplitude vary irregularly. In case of inner race faults, there is a significant amplitude modulation in the vibration waveform without UMP. With the introduction of the UMP, the amplitude in the originally low amplitude range (formed by amplitude modulation) increases significantly. In general, the effect of UMP on rotor vibration in the time domain is mainly shown as an increase in amplitude and jitter.

Vibration Characteristics Analysis
This is performed for the purpose of studying the effect of UMP on the rotor dynamic behavior. The simulation time is set to 8 s, and the coupled model is used to simulate the normal, outer race fault, and inner race fault states of the induction motor in two cases: without magnetic pull and with magnetic pull. The rotor vibration time domain in Figure 7 shows the effect of whether or not UMP is introduced. The simulation results show that the healthy motor rotor has stationary vibration without UMP. In the presence of UMP, the vibration amplitude of the healthy motor rotor increases, and modulation occurs. The spalling located on the outer race of the bearing makes the rotor vibration waveform have regular shocks, while the presence of UMP seems to make the shock amplitude vary irregularly. In case of inner race faults, there is a significant amplitude modulation in the vibration waveform without UMP. With the introduction of the UMP, the amplitude in the originally low amplitude range (formed by amplitude modulation) increases significantly. In general, the effect of UMP on rotor vibration in the time domain is mainly shown as an increase in amplitude and jitter.

Vibration Characteristics Analysis
This is performed for the purpose of studying the effect of UMP on the rotor dynamic behavior. The simulation time is set to 8 s, and the coupled model is used to simulate the normal, outer race fault, and inner race fault states of the induction motor in two cases: without magnetic pull and with magnetic pull. The rotor vibration time domain in Figure  7 shows the effect of whether or not UMP is introduced. The simulation results show that the healthy motor rotor has stationary vibration without UMP. In the presence of UMP, the vibration amplitude of the healthy motor rotor increases, and modulation occurs. The spalling located on the outer race of the bearing makes the rotor vibration waveform have regular shocks, while the presence of UMP seems to make the shock amplitude vary irregularly. In case of inner race faults, there is a significant amplitude modulation in the vibration waveform without UMP. With the introduction of the UMP, the amplitude in the originally low amplitude range (formed by amplitude modulation) increases significantly. In general, the effect of UMP on rotor vibration in the time domain is mainly shown as an increase in amplitude and jitter.   mainly in the horizontal direction, which may be due to bearing radial clearance. As shown in Figure 11, the radial clearance allows only a portion of the balls to be in contact with both the inner and outer races, so not all of the balls provide restoring force. The component of the combined restoring force in the horizontal direction depends on the rotor position and fluctuates with the rotor position. At the same time, the vertical vibration of the rotor caused by the restoring force fluctuation is limited by the presence of gravity. As shown in Figure 8, in the absence of UMP, the rotor axis has a clear trajectory despite the nonlinear factor caused by the radial clearance. When there is UMP, the range of the axis trajectory is significantly expanded, and the shape is close to the irregular motion under the condition of no UMP.  10 show the rotor axis trajectories for 7~8 s in normal, outer race fault, and inner race fault states, respectively. Under normal conditions, rotor displacement exists mainly in the horizontal direction, which may be due to bearing radial clearance. As shown in Figure 11, the radial clearance allows only a portion of the balls to be in contact with both the inner and outer races, so not all of the balls provide restoring force. The component of the combined restoring force in the horizontal direction depends on the rotor position and fluctuates with the rotor position. At the same time, the vertical vibration of the rotor caused by the restoring force fluctuation is limited by the presence of gravity. As shown in Figure 8, in the absence of UMP, the rotor axis has a clear trajectory despite the nonlinear factor caused by the radial clearance. When there is UMP, the range of the axis trajectory is significantly expanded, and the shape is close to the irregular motion under the condition of no UMP.      In the axial trajectory for the bearing outer race fault condition shown in Figure 9, the defect located at the bottom of the outer race causes the axial displacement in the vertical direction. The density of the axial trajectory reflects the difference in speed of the ball as it moves in and out of the defect position during counterclockwise rotation. The analysis of the process from position a to position d in Figure 9 is shown in Figure 12. As the deformation of the ball entering the defective position gradually decreases, the bearing restoring force supporting the rotor in the vertical direction gradually decreases. With gravity, the rotor accelerates from position a to position b. Position b to c is the process of the ball leaving the defective position. In this process, the ball deformation increases and the vertical restoring force increases. The resultant forces of the vertical restoring force and gravity support the rotor to accelerate upward after preventing it from continuing to fall. Due to the process of counteracting the falling momentum of the rotor, the rotor rises at a lower speed when the balls start to leave the defect position. In addition, before the ball enters the defect position, the resultant force of the bearing restoring force has a rightward horizontal component, and the rotor has a rightward horizontal velocity component. As the ball rotates counterclockwise away from the defect, although the rotor will rise to the vertical position before entering the defect, the speed of the rotor horizontally to the right at this time will tend to compress the balls leaving the defect. As a result, a greater restoring force of the ball leaving the defective position acts on the rotor and produces a behavior similar to "ejection", forming an approximately parabolic axial trajectory from position In the axial trajectory for the bearing outer race fault condition shown in Figure 9, the defect located at the bottom of the outer race causes the axial displacement in the vertical direction. The density of the axial trajectory reflects the difference in speed of the ball as it moves in and out of the defect position during counterclockwise rotation. The analysis of the process from position a to position d in Figure 9 is shown in Figure 12. As the deformation of the ball entering the defective position gradually decreases, the bearing restoring force supporting the rotor in the vertical direction gradually decreases. With gravity, the rotor accelerates from position a to position b. Position b to c is the process of the ball leaving the defective position. In this process, the ball deformation increases and the vertical restoring force increases. The resultant forces of the vertical restoring force and gravity support the rotor to accelerate upward after preventing it from continuing to fall. Due to the process of counteracting the falling momentum of the rotor, the rotor rises at a lower speed when the balls start to leave the defect position. In addition, before the ball enters the defect position, the resultant force of the bearing restoring force has a rightward horizontal component, and the rotor has a rightward horizontal velocity component. As the ball rotates counterclockwise away from the defect, although the rotor will rise to the vertical position before entering the defect, the speed of the rotor horizontally to the right at this time will tend to compress the balls leaving the defect. As a result, a greater restoring force of the ball leaving the defective position acts on the rotor and produces a behavior similar to "ejection", forming an approximately parabolic axial trajectory from position c to d in Figure 9. The range of axial trajectories is extended by adding UMP, which may be due to the fact that UMP always tends to pull the rotor away from the equilibrium position.
vertical restoring force increases. The resultant forces of the vertical restoring force and gravity support the rotor to accelerate upward after preventing it from continuing to fall. Due to the process of counteracting the falling momentum of the rotor, the rotor rises at a lower speed when the balls start to leave the defect position. In addition, before the ball enters the defect position, the resultant force of the bearing restoring force has a rightward horizontal component, and the rotor has a rightward horizontal velocity component. As the ball rotates counterclockwise away from the defect, although the rotor will rise to the vertical position before entering the defect, the speed of the rotor horizontally to the right at this time will tend to compress the balls leaving the defect. As a result, a greater restoring force of the ball leaving the defective position acts on the rotor and produces a behavior similar to "ejection", forming an approximately parabolic axial trajectory from position c to d in Figure 9. The range of axial trajectories is extended by adding UMP, which may be due to the fact that UMP always tends to pull the rotor away from the equilibrium position. As shown in Figure 10, when the bearing has an inner race fault, the balls can move in and out of the defective position at almost any angle because the spalled inner race rotates with the rotor. Therefore, a dynamic behavior similar to that described above for the outer race failure case can occur at all angles. Depending on the defect position angle, As shown in Figure 10, when the bearing has an inner race fault, the balls can move in and out of the defective position at almost any angle because the spalled inner race rotates with the rotor. Therefore, a dynamic behavior similar to that described above for the outer race failure case can occur at all angles. Depending on the defect position angle, the specific behavior may vary in part, as reflected in the more haphazard parabolic-like trajectory in Figure 10. It is worth noting that the overall shape in the dense area of the axial trajectory is close to the trajectory in the normal case, regardless of whether UMP is introduced.
Summarizing the rotor motion trajectories in the three bearing states, the simulation results of introducing UMP often cause the expansion in displacement range and the increase in disorder. cos β = 4.9319 f rm . Figure 13 shows the simulated displacement spectrum of the rotor in the Y direction without UMP. The rotor displacement spectrum curve in the case of healthy bearings is smooth and has prominent amplitudes only at the variable compliance frequency ( f VC , equal to the outer race fault frequency f om ≈ 56.43Hz) and its multiples (2 f VC , 3 f VC , · · · ). In the case of outer race fault, the amplitude at the outer race fault frequency ( f om ) and its multiples (2 f om , 3 f om , · · · ) increase significantly, while the amplitude on the other frequency components is also higher than that of the healthy bearing case. The inner race fault spectrum is significantly modulated by the rotational speed. In the marked red area, it can be seen that besides the inner race fault frequency ( f im ) and its multiples (2 f im , · · · ), the spectrum has the same prominent amplitude at the rotational frequency ( f rm ) and its multiples (2 f rm , 3 f rm , · · · ). With the addition of UMP, the vibration spectrum of the coupled model showed a more significant change. As shown in the areas highlighted in red in Figure 14, the speed modulation with the presence of UMP appears in the spectrum for the healthy and outer race fault cases. The spectrum of the inner race fault is less affected.
significantly modulated by the rotational speed. In the marked red area, it can be seen that besides the inner race fault frequency ( im f ) and its multiples ( im 2 , f  ), the spectrum has the same prominent amplitude at the rotational frequency ( rm f ) and its multiples ( rm rm 2 ,3 , f f  ). With the addition of UMP, the vibration spectrum of the coupled model showed a more significant change. As shown in the areas highlighted in red in Figure 14, the speed modulation with the presence of UMP appears in the spectrum for the healthy and outer race fault cases. The spectrum of the inner race fault is less affected. Figure 13. Rotor displacement spectrum in Y direction without unbalanced magnetic pull. Figure 13. Rotor displacement spectrum in Y direction without unbalanced magnetic pull. Experiments were conducted in the experimental bench in Figure 6, and the acceleration signal of the end shield on the motor drive end was measured. Figure 15 shows the time domain signal of motor end shield acceleration with outer race faulty bearing and normal bearing. Obviously, the shocks caused by the bearing faults make the peak acceleration in the fault condition higher than normal. From the perspective of frequency domain, the bearing outer race fault frequency ( om f ) and its multiples ( om om 2 ,3 , f f  ) can be seen in the vibration spectrum of the end shield at the drive end shown in Figure 16. The peaks in the spectrum at the industrial frequency ( s f ) and its multiples ( s s 2 ,3 , f f  ) may be caused by poor grounding. It is worth noting that the spectral curves in the case of outer race fault have significant peaks at rm f , rm 2 f , rm 3 f , and rm 4 f , which may correspond to the speed modulation caused by the UMP in Figure 14. In addition, the average amplitude located in the red area is slightly higher than the surrounding area, which is also Experiments were conducted in the experimental bench in Figure 6, and the acceleration signal of the end shield on the motor drive end was measured. Figure 15 shows the time domain signal of motor end shield acceleration with outer race faulty bearing and normal bearing. Obviously, the shocks caused by the bearing faults make the peak acceleration in the fault condition higher than normal. From the perspective of frequency domain, the bearing outer race fault frequency ( f om ) and its multiples (2 f om , 3 f om , · · · ) can be seen in the vibration spectrum of the end shield at the drive end shown in Figure 16. The peaks in the spectrum at the industrial frequency ( f s ) and its multiples (2 f s , 3 f s , · · · ) may be caused by poor grounding. It is worth noting that the spectral curves in the case of outer race fault have significant peaks at f rm ,2 f rm ,3 f rm , and 4 f rm , which may correspond to the speed modulation caused by the UMP in Figure 14. In addition, the average amplitude located in the red area is slightly higher than the surrounding area, which is also similar to the spectral pattern in Figure 14.
ation signal of the end shield on the motor drive end was measured. Figure 15 shows the time domain signal of motor end shield acceleration with outer race faulty bearing and normal bearing. Obviously, the shocks caused by the bearing faults make the peak acceleration in the fault condition higher than normal. From the perspective of frequency domain, the bearing outer race fault frequency ( om f ) and its multiples ( om om 2 ,3 , f f  ) can be seen in the vibration spectrum of the end shield at the drive end shown in Figure 16. The peaks in the spectrum at the industrial frequency ( s f ) and its multiples ( s s 2 ,3 , f f  ) may be caused by poor grounding. It is worth noting that the spectral curves in the case of outer race fault have significant peaks at rm f , rm 2 f , rm 3 f , and rm 4 f , which may correspond to the speed modulation caused by the UMP in Figure 14. In addition, the average amplitude located in the red area is slightly higher than the surrounding area, which is also similar to the spectral pattern in Figure 14.  the air gap variation characteristics caused by the bearing faults can be reflected in the electromagnetic model part of the coupled model. Figure 18 shows the experimentally measured stator current spectrum of the induction motor. The supply frequency (20 Hz) and its multiples are the main components of the spectrum. Affected by noise, only the characteristic frequencies 3 f om − f s and 5 f om − f s have obvious amplitude gain, and the characteristic frequencies f om − f s , f om + f s , 3 f om + f s , 4 f om + f s , and 5 f om + f s have weak amplitude gain. The comparison of simulation and experimental results verifies the correlation between the coupled model and the characteristic frequency formulation to a certain extent.

Conclusions
The simulation results show that the superposition of nonlinearities in the dynamical and electromagnetic models can induce complex dynamical behaviors reflected in the axial trajectories. These dynamical behaviors may be valuable for further analysis of the nonlinear and chaotic properties of the coupled model. Thus, although this paper is not a direct discussion of entropy and nonlinear mathematical problems, we still hope it can provide some reference for the extension of nonlinear theory in applications.
The construction of a model with effective coupling between dynamic model and electromagnetic model is beneficial for accurate analysis of the fault characteristics of the motor. Based on the velocity and air gap as the basic coupling paths, this study proposes an electromagnetic-dynamic coupled modeling method for induction motors considering UMP. In addition to the nonlinear factors such as bearing radial clearance and additional radial clearance, the rotor-bearing dynamics model based on the Hertzian contact theory introduces the effect of UMP caused by the radial displacement of the rotor. The change in air gap due to rotor vibration will affect the inductance calculation between the stator and rotor, thus serving as a path for the dynamics model to affect the electromagnetic model. The induction motor electromagnetic model based on the MCC theory feeds the rotor velocity directly into the dynamic model, while indirectly influencing the dynamic model through the UMP. In simulation results, the introduction of UMP makes the rotor show some special dynamical behaviors, and the amplitude and disorder of the vibration increases. In the vibration spectrum of the simulated signal, the UMP excites a speed modulation in normal and outer race fault conditions similar to that in the inner race fault case. The simulated spectrum curve located in the speed modulation range has a high amplitude, which is also reflected in the low frequency interval of the experimental signal spectrum. In the stator current simulation results of the coupled model, the frequency component associated with the fault can be found, which indicates the validity of the coupling path from the dynamic model to the electromagnetic model by air gap.
The coupled induction motor model integrating the rotor-bearing dynamic model and the MCC model has the ability to simulate a variety of mechanical and electrical faults, providing a more comprehensive approach for studying the dynamic behavior of motor rotors and fault characterization under different fault states. A correctly constructed coupled model is able to effectively simulate the motor state at different fault levels, which provides the possibility to track and predict the fault characteristics at each fault stage. Considering the industry's demand for confidence in the prediction data, the coupled model constructed entirely based on the mechanism has the prerequisites to provide reliable results.
Due to the experimental conditions, this paper has been experimentally verified only for normal and outer race fault. In addition, because of the lack of measured data on the small rotor displacements, the roughness of the model validation work must be acknowledged. In further work, we will improve the experimental conditions and try to simulate the complex rotational motion of ball fault so as to incorporate the experiment of inner race fault and the simulation of ball fault into this study. In order to further verify the rotor dynamics in the simulation results, the measurement of the radial displacement of the rotor is also forthcoming. Restoring force of the j-th bearing in the X direction F yj−restoring (j = 1, 2) Restoring force of the j-th bearing in the Y direction F x−magnetic Unbalanced magnetic pull of the rotor in the X direction F y−magnetic Unbalanced magnetic pull of the rotor in the Y direction g Gravitational acceleration D b Bearing pitch diameter n ball Number of bearing balls R e End ring resistance i ri (i = 1, 2, . . . , n b ) Current through the i-th rotor circuit λ ri (i = 1, 2, . . . , n b ) Flux linkage through i-th rotor circuit i e Current through the end ring λ e Flux linkage through the end ring L ss Stator self-inductance matrix L rr Rotor self-inductance matrix L sr Stator-rotor mutual inductance matrix L rs Rotor-stator mutual inductance matrix L sisj (i = 1, 2, 3; j = 1, 2, 3) Mutual inductance between stator phase i and stator phase j L sirj (i = a, b, c; j = 1, 2, . . . , n b ) Mutual inductance between stator phase i and rotor circuit j L rirj (i = 1, 2, . . . , n b , j = 1, 2, . . . , n b ) Mutual inductance between the i-th and j-th rotor circuits L b Rotor bar leakage inductance L e End ring leakage inductance