Dynamics Modeling and Analysis of Rolling Bearings Variable Stiffness System with Local Faults

: By analyzing the support of load-carrying rolling elements when the rolling elements fall into the fault position, the dynamics model of a rolling bearing variable stiffness system with local faults is proposed, considering the retention factor of the contact deformation. Then, this paper researches the change of effective contact stiffness, contact deformation, contact force


Introduction
With the continuous advancement of industrial technology, rotating machinery is developing in the direction of automation, integration, intelligence, high speed, and precision, which has led to an increase in the fault rate. As a key component widely used in rotating machinery, rolling bearing is prone to failure under harsh working conditions [1]. The unexpected fault of bearings can lead to a sudden collapse of a machine or system, and it may result in significant economic losses or even casualties [2]. Bearing fault diagnosis is mainly performed in two aspects; on the one hand, the fault information is extracted from the vibration response signal, which is mainly extracted and diagnosed by various learning algorithms. For example, the deep learning approach [3,4], digital twin-driven approach [5], optimized adaptive deep belief network [6], etc. [7,8]. With the increase in the use of smart machines, the detection and diagnosis of mechanical faults by these methods are increasing day by day. On the other hand, bearing dynamics models are mainly studied. Rolling bearing fault dynamics is used as a comprehensive and accurate method to predict the vibration characteristics of rolling bearings with various faults and to provide in-depth guidance for detection and diagnosis applications [9]. Fault diagnosis mechanisms and Machines 2023, 11, 609 2 of 22 methods for ball bearings are an ongoing research focus [10]. The dynamic models of shaft bearing systems have been developed for theoretical studies [9,11].
Cui et al. [12][13][14][15][16] developed a rolling bearing dynamics model by introducing the circumferential and radial dimensional parameters of two-dimensional faults into the displacement excitation function, estimated the magnitude of the faults by analyzing the time interval characteristics of the impact response, predicted the effect of local defects on ball bearing vibration, and also studied the vibration of deep groove ball bearings with single and multiple defects on the inner and outer race surfaces. Wu et al. [17][18][19][20] described the geometric displacement of the rolling body through a two-dimensional fault with a segmentation function and used this displacement as displacement excitation to establish a rolling bearing dynamics model, observed the relationship between vibration response and fault size, considered the coupling of rolling elements and segmentation effects, brought the acceleration response of the model more in line with the actual situation, and also considered the influence of the dent shoulder on the vibration, and derived a more close to real impulse characteristics caused by a practical dent. To build the rolling bearing dynamics model, Qui et al. [21,22] use the segmented function to introduce the three-dimensional geometric parameters of the fault into the displacement excitation. They analyze the variation of contact force of rolling elements under different contact types and the relationships between the fault size and the system vibration response. Zhang et al. [23,24] considered the dynamic lubrication conditions in an elastic fluid, described the geometric displacement of a rolling body as it passes through a fault using a segmentation function, used this displacement as a displacement excitation, developed a rolling bearing dynamics model, introduced the transient collision force excited by the strike of the rolling element on the trailing edge of the spall area, and analyzed the double pulse time interval, and through additional deflection and multi-impact theories, it was found that the location and the number of impulses due to varying compliance strongly depend on multiple factors, and mainly on the values of applied load and shaft rotational speed. Gao et al. [25,26] coupled displacement-excited rolling bearings in rail vehicles and rotor systems, respectively. Petersen et al. [27][28][29] established the system dynamics equations based on the displacement excitation function and analyzed the characteristics of the system stiffness and contact force change under different fault dimensions, i.e., the stiffness decreases in the loading direction and increases in the unloading direction, and also proposed a method for accounting for the finite rolling element size, which means that the time-frequency characteristics of the low-frequency event that occurs when a rolling element enters the defect entry and the multiple high-frequency events that occur when it exits the defect can be predicted more accurately. Sarabjeet et al. [30] used the explicit dynamics finite element software package LS-DYNA to build a rolling bearing fault model, and after an in-depth analysis of the numerically estimated dynamic contact forces between the rolling elements and the raceways of a bearing, it was found that the re-stressing of the rolling elements that occurs near the end of a raceway defect generates a burst of multiple short-duration force impulses, and the contact forces and accelerations generated on the exit of the rolling elements out of the defect compared to when they strike the defective surface are much higher.
The above-mentioned rolling bearing faulty dynamic models can be used in the analysis and diagnosis of system faults. In these dynamic models, the maximum displacement of the impact excitation is calculated by the fault size and the bearing geometric parameters. Additionally, the contact deformation and force are calculated by the displacement when the rolling elements fall into the fault position. Because of the rigid nature of the rolling elements, the deformation is generally a low order of magnitude. The maximum displacement is far greater than the contact deformation, so there is a large deviation.
In general, the radial load of the bearing is balanced by the contact force of the multiple load-bearing rolling elements. When a rolling element falls into the fault position and loses all or part of its load-carrying capacity, the radial load is redistributed between the load-bearing rolling elements. Additionally, the contact deformation and contact force of each load-bearing rolling element will be increased. When the rolling element loses its load capacity, that is, its contact stiffness, the total contact stiffness of the system will be changed suddenly. Then, the shock of the system is caused.
Based on Hertz's contact theory, the retention factor of contact deformation will be defined in this paper. During the rotation of different rolling elements with the cage, this paper analyzes the variation of the effective contact stiffness, contact force, and total effective stiffness of the system in the direction of radial load. A dynamic model of a single degree of freedom variable stiffness of a rolling bearing system with a single local damage fault is established, and simulation analysis and experimental verification are carried out.

Establish a Dynamic Model of the System with Fault
The deep groove ball bearing structure under radial load is shown in Figure 1a, the numbers in the figure are rolling element serial numbers, simplifying the elastic contact of a single rolling element with the inner and outer races into springs and dampers, as shown in Figure 1b. In Figure 1, O o and O i are the centers of the outer and inner races, β 0 is the initial position angle of the fault location, ψ is the azimuth at which the rolling elements are located, Q max is the maximum contact force in the load zone, f m is the cage rotation frequency, f s is the inner race rotation frequency, F r is an external radial load, the amplitude and direction are fixed, and δ r is the maximum relative displacement between the inner and outer races. Assume that the x-axis is in the radial load direction and the positive direction is the same as the radial load direction and the sensor is set on the outer race in the positive direction of the x-axis. In general, the radial load of the bearing is balanced by the contact force of the multiple load-bearing rolling elements. When a rolling element falls into the fault position and loses all or part of its load-carrying capacity, the radial load is redistributed between the load-bearing rolling elements. Additionally, the contact deformation and contact force of each load-bearing rolling element will be increased. When the rolling element loses its load capacity, that is, its contact stiffness, the total contact stiffness of the system will be changed suddenly. Then, the shock of the system is caused.
Based on HerR's contact theory, the retention factor of contact deformation will be defined in this paper. During the rotation of different rolling elements with the cage, this paper analyzes the variation of the effective contact stiffness, contact force, and total effective stiffness of the system in the direction of radial load. A dynamic model of a single degree of freedom variable stiffness of a rolling bearing system with a single local damage fault is established, and simulation analysis and experimental verification are carried out.

Establish a Dynamic Model of the System with Fault
The deep groove ball bearing structure under radial load is shown in Figure 1a, the numbers in the figure are rolling element serial numbers, simplifying the elastic contact of a single rolling element with the inner and outer races into springs and dampers, as shown in Figure 1b. In Figure 1, and are the centers of the outer and inner races, is the initial position angle of the fault location, is the azimuth at which the rolling elements are located, is the maximum contact force in the load zone, is the cage rotation frequency, is the inner race rotation frequency, is an external radial load, the amplitude and direction are fixed, and is the maximum relative displacement between the inner and outer races. Assume that the x-axis is in the radial load direction and the positive direction is the same as the radial load direction and the sensor is set on the outer race in the positive direction of the x-axis.  The situation when the rolling element falls into the fault position is shown in Figure 2. In Figure 2, h s is the maximum theoretical distance that the rolling element can fall when it falls into the fault position, r is the radius of the rolling element, r o is the radius of the outer raceway, r i is the radius of the inner raceway, b is the fault width, and h d is the fault depth.  In Figure 2, ℎ is the maximum theoretical distance that the rolling element can fall when it falls into the fault position, is the radius of the rolling element, is the radius of the outer raceway, is the radius of the inner raceway, is the fault width, and ℎ is the fault depth.
In the case of < , the maximum theoretical geometric distance of the collapse when the rolling element falls into the fault position is ℎ = ℎ ℎ ≥ ℎ Bottoming out # + # Inner race ℎ < ℎ # − # Outer race ℎ < ℎ Not bottoming out According to the deformation coordination conditions, the radial contact deformation at any azimuth is given by [31] where is the maximum value of contact deformation at = 0° and = is the load distribution range coefficient [15,23].
where C is the radial clearance and is the maximum relative radial displacement of the inner and outer races at = 0° [15,31]. In the case of b < r, the maximum theoretical geometric distance of the collapse when the rolling element falls into the fault position is In the formula, where subscripts i, o correspond to the inner and outer races, respectively. The central angle corresponding to the fault width is Outer race (4) According to the deformation coordination conditions, the radial contact deformation at any azimuth ψ is given by [31] where δ max is the maximum value of contact deformation at ψ = 0 • and ε is the load distribution range coefficient [15,23].
where P d is the radial clearance and δ r is the maximum relative radial displacement of the inner and outer races at ψ = 0 • [15,31].

of 22
Under the action of radial load F r , the angular range of the load zone, is [24,31] ψ L = ±arccos P d 2δ r (8) The relationship between the contact deformation δ ψ and the contact force Q(ψ) is [14][15][16]28] where K is the radial contact stiffness coefficient at a single rolling element [21,22].
where K ρi and K ρo are the radial contact stiffness coefficients of the rolling elements and the inner and outer raceways, please refer to the literature [31] for specific calculations. Therefore, the contact force Q(ψ) can be expressed as [15,23] For ball bearings subjected to a single radial load [31] where Z is the number of rolling elements and α is the contact angle. When the rolling element with an angle of ψ falls into the fault position, according to the geometric relationship between the amount of the contact deformation δ ψ and the fall distance h s , the retention factor of contact deformation is (1) At h s < δ ψ , then 0 < γ < 1, that is, the contact deformation at the rolling element is partially released, providing a partially effective contact load. (2) At δ ψ ≤ h s , then γ = 0, that is, the contact deformation at the rolling element is all released, and the effective contact load cannot be provided.
The change of retention factor γ with distance h s is shown in Figure 3.
Machines 2023, 11, x FOR PEER REVIEW 5 of 23 Under the action of radial load , the angular range of the load zone, is [24,31] The relationship between the contact deformation : and the contact force > ? is [14][15][16]28] > ? = G : H.J (9) where K is the radial contact stiffness coefficient at a single rolling element [21,22].
where G M and G M are the radial contact stiffness coefficients of the rolling elements and the inner and outer raceways, please refer to the literature [31] for specific calculations. Therefore, the contact force > ? can be expressed as [15,23] For ball bearings subjected to a single radial load [31] = R 4.37 VcosW where Z is the number of rolling elements and W is the contact angle. When the rolling element with an angle of falls into the fault position, according to the geometric relationship between the amount of the contact deformation : and the fall distance ℎ , the retention factor of contact deformation is (1) At ℎ < : , then 0 < Z < 1, that is, the contact deformation at the rolling element is partially released, providing a partially effective contact load. (2) At : ≤ ℎ , then Z = 0, that is, the contact deformation at the rolling element is all released, and the effective contact load cannot be provided.
The change of retention factor Z with distance ℎ is shown in Figure 3.   Figure 3 shows that with the increase in h s , the release of the deformation increases gradually, and the retention factor γ decreases gradually. When h s ≥ δ ψ , the amount of all deformation is released, the rolling elements cannot provide effective support, and the retention factor γ is reduced to 0.
From Equation (9), it can be seen that the contact load when a rolling element with an angle of ψ falls into the fault position is Let the No.1 rolling element be located in this position when ψ = 0 • , as shown in Figure 4a. The angle at which the i-th rolling element rotates with the cage is [20,21] where the cage's rotation angle is θ m = 2π f m t.
Machines 2023, 11, x FOR PEER REVIEW 6 of 23 Figure 3 shows that with the increase in ℎ , the release of the deformation increases gradually, and the retention factor Z decreases gradually. When ℎ ≥ : , the amount of all deformation is released, the rolling elements cannot provide effective support, and the retention factor Z is reduced to 0.
From Equation (9), it can be seen that the contact load when a rolling element with an angle of falls into the fault position is Let the No.1 rolling element be located in this position when = 0°, as shown in Figure 4a. The angle at which the i-th rolling element rotates with the cage is [20,21] where the cage's rotation angle is ^ = 2_ `.  where mod>•? is the remainder operation; parameter a determines whether the rolling element is located in the load zone. When the outer race fails, the contact stiffness at the i-th rolling element is When the inner race fails, the rotation angle of the fault point with the inner race is The contact stiffness at the i-th rolling element is To satisfy the static equilibrium relationship between the system and the radial load , must be equal to the sum of the components of loads of each rolling element [15].
where is the azimuth of the i-th rolling element. Cause where mod(·) is the remainder operation; parameter λ determines whether the rolling element is located in the load zone. When the outer race fails, the contact stiffness at the i-th rolling element is When the inner race fails, the rotation angle of the fault point with the inner race is [23] θ s = 2π f s t + β 0 The contact stiffness at the i-th rolling element is To satisfy the static equilibrium relationship between the system and the radial load F r , F r must be equal to the sum of the components of loads of each rolling element [15].
where ψ i is the azimuth of the i-th rolling element. The static equilibrium relationship of the system when there is no fault can be expressed as F r = k eqa δ 1.5 max (21) where k eqa is the total effective stiffness in the x-direction when there is no fault.
The effective contact stiffness of the i-th rolling element in the x-direction can be expressed as In particular, when the radial clearance P d is 0 and a rolling element is located directly below the radial load F r , there is just an odd number of rolling elements to support, as shown in Figure 4a, and the total effective stiffness of the system is When the j-th rolling element falls into the outer or inner race fault position, the static equilibrium relationship of the system is F r = k eqb δ max 1.5 (25) where δ max is the nominal maximum contact deformation in the load zone at ψ = 0 • ; k eqb is the total effective stiffness in the x-direction in case of fault.
The effective contact stiffness of the j-th rolling element in the x-direction can be expressed as In particular, when the radial clearance P d is 0, the j-th rolling element and point of the outer or inner race fault are directly below the radial load F r (ψ j = 0 • ), and the total effective stiffness of the system is Further, if δ j ≤ h s , then an even number of rolling elements are supported, as shown in Figure 4c, and the total effective stiffness of the system is Therefore, the total effective stiffness of the system is Take a deep groove ball bearing with nine rolling elements in Figure 4 as an example to illustrate the difference between Equations (24) and (29). Figure 4a shows that the normal bearing with exactly one rolling element located directly below the radial load F r and the number of load-carrying rolling elements comprises odd numbers. The rolling elements are rotated with the cage to the position of Figure 4b, and the number of load-carrying rolling elements is an even number; the azimuth of the No. 9 and No. 1 rolling elements on both sides of the radial load F r are +π/Z and −π/Z, respectively.
In Figure 4c, the No. 1 rolling element, located below the radial load F r , falls into the inner or outer race fault position and δ max < h s . When the contact deformation is fully released, the number of load-carrying rolling elements suddenly changes from an odd number to an even number, and the azimuth of the No.9 and No.2 rolling elements on both sides of the radial load F r are +2π/Z and −2π/Z, respectively. In this case, the lack of support for the No. 1 rolling element that provides the largest radial stiffness causes the total effective stiffness to decrease greatly so that the contact deformation of other load-carrying rolling elements increases to rebalance the radial load F r .
The system is simplified in the x-direction to a single-degree-of-freedom system with variable stiffness under radial load F r , as shown in Figure 5.
to illustrate the difference between Equations (24) and (29). Fig  normal bearing with exactly one rolling element located directly b and the number of load-carrying rolling elements comprises odd elements are rotated with the cage to the position of Figure 4b, an carrying rolling elements is an even number; the azimuth of the N elements on both sides of the radial load are +_/V and −_/V, In Figure 4c, the No. 1 rolling element, located below the radi inner or outer race fault position and v < ℎ . When the conta released, the number of load-carrying rolling elements suddenly number to an even number, and the azimuth of the No.9 and No both sides of the radial load are +2_/V and −2_/V, respective of support for the No. 1 rolling element that provides the largest ra total effective stiffness to decrease greatly so that the contact defo carrying rolling elements increases to rebalance the radial load .
The system is simplified in the x-direction to a single-degree-o variable stiffness under radial load , as shown in Figure 5. In Figure 5, i mn is the coefficient of the total effective stiffn indicates that the value is variable, ~m n is the damping coefficien mass of the outer or inner race.
In the case where there is no force on the bearing, the system established with the variable x as the origin when the inner and ou The differential equation of motion for the system is Combining the above formulas, the differential Equation (3 fourth-order Runge-Kuta method to obtain the vibration response

Bearing Parameters and Conditions
Taking the deep groove ball bearing of SKF6205 as an parameters are shown in Table 1. In Figure 5, k eq is the coefficient of the total effective stiffness, the slanted arrow indicates that the value is variable, c eq is the damping coefficient [27,32], and m is the mass of the outer or inner race.
In the case where there is no force on the bearing, the system coordinate system is established with the variable x as the origin when the inner and outer races are concentric.
The differential equation of motion for the system is Combining the above formulas, the differential Equation (31) is solved using the fourth-order Runge-Kuta method to obtain the vibration response of the system.

Bearing Parameters and Conditions
Taking the deep groove ball bearing of SKF6205 as an example, the bearing parameters are shown in Table 1. The initial condition of the rolling bearing is shown in Figure 4c. At this time, the inner race rotates counterclockwise, and the initial angle of the inner and outer race faults β 0 is 0 • . The No. 1 rolling element is located below the radial load F r , and δ max < h s , the radial clearance P d is 0, and the angle of the load zone ψ L is ±90 • .

Simulation of Total Effective Stiffness k eq
From the simulation of Equations (22) and (30), the simulation of the change of the total effective stiffness of the system with the angle of rotation of the cage θ m during faultfree, inner race, and outer race faults is obtained as shown in Figure 6. The numbers in Figure 6b,c indicate the serial number of the rolling elements through the fault location. The initial condition of the rolling bearing is shown in Figure 4c. At this time, the inner race rotates counterclockwise, and the initial angle of the inner and outer race faults is 0 . The No. 1 rolling element is located below the radial load , and v < ℎ , the radial clearance C is 0, and the angle of the load zone D is ±90°.

Simulation of Total Effective Stiffness i mn
From the simulation of Equations (22) and (30), the simulation of the change of the total effective stiffness of the system with the angle of rotation of the cage ^ during faultfree, inner race, and outer race faults is obtained as shown in Figure 6. The numbers in Figure 6b,c indicate the serial number of the rolling elements through the fault location. As shown in Figure 6a, when there is no fault, each rolling element rotates with the cage, so that the number of carrying rolling elements in the x-direction changes periodically from odd numbers (Figure 4a) to even numbers (Figure 4b), and then to odd numbers (Figure 4a). Additionally, the total effective stiffness of the system i mn fluctuates periodically within a small range. Figure 6b shows the change in the total effective stiffness i mn of the system when the outer race has fault. The outer race fault is located in the load zone and the relative radial load , and the position of the load zone is fixed. When each rolling element falls into the fault position in turn, the overall support situation of other rolling elements is the same. As shown in Figure 7, the total effective stiffness is reduced from i mno to i mnu , and the i mnu is equal in amplitude. The value of i mnu is related to the location angle of the fault , independent of the angle of rotation of the inner race. As shown in Figure 6a, when there is no fault, each rolling element rotates with the cage, so that the number of carrying rolling elements in the x-direction changes periodically from odd numbers (Figure 4a) to even numbers (Figure 4b), and then to odd numbers ( Figure 4a). Additionally, the total effective stiffness of the system k eq fluctuates periodically within a small range. Figure 6b shows the change in the total effective stiffness k eq of the system when the outer race has fault. The outer race fault is located in the load zone and the relative radial load F r , and the position of the load zone is fixed. When each rolling element falls into the fault position in turn, the overall support situation of other rolling elements is the same. As shown in Figure 7, the total effective stiffness is reduced from k eqa to k eqb , and the k eqb is equal in amplitude. The value of k eqb is related to the location angle of the fault β 0 , independent of the angle of rotation of the inner race.  Figure 6b shows the change in the total effective stiffness i mn of the system when the inner race has fault. Because the fault rotates with the inner race, its position varies periodically with respect to the radial load and the position of the load zone. When the rolling element and the fault point are in the load zone and meet, it causes the contact stiffness at the rolling element to change, resulting in the total effective stiffness decreasing from i mno to i mnu . The value of i mnu is related to the position of the rolling element when it meets the fault point in the load zone. The i mnu is variable in amplitude and it is modulated by the rotation of the inner race.

Simulation of Stiffness, Contact Force, and Contact Deformation of Each Rolling Element in Case of Outer Race Fault
From the simulation of Equations (23) and (27), it is obtained that the change of effective contact stiffness of each rolling element in the x-direction during the cage rotates one turn in case of the outer race fault as shown in Figure 7. Each rolling element calculates its angle of rotation with the cage starting from the initial azimuth. The key position of  Figure 6b shows the change in the total effective stiffness k eq of the system when the inner race has fault. Because the fault rotates with the inner race, its position varies periodically with respect to the radial load F r and the position of the load zone. When the rolling element and the fault point are in the load zone and meet, it causes the contact stiffness at the rolling element to change, resulting in the total effective stiffness decreasing from k eqa to k eqb . The value of k eqb is related to the position of the rolling element when it meets the fault point in the load zone. The k eqb is variable in amplitude and it is modulated by the rotation of the inner race.

Simulation of Stiffness, Contact Force, and Contact Deformation of Each Rolling Element in Case of Outer Race Fault
From the simulation of Equations (23) and (27), it is obtained that the change of effective contact stiffness of each rolling element in the x-direction during the cage rotates one turn in case of the outer race fault as shown in Figure 7. Each rolling element calculates its angle of rotation with the cage starting from the initial azimuth. The key position of each rolling element entering and exiting the load zone and falling into the fault position is shown in Table 2. Initial azimuth ψ 0i /  Figure 7a shows that the No. 1 rolling element falls into the fault position and is unable to provide support at this time. After the No. 1 rolling element is out of the fault with the rotation of the cage, its azimuth gradually increases, making the effective contact stiffness provided gradually decrease. The No. 1 rolling element rotates 90 • with the cage and then exits the load zone, and the effective contact stiffness drops to 0 at this time. However, the initial azimuth of the No. 9 rolling element is 320 • ; its azimuth gradually decreases with the rotation of the cage, making the effective contact stiffness provided gradually increase. The No. 9 rolling element falls into the fault position when it rotates 40 • with the cage, the effective contact stiffness drops to 0 at this time. After the No. 9 rolling element becomes out of the fault with the rotation of the cage, its azimuth gradually increases, making the effective contact stiffness provided gradually decrease. The No. 9 rolling element rotates 90 • with the cage and then exits the load zone; the effective contact stiffness drops to 0 at this time, as shown in Figure 7i. The case of other rolling elements is similar and will not be repeated.
Each rolling element is 40 • apart, and they fall into the fault position in turn when the cage rotates counterclockwise. Therefore, the entire rolling bearing in the load zone is four rolling elements to provide support, and the support situation is identical; only the rolling element serial number is different, and the total effective stiffness k eqb is equal in this case. By superimposing the effective contact stiffness of all rolling elements in the x-direction, the total effective stiffness is obtained as shown in Figure 6b. Taking

Simulation of the Stiffness, Contact Force, and Contact Deformation of Each Rolling Element in Case of the Inner Race Fault
From the simulation of Equations (23) and (27), the change of stiffness of each rolling element as it rotates with the cage and meets the inner race fault is obtained as shown in Figure 10, the key position of each rolling element as it meets the fault is shown in Table  3, and the angle of entering and exiting the load zone is the same as Figure 7 and Table 2. As can be seen from Figure 8, the No. 5 rolling element is located below the radial load F r after rotating 200 • with the cage, and it is subjected to the largest load and contact deformation in the load zone; the No. 4 and No. 6 rolling elements are next, and the No. 3 and No. 7 rolling elements are the smallest.
As can be seen from Figure 9, when the No. 5 rolling element falls into the fault position, its contact deformation δ 5 is all released, and the contact deformation δ 5 and contact force Q 5 are 0. To rebalance the external radial load F r , the No. 4 and No. 6 rolling elements, which are closest to the No. 5 rolling element, become the most important load-carrying rolling elements. Additionally, their contact force and contact deformation increase greatly, and because the angle of the No. 4 rolling element and the radial load F r and the angle of the No. 6 rolling element and the radial load F r are equal, their increase is the same. The contact force and contact deformation of the No. 3 and No. 7 rolling elements, which are far from the No. 5 rolling element, also have an equal increase, but because the angle of the No. 3 rolling element and F r and the angle of the No. 7 rolling element and F r are large, their increase is smaller.

Simulation of the Stiffness, Contact Force, and Contact Deformation of Each Rolling Element in Case of the Inner Race Fault
From the simulation of Equations (23) and (27), the change of stiffness of each rolling element as it rotates with the cage and meets the inner race fault is obtained as shown in Figure 10, the key position of each rolling element as it meets the fault is shown in Table 3, and the angle of entering and exiting the load zone is the same as Figure 7 and Table 2. As can be seen from Figure 10, the change in stiffness in the case of the inner race fault is much more complex than in the case of the outer race fault, and the position of each rolling element that falls into the fault position is related to the inner race rotation. Each rolling element meets the inner race fault every 26.5 • as it rotates with the cage, the inner race rotation angle is 66.5 • at this time, and the number of rolling elements in the load zone is 4 or 5. If the inner race fault is rotated outside the load zone, although it meets the rolling element, they are not in contact, so there is no change in the effective contact stiffness. Like the No. 3,No. 4,No. 5,No. 8,and No. 9 rolling elements in Figure 10, the contact force and contact deformation are 0 at this time.   From the simulation of Equations (5), (13), and (16), the change in the contact deformation of each rolling element in the case of the inner race fault is obtained as shown in Figure 11. From the simulation of Equations (9), (23) and (27), the change in contact force is obtained as shown in Figure 12. As can be seen in Figures 11 and 12, due to the fault rotating with the inner race, the position of each rolling element when it falls into the fault position is different, and the contact deformation and contact force increase in other rolling elements in the load zone are also different.
in Figure 11. From the simulation of Equations (9), (23), and (27), the change in contact force is obtained as shown in Figure 12. As can be seen in Figures 11 and 12, due to the fault rotating with the inner race, the position of each rolling element when it falls into the fault position is different, and the contact deformation and contact force increase in other rolling elements in the load zone are also different.   Figure 13. When the No. 6 rolling element falls into the inner race fault position after rotating 131.1 • with the cage, its contact deformation δ 6 is all released, and its contact deformation δ 6 and contact force Q 6 are 0. To balance the external radial load F r , the angle between the No. 7 rolling element and the radial load F r is minimal; that is, near the center of the load area, its contact force and contact deformation increase greatly and become the most important load-carrying rolling elements. At the same time, the No.  Figure 13. When the No. 6 rolling element falls into the inner race fault position after rotating 131.1° with the cage, its contact deformation Š is all released, and its contact deformation Š and contact force Š are 0. To balance the external radial load , the angle between the No. 7 rolling element and the radial load is minimal; that is, near the center of the load area, its contact force and contact deformation increase greatly and become the most important load-carrying rolling elements. At the same time, the No. 5 and No. 8 rolling elements are at the edge of the load zone; their contact force and contact deformation also increase, but not much. Similarly, when rolling element No. 7 falls into the inner race fault position after rotating 157.6° with the cage, the rolling elements in the load zone at this time are No. 4,No. 5,No. 6,No. 7,and No. 8, and the changes in the contact deformation and the contact force are similar and will not be repeated.

Simulation and Experiment of System Response
Fault characteristic frequency of the rolling bearing inner and outer races [11,14] = 2 1 + cos (32)

Simulation and Experiment of System Response
Fault characteristic frequency of the rolling bearing inner and outer races [11,14] f where f i is the fault characteristic frequency of the inner race; f o is the fault characteristic frequency of the outer race.
The bearing experimental data of the rotor test bench are from Case Western Reserve University (USA). The model of the rolling bearing is SKF6205, and its parameters are shown in Table 1. In the case of the inner race fault, the rotor test stand motor speed is 1721 rpm, namely, 28.68 Hz, and the theoretical calculation of the characteristic frequency f i of the inner race fault is 155.35 Hz. In the case of the outer race fault, the rotor test stand motor speed is 1725 rpm, namely, 28.75 Hz, and the characteristic frequency f o is 103.03 Hz. The dimensions of the fault are ∅0.18 mm and 0.28 mm deep, located directly below the bearing. The fault bearing is installed at the drive end, and the sampling frequency is 48 kHz. The manufacturing and installation of the bearing will cause deviations between the measured frequency and the theoretical calculation. In addition, the rotor test bench of our university was used for the experiment, as shown in Figure 14. The model of the rolling bearing is 6000, and its main parameters are shown in Table 4, where the sensor is a KISTLER 8702B25. The motor speed at the time of failure is 900 rpm, i.e., 15.0 Hz, and the theoretical calculation of f i is 66.55 Hz and f o is 38.45 Hz. The dimensional width of the failure is 1 mm and 0.7 mm deep, as shown in Figure 15.      The theoretical simulation uses the variable stiffness model (VEM) of this paper and the more widely used two-degree-of-freedom displacement excitation model (DEM) and compares and analyzes the experimental data. The simulation and experimental results of the local fault of the rolling bearing outer race are shown in Figure 16. The experimental data and the VSM in Figure 16a,c are represented by the left axis, and the DEM is represented by the right axis. In Figure 16b, the experimental data and the VSM and the DEM are represented by the left coordinate axis, and in Figure 16d, the experimental data and the DEM are represented by the left coordinate axis, and the VSM is represented by the right coordinate axis. The theoretical simulation uses the variable stiffness model (VEM) of this paper and the more widely used two-degree-of-freedom displacement excitation model (DEM) and compares and analyzes the experimental data. The simulation and experimental results of the local fault of the rolling bearing outer race are shown in Figure 16. The experimental data and the VSM in Figure 16a,c are represented by the left axis, and the DEM is represented by the right axis. In Figure 16b, the experimental data and the VSM and the DEM are represented by the left coordinate axis, and in Figure 16d, the experimental data and the DEM are represented by the left coordinate axis, and the VSM is represented by the right coordinate axis.  Figure 16a,c shows the time domain waveform in each impact response position. The VSM and the DEM match with the experimental data. The amplitude of the VSM is basically the same as the experimental data, but the DEM is more different. In general, the change in each shock response amplitude is small, and the distance between two adjacent amplitudes is 1/ . Because the overall support of the bearing is the same after each rolling element falls into the fault position, the change in total effective stiffness and the change in contact stiffness when each rolling element falls into the fault position are  the same as the experimental data, but the DEM is more different. In general, the change in each shock response amplitude is small, and the distance between two adjacent amplitudes is 1/ f o . Because the overall support of the bearing is the same after each rolling element falls into the fault position, the change in total effective stiffness k eq and the change in contact stiffness when each rolling element falls into the fault position are also the same, as shown in Figures 6b and 7, resulting in the same amplitude of the system response under radial load. Figure 16b,d shows the envelope spectrum. There are significant outer race local fault characteristic frequencies f accompanied by frequency multiplication in the figure. For example, the characteristic frequencies of bearing 6000 are 38.5 Hz, 76.9 Hz, etc., and the characteristic frequencies of bearing SKF6205 are 103.5 Hz, 206.9 Hz, 310.5 Hz, etc. At each fault characteristic frequency position, the VSM and the DEM match with the experimental data, but the amplitude of the DEM in Figure 16b is too large, while the amplitude of the VSM in Figure 16d is slightly larger. Therefore, the results of the VSM are in better agreement with the experimental results and are superior to the DEM.
The simulation and experimental results of the local fault of the rolling bearing inner race are shown in Figure 17. The experimental data in Figure 17 are represented by the left axis, and the VSM and the DEM are represented by the right axis.  Figure 17a,c shows the time domain waveform in each impact response position. The VSM and the DEM match the experimental data. The amplitudes of the DEM and the VSM are larger than the experimental data, but the DEM is more different in Figure 17a. In general, each shock response amplitude varies drastically, and the distance between two adjacent amplitudes is 1/ . When each rolling element rotates with the cage and meets the inner race fault, the total effective stiffness is changed, the change in the total effective stiffness and the change in the contact stiffness when each rolling element falls into the fault position are different, as shown in Figures 6c and 8, and the magnitude of the change of total effective stiffness is modulated by the inner race rotation, resulting in a change in the system response amplitude as well.  Figure 17a,c shows the time domain waveform in each impact response position. The VSM and the DEM match the experimental data. The amplitudes of the DEM and the VSM are larger than the experimental data, but the DEM is more different in Figure 17a. In general, each shock response amplitude varies drastically, and the distance between two adjacent amplitudes is 1/ f i . When each rolling element rotates with the cage and meets the inner race fault, the total effective stiffness k eq is changed, the change in the total effective stiffness k eq and the change in the contact stiffness when each rolling element falls into the fault position are different, as shown in Figures 6c and 8, and the magnitude of the change of total effective stiffness k eq is modulated by the inner race rotation, resulting in a change in the system response amplitude as well. Figure 17b,d shows the envelope spectrum. There is a significant inner race rotation frequency f s accompanied by frequency multiplication as well as fault characteristic frequency f i and its frequency multiplication in the figure. The side frequency with obvious amplitude is distributed on both sides of the fault characteristic frequency f i and its frequency multiplication. For example, the characteristic frequencies of bearing 6000 are 66. 6 31 Hz are the inner race rotation frequencies. Therefore, the theoretical simulation and experimental results agree that there is a significant modulation of the response amplitude in the case of the inner race fault. At each fault characteristic frequency and its side frequency and rotation frequency positions, the VSM and the DEM match with the experimental data, but the amplitude of the DEM is significantly larger than that of the VSM and experimental data. Therefore, the results of the VSM are in better agreement with the experimental results and are better than the DEM.
The degree of agreement between the simulated and experimental data can be expressed in terms of average coherence: where N is the number of time domain peak samples (N = 30), A Ei is the peak of experimental data, and A i is the peak of the VSM or DEM simulation data. From Table 5, the amplitude of the VSM is more consistent with the experimental data, and the consistency of the outer ring is the best, with a minimum of 49%; the inner ring is slightly worse, with a minimum of 241%. However, both are significantly better than the DEM, indicating that the VSM is more consistent with the reality. In summary, when the rolling elements fall into the fault positions of the inner and outer races, the total effective stiffness of the system will change abruptly. Therefore, to re-balance the external radial load, the contact deformation and contact force of each load-carrying rolling element in the load zone change, and finally, this causes the dynamic system vibration. As the fault position of the outer race is unchanged relative to the radial load and the position of the load zone, the amplitude of the total effective stiffness reduction of the system is alike when each rolling element falls into the fault position with the rotation of the cage. Therefore, the amplitude of the outer loop response is alike, and there is a significant characteristic frequency f o of outer race faults and its frequency multiplication in the fault characteristic spectrum. However, the position of the inner race fault relative to the radial load and the load zone changes periodically with the rotation of the inner race; the rolling elements in the load zone will change the total effective stiffness of the system due to the inner race fault. Therefore, the total effective support stiffness of the system is affected by the rotation of the inner race, which further leads to the amplitude of the inner loop response also being affected. In the fault characteristic spectrum of the inner race, there is a significant rotation frequency of the inner race f s accompanied by frequency multiplication as well as a fault characteristic frequency of the inner race f i and its frequency multiplication and side frequency.

Conclusions
By setting up the dynamics modeling of the rolling bearing variable stiffness system with local fault, this paper researches the change in effective contact stiffness, contact deformation, contact forces, and total effective stiffness of the rolling elements. Some conclusions are as follows.
(1) Based on the retention factor of contact deformation, the single-degree-of-freedom variable stiffness model of rolling bearings with fault is proposed, and the dynamics modeling of the variable stiffness of rolling bearings with fault is established. (2) The contact stiffness of the rolling elements abruptly decreases when the rolling elements fall into the fault position. The contact deformation and contact force of the load-carrying rolling elements in the load zone increases, rebalancing the external radial load while causing a sudden reduction in the total effective stiffness, resulting in the vibration of the system. (3) When different rolling elements fall into the outer race fault position, the change in the total effective stiffness and the system response are equal in magnitude. Additionally, there is significant outer race fault characteristic frequency and its frequency multiplication in the fault characteristic spectrums. When different rolling elements fall into the inner race fault position, the total effective stiffness is modulated by the inner race rotation and varies dramatically, resulting in the amplitude of the system time domain vibration response also being modulated by the inner race rotation and varying dramatically. Additionally, there are significant inner race rotational frequencies and their frequency multiplications, inner race fault characteristic frequencies and their frequency multiplication, and side frequency in the fault characteristic spectrums. (4) The VSM is more consistent with the experiment and provides some theoretical basis for the effective diagnosis of rolling bearing faults. The contact deformation retention factor is only for rectangular or circular faults, which has some limitations. In the subsequent research, the VSM will be applied to other types of rolling bearings to expand the application scope of the VSM. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

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