Mechanism and Characteristics of Global Varying Compliance Parametric Resonances in a Ball Bearing

: Varying compliance (VC) is an unavoidable form of parametric excitation in rolling bearings and can a ﬀ ect the stability and safety of the bearing and its supporting rotor system. To date, we have investigated VC primary resonance in ball bearings, and in this paper other parametric VC resonance types are addressed. For a classical ball bearing model with Hertzian contact and clearance nonlinearities between the rolling elements and raceway, the harmonic balance and alternating frequency / time domain (HB–AFT) method and Floquet theory are adopted to analyze the VC parametric resonances and their stabilities. It is found that the 1 / 2-order subharmonic resonances, 2-order superharmonic resonances, and various VC combination resonances, such as the 1-order and 2-order summed types, can be excited, thus resulting in period-1, period-2, period-4, period-8, period-35, quasi-period, and even chaotic VC motions in the system. Furthermore, the bifurcation and hysteresis characteristics of complex VC resonant responses are discussed, in which cyclic fold, period doubling, and the second Hopf bifurcation can occur. Finally, the global involution of VC resonances around bearing clearance-free operations (i.e., adjusting the bearing clearance to zero or one with low interference) are provided. The overall results extend the investigation of VC parametric resonance cases in rolling bearings.


Introduction
Faraday [1] previously recognized the phenomenon of parametric resonance when he studied the response frequency of surface waves in a fluid-filled cylinder with vertical excitations. Researchers have subsequently carried out extensive investigations of parametric resonances for various nonlinear vibration systems in the field of neutrino oscillations [2], quantum fields [3], ship dynamics [4], Brownian dynamics of water droplets [5], bubbles oscillating in water [6], micro-electro mechanical systems [7,8], and rotor dynamics [9,10]. Therefore, the topic of parametric resonance is a basic scientific problem in the field of nonlinear dynamics.
Rolling bearings are widely used in rotating machineries due to their advantages of lower friction compared to sliding bearings and ease with which they can be selected from a manufacturer's catalogue [11]. There are two major effects of a healthy rolling bearing on machine vibration [12,13]: spring characteristics and the excitation force on the supported structure. In the dynamic design of the bearing-rotor system, the excitation characteristics of rolling bearings should be given adequate consideration. Varying compliance (VC) excitation is induced when rolling elements move along the bearing raceway with an external load, which results in time-varying bearing stiffness in the rotor system. Thus, VC excitation is an inevitable case of parametric excitation [9,14].
Sunnersjö [9] found that VC excitation can result in periodic and irregular non-periodic vibrations in cylindrical roller bearings using digital simulation and experiments. Fukata et al. [14] identified that VC vibration may exhibit subharmonic, quasi-periodic, chaotic, and beat motions at critical speeds, considering a two-degree-of-freedom ball bearing model with Hertzian contact nonlinearity. Mevel and Guyader [15,16] studied the routes to chaos in ball bearings and noted the chaos due to period doubling in the resonant speed range. Tomović [17] found that the VC vibration amplitude of rolling bearings is directly related to the bearing deflection in two specific boundary positions of the rolling elements, and a simplified rolling bearing model for description of VC vibrations was presented. Regarding VC hysteresis characteristics, Sankaravelu et al. [18] found a hysteretic amplitude-frequency response in ball bearing VC vibration using an early shooting method. Ghafari et al. [19] investigated the hysteresis behaviors using static bifurcation analysis and provided closed-form expressions of VC periodic solutions using the generalized averaging method. Based on Hertzian contact resonance [20,21], Zhang et al. [22][23][24][25] investigated the mechanism, influence, and suppression of VC resonant hysteresis in ball bearings. Theoretical and experimental results showed that the primary contact resonance exhibits soft spring behavior in the direction exerting constant force, and the coexistence of soft and hard spring characteristics in the direction with zero external force. In addition, Zhang et al. [26] noted a case of complex coupling VC resonance caused by the superposition of VC parametric primary resonance and 1:2 internal resonance. The above research is mainly focused on the primary VC resonant responses. However, it should be noted that a small parametric excitation may produce remarkable responses, even if the excitation frequency significantly differs from the linear natural frequencies of the system [27]. Other parametric resonance types may also exist in a parametric system, such as superharmonic, subharmonic, and combined cases.
Rapid, flexible, and automated development requires higher stability and safety characteristics in rotating machineries, which promotes the study of super-and subharmonic resonances, and their combinations, in rolling bearing-rotor systems. Tiwari et al. [28,29], Harsha et al. [30,31], and Bai et al. [32] presented the amplitude and stability of VC responses for a rolling bearing-rotor system in a large speed range; however, the mechanism of the resonant amplitude has been little addressed. Ehrich [33] observed high-order subharmonic vibration responses under low damping and strong nonlinearity conditions in high speed rotating machinery with asymmetric bearing supports, based on the numerical integration of a finite difference formulation. Yoshida et al. [34] showed the soft and hard hysteretic 1/2-order subharmonic resonances in different degrees of freedom of a ball bearing-rotor system using the multi-scale method. Although this study also found combination resonances between the 1/2-order subharmonic resonance intervals, this has not been analyzed in depth. Bai et al. [35] established a 6DOF nonlinear ball bearing-rotor model and found that the occurrence of 2-order subharmonic resonance is provoked by the coupling of Hertzian contact and internal clearance nonlinear factors of ball bearings. Wu et al. [36] noted that the subharmonic hysteretic resonance become significant and the resonance amplitude moved to lower frequencies as the assembly gap increased for an angular contact high-speed ball bearing system. However, these studies [33][34][35][36] relate to the case of external excitations from unbalanced rotors, for which a small excitation produces a large response only if the frequency of the external excitation is close to that of the linear natural frequency of the system, and it is essentially different from that of the parametric excitation system. Therefore, a global mechanism analysis of VC parametric resonances in rolling bearings is necessary.
Bearing clearance is a basic parameter of rolling bearings, and has a significant effect on bearing life, installation, and bearing stiffness [11,12]. Yammamoto [37] previously showed that the resonant frequency and amplitude can be reduced to some extent as the clearance adjusts. Oswald et al. [38] Appl. Sci. 2020, 10, 7849 3 of 28 noted that the bearing life declines gradually with positive clearance and rapidly with increasing negative clearance, and the life can be maximized for a small negative operating clearance. In the field of bearings, it is generally considered that adverse vibration and noise behaviors can be reduced by adjusting the ball bearings to zero clearance [39]. In particular, a clearance-free operation is usually desirable for precision or high-speed machinery, such as machine tools and turbines [11]. Therefore, the study of VC parametric resonances in clearance-free operations for rolling bearing systems is valuable.
The harmonic balance and alternating frequency/time domain (HB-AFT) method is an effective means to obtain the steady-state harmonic solutions of a general nonlinear system. It can be used to quickly obtain the frequency domain information of the nonlinear term for the harmonic balance process of the AFT discretization technology [40]. Kim et al. [41] and Guskov et al. [42] developed the generalized HB-AFT method to analyze the quasi-periodic response characteristics of multiple degrees-of-freedom and multi-frequency excitation systems. Zhang et al. [22] embedded continuation and Hsu's Floquet stability analysis technology into the HB-AFT method, which can automatically track the periodic response and bifurcation characteristics of nonlinear systems. Recently, Chen et al. [43] developed a harmonic balance method combined with the model incremental transfer matrix of the system, which is highly convenient for dealing with systems with multiple degrees-of-freedom.
The motivation of the present work is to provide a better understanding of the VC parametric resonances in a ball bearing. The mechanisms of various new cases of parametric resonances are provided mainly using the HB-AFT method with the Floquet theory, and their detailed response characteristics, such as time series, orbit, stability, and hysteresis, are investigated. The global amplitude-frequency responses around the bearing clearance-free operations are examined in detail and predicted.

Description of Ball Bearing VC Vibration
Considering the nonlinearities of VC parametric excitation, bearing clearance, and the Hertzian contact between the rolling elements and raceway, a classic two degrees-of-freedom ball bearing model is studied here. This model has been verified for qualitative and quantitative analysis of the response characteristics of radial VC motions in ball bearings [14][15][16]24], and the JIS6306 specification shown in Table 1 is adopted in this article. As shown in Figure 1, the instant angular location of the ith ball is denoted as: with ball passing frequency (i.e., the cage angular velocity) of: Then, the contact deformation between the ith ball and bearing raceways is expressed as: where ω s and 2δ 0 denote the rotor shaft speed and the radial internal bearing clearance, respectively. Assuming only Hertzian contact deformations are considered, the nonlinear restoring forces F x and F y are: where G[·] is a Heaviside function, taking 1 if δ i ≥ 0 represents the condition of contact, or else taking 0, which denotes the contact loss between the rolling element and raceway. The VC motion of a rigid Jeffcott rotor-ball bearing system can be expressed as: Due to the Hertzian contact deformation relationship, the system has hard spring characteristics [44][45][46]. The approximate linearized static stiffness of System (5) in the direction of static load W s can be expressed as [46]: where δ Ws denotes the total deformation along the direction of W s . Equation (6) is clearly not applicable for the estimation of the dynamic stiffness of a ball bearing system during operation. The dynamic stiffness characteristics of rolling bearings have a significant impact on the nonlinear response characteristics of their supporting systems [12,22,26]. According to the restoring forces of Equation (4), the radial linearized dynamic stiffness of the system can be approximated as [15]: Then, the dynamic resonant frequencies of the system can be expressed as: and the equivalent resonant frequencies ω x0 , ω y0 can be estimated from the arithmetic mean of Equation (8) [22]. The VC parametric excitation frequency occurs due to the ball passage frequency, which can be expressed as: Many cases of resonance might exist in a parametric excitation system with many degrees-of-freedom [27,47,48], and the conditions of possible excited VC parametric primary, subharmonic, and superharmonic resonances are: where the VC parametric primary resonance can occur when n(m) = 1; 1/h-order subharmonic resonances can emerge when n(m) = h; the h-order superharmonic resonance may be excited when n(m) = 1/h. In addition, VC parametric combination resonance might occur in (but is not limited to) the following conditions: Ω VC ≈ n · ω xx ± m · ω yy ± 1/ j · ω xx ± 1/k · ω yy (11) Figure 1. Schematic diagram of (a) the ball bearing system, and (b) its two degrees-of-freedom spring model.

Methodology
For a periodically excited nonlinear vibration system, a periodic motion is called the period-M motion if the period of the response is M times the excited period, and the M-order subharmonic resonance interval must contain period-M solutions. The process of deriving the VC period-M solution and the stability analysis for the studied System (5) is given below.

HB-AFT Method for Obtaining VC Period-M Solution
To improve the accuracy and efficiency of obtaining the VC period-M solution, one can introduce a dimensionless displacement X = [X, Y] T = [x/δ0， y/δ0] T , and dimensionless time τ = Nb·Ω·t/M; then, Equation (5) can be expressed as: F X X (12) Here, let the control parameter λ = Ω for tracing the amplitude-frequency response curves, and for System (12) the period TVC of the VC period-M solution is 2π. Then, X(τ) and F(·) can be expressed as: (13) According to the process of harmonic balance, insert Equation (13) into Equation (12) and obtain the following 4K + 2 algebraic relationships: ( , , ) 0 λ = g P Q (14) where the harmonic coefficients of X(τ) and F(·) are :   T  T  0  1  1  2  2   0  1  1  2  2 , , , , , , , , , , , , , , To obtain P from Equation (14), it should first formulate Q by P, because Equation (14) contains 8K + 4 unknown harmonic coefficients of P and Q, and this process can be completed by the following AFT technology.
The discrete sample of X(τ) is given by the inverse discrete Fourier transform as: (16) and the discrete time domain sample of Equation (12) is: Schematic diagram of (a) the ball bearing system, and (b) its two degrees-of-freedom spring model.

Methodology
For a periodically excited nonlinear vibration system, a periodic motion is called the period-M motion if the period of the response is M times the excited period, and the M-order subharmonic resonance interval must contain period-M solutions. The process of deriving the VC period-M solution and the stability analysis for the studied System (5) is given below.

HB-AFT Method for Obtaining VC Period-M Solution
To improve the accuracy and efficiency of obtaining the VC period-M solution, one can introduce a dimensionless displacement X = [X, Y] T = [x/δ 0 , y/δ 0 ] T , and dimensionless time τ = N b ·Ω·t/M; then, Equation (5) can be expressed as: Here, let the control parameter λ = Ω for tracing the amplitude-frequency response curves, and for System (12) the period T VC of the VC period-M solution is 2π. Then, X(τ) and F(·) can be expressed as: According to the process of harmonic balance, insert Equation (13) into Equation (12) and obtain the following 4K + 2 algebraic relationships: where the harmonic coefficients of X(τ) and F(·) are: To obtain P from Equation (14), it should first formulate Q by P, because Equation (14) contains 8K + 4 unknown harmonic coefficients of P and Q, and this process can be completed by the following AFT technology.
Appl. Sci. 2020, 10, 7849 6 of 28 The discrete sample of X(τ) is given by the inverse discrete Fourier transform as: 16) and the discrete time domain sample of Equation (12) is: Then, by discrete Fourier transform, Q can be expressed as: N is the numbers of samplings in the time domain, and At this moment, Q is formulated by P as shown in Equation (18). Taking P as an unknown variable, one can obtain the harmonic solution of X(τ) by solving the nonlinear algebraic Equation (14). In order to avoid the singular case at the turning point and track P automatically, the arc-length continuation is introduced by transforming Equation (14) into following form [22]: where s denotes the arc-length parameter of curve (P(s), λ(s)).

Floquet Stability Analysis of VC Period-M Solution
In the following section, the detailed process of stability analysis for VC period-M solution U * (τ) obtained in Section 3.1 will be given by using Floquet stability theory [51].
Introducing U = [X, X , Y, Y ] T = [X 1 , X 2 , X 3 , X 4 ] T , transform System (12) into: where For Equation (20), superpose ∆U to perturb the VC period-M solution U * (τ), i.e., and retain only linear terms in the disturbance as: Then, the local linear stability of U * (τ) can be provided by the stability information of Equation (22). As U * (τ) is a period-M solution of the System (20), ∂H(τ, U * (τ))/∂U * (τ) is also a periodic variable matrix, and the stability of the VC period-M solution U * (τ) can be determined by the eigenvalues of the Floquet monodromy matrix of Equation (22). The Hsu method [22] is used to calculate the Floquet multiplier λ m and, according to the three means of λ m passing the unit circle of a complex plane, the period-M solution U * (τ) can lose stability as follows [51]: • A transcritical, symmetry-breaking, or cyclic-fold bifurcation may originate when the leading multiplier crosses through the unit circle at +1. Herein, the cyclic-fold bifurcation emerges only at the turning point of a solution branch.

•
A period-doubling bifurcation occurs when the leading multiplier crosses the unit circle at −1. • A secondary Hopf bifurcation occurs when a complex pair of multipliers moves out of the unit circle.

Results
4.1. The VC 2-Order Superharmonc Resonance, 1/2-Order Subharmonic Resonance and the Case of ω x0 + ω y0 = Ω vc Combination Resonance For δ 0 = 1.5 µm, Figure 2 shows the VC period-1 and period-2 solution branches traced by the HB-AFT method with arc-length continuation. Predicting from Equations (8) and (9), the equivalent resonant frequencies of the two degrees of freedom of System (5) are about ω x0 = 1908.14 rad/s and ω y0 = 1297.36 rad/s. Thus, the amplitudes in the R 1 and R 2 frequency ranges, including ω x0/ N b = 238.52 rad/s and ω y0/ N b = 162.17 rad/s, are the VC primary resonances excited in the x and y directions of the system, respectively. These two primary resonances exhibit soft spring and soft-hard spring characteristics in the vertical and horizontal directions, respectively, which agrees with previous theoretical and experimental studies [22,24]. Moreover, 1/2-order subharmonic resonances in two degrees of freedom are excited in R 3 and R 4 intervals, where the subharmonic resonant amplitudes are around 2·ω x0/ N b = 477.04 rad/s and 2·ω y0/ N b = 324.34 rad/s. For δ0 = 1.5 μm, Figure 2 shows the VC period-1 and period-2 solution branches traced by the HB-AFT method with arc-length continuation. Predicting from Equations (8) and (9)  As shown in Figure 3a, the resonant amplitude in the R4 range is excited by supercritical perioddoubling bifurcations at points A1 (see Table 2) and A2, where the leading Floquet multipliers cross the unit circle at −1. The response characteristics of the period-1 and period-2 motions before and after period-doubling transitions at point A1 are presented in Figures 4 and 5, and the response amplitude of y(τ) clearly increases after the bifurcation. Then, the leading Floquet multiplier of the stable VC period-2 solution crosses the +1 axis at A3 and A4 (see Table 3), and cyclic fold bifurcation occurs at these two points. This causes significant resonant hysteresis and jumping behaviors in the 1/2-order subharmonic resonance in the horizontal direction. For example, as the control parameter Ω increases, the stable period-2 branch suddenly jumps to the stable period-1 branch from point A4 to A5. For the 1/2-order subharmonic resonance as Ω around 2·ωx0/Nb = 477.04 rad/s (i.e., R3 range), as shown in Figure 3b, supercritical and subcritical period-doubling bifurcations trigger the VC period-2 resonant responses at points B1 and B2, respectively. At this time, the subharmonic resonant frequency-response curve bends to the left, and this soft spring characteristic is consistent with the hysteresis characteristics of the primary resonance in the vertical direction as Ω around ωx0/Nb = 238.52 rad/s. As shown in Figure 3a, the resonant amplitude in the R 4 range is excited by supercritical period-doubling bifurcations at points A 1 (see Table 2) and A 2 , where the leading Floquet multipliers cross the unit circle at −1. The response characteristics of the period-1 and period-2 motions before and after period-doubling transitions at point A 1 are presented in Figures 4 and 5, and the response amplitude of y(τ) clearly increases after the bifurcation. Then, the leading Floquet multiplier of the stable VC period-2 solution crosses the +1 axis at A 3 and A 4 (see Table 3), and cyclic fold bifurcation occurs at these two points. This causes significant resonant hysteresis and jumping behaviors in the 1/2-order subharmonic resonance in the horizontal direction. For example, as the control parameter Ω increases, the stable period-2 branch suddenly jumps to the stable period-1 branch from point A 4 to A 5 . For the 1/2-order subharmonic resonance as Ω around 2·ω x0/ N b = 477.04 rad/s (i.e., R 3 range), as shown in Figure 3b, supercritical and subcritical period-doubling bifurcations trigger the VC period-2 resonant responses at points B 1 and B 2 , respectively. At this time, the subharmonic resonant frequency-response curve bends to the left, and this soft spring characteristic is consistent with the hysteresis characteristics of the primary resonance in the vertical direction as Ω around ω x0/ N b = 238.52 rad/s.     Table 3. Period-2 motion Floquet multipliers λm around A4 (counterclockwise, see Figure 3a).  Moreover, the stable period-2 solution branch undergoes supercritical period-doubling bifurcation at B3, generating a stable period-4 branch B3-B4. Then, a period-8 solution branch B4-B5 is excited by supercritical period-doubling bifurcation at point B4., where point B5 is the cycle-doubling bifurcation point. In the above bifurcation process, the strongly excited resonant period-2 and period-4 motion modes in the x direction become saturated at B3 and B4 points, and the energy spills over   Table 3. Period-2 motion Floquet multipliers λ m around A 4 (counterclockwise, see Figure 3a).  Moreover, the stable period-2 solution branch undergoes supercritical period-doubling bifurcation at B 3 , generating a stable period-4 branch B 3 -B 4 . Then, a period-8 solution branch B 4 -B 5 is excited by supercritical period-doubling bifurcation at point B 4 ., where point B 5 is the cycle-doubling bifurcation point. In the above bifurcation process, the strongly excited resonant period-2 and period-4 motion modes in the x direction become saturated at B 3 and B 4 points, and the energy spills over into the period-4 and period-8 motion modes in the y direction as Ω decreases, respectively. This typical saturation phenomenon has been discussed in [26,27]. In short, the 1/2-order subharmonic resonances are aroused by the period-doubling bifurcation instability, and they contain complex hysteresis characteristics. Finally, as shown in Figure 6, the above results agree well with the numerical bifurcation diagrams from the local-maxima method [52].
into the period-4 and period-8 motion modes in the y direction as Ω decreases, respectively. This typical saturation phenomenon has been discussed in [26,27]. In short, the 1/2-order subharmonic resonances are aroused by the period-doubling bifurcation instability, and they contain complex hysteresis characteristics. Finally, as shown in Figure 6, the above results agree well with the numerical bifurcation diagrams from the local-maxima method [52]. Furthermore, another resonant amplitude is excited around ωx0/Nb + ωy0/Nb = 400.69 rad/s (i.e., the resonant range R5 between intervals R3 and R4 presented in Figure 2), and this case of resonance results in strong coupling to the two degrees-of-freedom of the system as shown in Figure 6. In terms of the bifurcation characteristics concerned, as a pair of complex conjugate Floquet multipliers move out of the unit circle (see Table 4), subcritical and supercritical secondary Hopf bifurcations occur at points B6 and B7, respectively, leading quasi-period responses excited at point B7. The response characteristics of the quasi-period motion at Ω = 400 rad/s is presented in Figure 7; as there are incommensurable frequency components, the Poincare mapping of the response orbit tends to a closed point set. Herein, the values of the main incommensurable dimensionless VC frequencies p and q are about 0.5936 and 0.4074, respectively, and p + q = 1, where the value 1 is the dimensionless VC excitation frequency ΩVC = Nb·Ω. Moreover, for the original time scale, the corresponding main frequency components of the response are Nb·Ω·p ≈ 1899.52 rad/s and Nb·Ω·q ≈ 1303.68 rad/s, which are close to the predicted equivalent resonant frequencies ωx0 = 1908.14 rad/s and ωy0 = 1297.36 rad/s of the two degrees-of-freedom of System (5). This indicates that the case of ωx0 + ωy0 = ΩVC combination resonance (i.e., the 1-order summed type of combination resonance) [39,53] occurs in range R5, and quasi-period resonant responses are excited. In general, the typical combination resonance of multiple degrees-of-freedom systems often causes quasi-period motions [27,54,55]. Therefore, it can be considered that the combination resonance is one of the physical mechanisms of the generation of VC quasi-period motions in rolling bearing systems. In addition, due to the nonlinearities from Furthermore, another resonant amplitude is excited around ω x0 /N b + ω y0/ N b = 400.69 rad/s (i.e., the resonant range R 5 between intervals R 3 and R 4 presented in Figure 2), and this case of resonance results in strong coupling to the two degrees-of-freedom of the system as shown in Figure 6. In terms of the bifurcation characteristics concerned, as a pair of complex conjugate Floquet multipliers move out of the unit circle (see Table 4), subcritical and supercritical secondary Hopf bifurcations occur at points B 6 and B 7 , respectively, leading quasi-period responses excited at point B 7 . The response characteristics of the quasi-period motion at Ω = 400 rad/s is presented in Figure 7; as there are incommensurable frequency components, the Poincare mapping of the response orbit tends to a closed point set. Herein, the values of the main incommensurable dimensionless VC frequencies p and q are about 0.5936 and 0.4074, respectively, and p + q = 1, where the value 1 is the dimensionless VC excitation frequency Ω VC = N b ·Ω. Moreover, for the original time scale, the corresponding main frequency components of the response are N b ·Ω·p ≈ 1899.52 rad/s and N b ·Ω·q ≈ 1303.68 rad/s, which are close to the predicted equivalent resonant frequencies ω x0 = 1908.14 rad/s and ω y0 = 1297.36 rad/s of the two degrees-of-freedom of System (5). This indicates that the case of ω x0 + ω y0 = Ω VC combination resonance (i.e., the 1-order summed type of combination resonance) [39,53] occurs in range R 5 , and quasi-period resonant responses are excited. In general, the typical combination resonance of multiple degrees-of-freedom systems often causes quasi-period motions [27,54,55]. Therefore, it can be considered that the combination resonance is one of the physical mechanisms of the generation of VC quasi-period motions in rolling bearing systems. In addition, due to the nonlinearities from Hertzian contact and bearing clearance, this combination resonance has soft hysteresis characteristics, and the quasi-period response jumps to the VC period-1 motion branch as Ω decreases to 388.20 rad/s.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 28 Hertzian contact and bearing clearance, this combination resonance has soft hysteresis characteristics, and the quasi-period response jumps to the VC period-1 motion branch as Ω decreases to 388.20 rad/s.
For certain applications, such as high-speed or high-precision machines, to increase bearing stiffness and reduce the run-out and noise due to elastic deformation or clearance [11,39], it is possible to eliminate clearance and introduce small interference (negative clearance) in the bearings. Figure 8 presents frequency-response curves of the VC parametric response for δ0 = −1.5 μm. In this case, the equivalent resonant frequencies of the two degrees-of-freedom of System (5) are about ωx0 = 1959.49 rad/s and ωy0 = 2105.55 rad/s. Compared to the case of δ0 = 1.5 μm, the resonant excitation frequency is higher in the horizontal direction than in the vertical direction, and this frequency exchange, due to the adjusting of bearing clearance, was previously illustrated by Fang [56] and Zhang [25], respectively. Herein, the VC primary resonances as Ω rotates around ωx0/Nb = 244.94 rad/s and ωy0/Nb = 263.19 rad/s, 1/2-order subharmonic resonances as Ω rotates around 2·ωx0/Nb = 489.88 rad/s and 2·ωy0/Nb = 526.38 rad/s, and the case of ωx0 + ωy0 = ΩVC combination resonance as Ω rotates around ωx0/Nb + ωy0/Nb = 508.13 rad/s are also excited. The difference compared to the case of δ0 = 1.5 μm is that, as shown in Figure 9, the primary and 1/2-order subharmonic resonances of the system in vertical For certain applications, such as high-speed or high-precision machines, to increase bearing stiffness and reduce the run-out and noise due to elastic deformation or clearance [11,39], it is possible to eliminate clearance and introduce small interference (negative clearance) in the bearings. Figure 8 presents frequency-response curves of the VC parametric response for δ 0 = −1.5 µm. In this case, the equivalent resonant frequencies of the two degrees-of-freedom of System (5) are about ω x0 = 1959.49 rad/s and ω y0 = 2105.55 rad/s. Compared to the case of δ 0 = 1.5 µm, the resonant excitation frequency is higher in the horizontal direction than in the vertical direction, and this frequency exchange, due to the adjusting of bearing clearance, was previously illustrated by Fang [56] and Zhang [25], respectively. Herein, the VC primary resonances as Ω rotates around ω x0/ N b = 244.94 rad/s and ω y0/ N b = 263.19 rad/s, 1/2-order subharmonic resonances as Ω rotates around 2·ω x0/ N b = 489.88 rad/s and 2·ω y0/ N b = 526.38 rad/s, and the case of ω x0 + ω y0 = Ω VC combination resonance as Ω rotates around ω x0 /N b + ω y0/ N b = 508.13 rad/s are also excited. The difference compared to the case of δ 0 = 1.5 µm is that, as shown in Figure 9, the primary and 1/2-order subharmonic resonances of the system in vertical directions (i.e., the R 1 and R 3 ranges) exhibit soft-hard and hard spring characteristics, respectively, and the bifurcations excited by the combination resonance at the points C 1 and C 2 are all supercritical Hopf types. These bifurcation characteristics are clearly illustrated by the numerical bifurcation diagrams in Figure 10. In addition, the VC 2-order superharmonic resonances of the two degrees-of-freedom of the system are strongly excited as Ω rotates around 1/2·ω x0/ N b = 122.47 rad/s and 1/2·ω y0/ N b = 131.60 rad/s (i.e., the R 6 and R 7 ranges in Figure 8), respectively. Figure 11 presents the response characteristics of the 2-order superharmonic resonance in the y direction of the system, in which the superharmonic frequency component dominates the period-1 response of y(τ).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 28 and the bifurcations excited by the combination resonance at the points C1 and C2 are all supercritical Hopf types. These bifurcation characteristics are clearly illustrated by the numerical bifurcation diagrams in Figure 10. In addition, the VC 2-order superharmonic resonances of the two degrees-offreedom of the system are strongly excited as Ω rotates around 1/2·ωx0/Nb = 122.47 rad/s and 1/2·ωy0/Nb = 131.60 rad/s (i.e., the R6 and R7 ranges in Figure 8), respectively. Figure 11 presents the response characteristics of the 2-order superharmonic resonance in the y direction of the system, in which the superharmonic frequency component dominates the period-1 response of y(τ).  Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 28 and the bifurcations excited by the combination resonance at the points C1 and C2 are all supercritical Hopf types. These bifurcation characteristics are clearly illustrated by the numerical bifurcation diagrams in Figure 10. In addition, the VC 2-order superharmonic resonances of the two degrees-offreedom of the system are strongly excited as Ω rotates around 1/2·ωx0/Nb = 122.47 rad/s and 1/2·ωy0/Nb = 131.60 rad/s (i.e., the R6 and R7 ranges in Figure 8), respectively. Figure 11 presents the response characteristics of the 2-order superharmonic resonance in the y direction of the system, in which the superharmonic frequency component dominates the period-1 response of y(τ).   Figure 12 shows the power spectra of the dynamic displacement parameter (DDP) [9]: DDP = x peak-to-peak + y peak-to-peak , where x peak-to-peak = x(t) max − x(t) min , and y peak-to-peak = y(t) max − y(t) min . To reflect the hysteresis characteristics, here the larger DDP values are selected from the two DDP values obtained as the controlled parameter Ω runs up and down. It is clear that the system contains many types of VC parametric resonances. From Section 4.1, some resonance types of the two degrees-of-freedom of the system have been confirmed, such as the primary resonances (i.e., R 1 and R 2 ), 1/2-order subharmonic resonances (i.e., R 3 and R 4 ), 2-order superharmonic resonances (i.e., R 6 and R 7 ), and the case of ω x0 + ω y0 = Ω VC combination resonance (i.e., R 5 ). Moreover, the cases of 1/2·(ω x0 + ω y0 ) = Ω VC combination resonance (i.e., the 2-order summed type of combination resonance [39,57]) and 1/2·(1/4·ω x0 + ω y0 ) = Ω VC combination resonance (labelled as R 8 and R 9 , respectively) were found in our recent work [25]. As shown in Figure 13, these resonant locations can be predicted well by the equivalent resonant frequencies ω x0 and ω y0 from the arithmetic mean of Equation (8) at Ω = 700 rad/s, where we assume that the system is far from any resonance when Ω takes 700 rad/s to achieve an almost linearized dynamic stiffness (7). In addition, there are the other cases of 1/2·(1/2·ω x0 + ω y0 ) = Ω VC and 1/2·(ω x0 + 1/2·ω y0 ) = Ω VC combination resonances (labelled as R 10 and R 11 , respectively) excited as the bearing clearance δ 0 increases beyond 2.0 µm.

Global Response Characteristics of VC Resonances around the Bearing Clearance-Free Operations
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 28 Figure 11. For δ0 = −1.5 μm, response characteristics of the VC superharmonic motion when Ω = 131.6 rad/s, (a) orbit (black line) and its Poincare mapping (red dots), (b) time series and their frequency spectra. Figure 12 shows the power spectra of the dynamic displacement parameter (DDP) [9]: DDP = xpeak-to-peak + ypeak-to-peak, where xpeak-to-peak = x(t)max − x(t)min, and ypeak-to-peak = y(t)max − y(t)min. To reflect the hysteresis characteristics, here the larger DDP values are selected from the two DDP values obtained as the controlled parameter Ω runs up and down. It is clear that the system contains many types of VC parametric resonances. From Section 4.1, some resonance types of the two degrees-of-freedom of the system have been confirmed, such as the primary resonances (i.e., R1 and R2), 1/2-order subharmonic resonances (i.e., R3 and R4), 2-order superharmonic resonances (i.e., R6 and R7), and the case of ωx0 + ωy0 = ΩVC combination resonance (i.e., R5). Moreover, the cases of 1/2·(ωx0 + ωy0) = ΩVC combination resonance (i.e., the 2-order summed type of combination resonance [39,57]) and 1/2·(1/4·ωx0 + ωy0) = ΩVC combination resonance (labelled as R8 and R9, respectively) were found in our recent work [25]. As shown in Figure 13, these resonant locations can be predicted well by the equivalent resonant frequencies ωx0 and ωy0 from the arithmetic mean of Equation (8) at Ω = 700 rad/s, where we assume that the system is far from any resonance when Ω takes 700 rad/s to achieve an almost linearized dynamic stiffness (7). In addition, there are the other cases of 1/2·(1/2·ωx0 + ωy0) = ΩVC and 1/2·(ωx0 + 1/2·ωy0) = ΩVC combination resonances (labelled as R10 and R11, respectively) excited as the bearing clearance δ0 increases beyond 2.0 μm.  As shown in Figures 12 and 13, the equivalent resonant frequency ωy0 quickly falls relative to ωx0.

Global Response Characteristics of VC Resonances around the Bearing Clearance-Free Operations
This results in the resonant amplitudes R2 and R4 shifting to a significantly lower frequency, and the locations of primary resonance R1 and the 1/2-order subharmonic resonance R4 to be closer. For δ0 = 4.0 μm, the equivalent resonant frequencies of the two degrees-of-freedom of System (5) are about  Figure 12 shows the power spectra of the dynamic displacement parameter (DDP) [9]: DDP = xpeak-to-peak + ypeak-to-peak, where xpeak-to-peak = x(t)max − x(t)min, and ypeak-to-peak = y(t)max − y(t)min. To reflect the hysteresis characteristics, here the larger DDP values are selected from the two DDP values obtained as the controlled parameter Ω runs up and down. It is clear that the system contains many types of VC parametric resonances. From Section 4.1, some resonance types of the two degrees-of-freedom of the system have been confirmed, such as the primary resonances (i.e., R1 and R2), 1/2-order subharmonic resonances (i.e., R3 and R4), 2-order superharmonic resonances (i.e., R6 and R7), and the case of ωx0 + ωy0 = ΩVC combination resonance (i.e., R5). Moreover, the cases of 1/2·(ωx0 + ωy0) = ΩVC combination resonance (i.e., the 2-order summed type of combination resonance [39,57]) and 1/2·(1/4·ωx0 + ωy0) = ΩVC combination resonance (labelled as R8 and R9, respectively) were found in our recent work [25]. As shown in Figure 13, these resonant locations can be predicted well by the equivalent resonant frequencies ωx0 and ωy0 from the arithmetic mean of Equation (8) at Ω = 700 rad/s, where we assume that the system is far from any resonance when Ω takes 700 rad/s to achieve an almost linearized dynamic stiffness (7). In addition, there are the other cases of 1/2·(1/2·ωx0 + ωy0) = ΩVC and 1/2·(ωx0 + 1/2·ωy0) = ΩVC combination resonances (labelled as R10 and R11, respectively) excited as the bearing clearance δ0 increases beyond 2.0 μm.  As shown in Figures 12 and 13, the equivalent resonant frequency ωy0 quickly falls relative to ωx0.

Global Response Characteristics of VC Resonances around the Bearing Clearance-Free Operations
This results in the resonant amplitudes R2 and R4 shifting to a significantly lower frequency, and the locations of primary resonance R1 and the 1/2-order subharmonic resonance R4 to be closer. For δ0 = 4.0 μm, the equivalent resonant frequencies of the two degrees-of-freedom of System (5) are about As shown in Figures 12 and 13, the equivalent resonant frequency ω y0 quickly falls relative to ω x0 . This results in the resonant amplitudes R 2 and R 4 shifting to a significantly lower frequency, and the locations of primary resonance R 1 and the 1/2-order subharmonic resonance R 4 to be closer. For δ 0 = 4.0 µm, the equivalent resonant frequencies of the two degrees-of-freedom of System (5) are about ω x0 = 1881.47 rad/s and ω y0 = 1047.21 rad/s, and the amplitudes of R 1 and R 2 as Ω rotates around ω x0/ N b = 235.18 rad/s and ω y0/ N b = 130.90 rad/s in Figure 14 are the two VC primary resonances, respectively. Herein, the primary resonance R 1 and the 1/2-order subharmonic resonance R 4 are close to a slight superposition.
The case of ω x0 + ω y0 = Ω VC combination resonance is triggered from the subcritical and supercritical secondary Hopf bifurcations occurring at points E 5 and E 6 , and bends to the left, exhibiting strong soft spring characteristics (see Figure 15). There are abundant quasi-period responses (e.g., see Figure 16), tori doubling (e.g., see Figure 17) [22,51], and even chaotic motion (e.g., see Figure 18) excitements. Moreover, VC period responses such as period-35 and period-8 motions are also excited by this combination resonance as shown in Figures 19 and 20, respectively. Based on the above analysis, it is clear that the VC parametric combination resonance can result in various response types to the system for a bigger bearing clearance δ 0 . One reason is that, due to the coupling nonlinearities from the bearing clearance and the Hertzian contact of System (5), the amplitude of the combination resonance strongly bends to the left with a large interval of Ω, and this causes the main frequency components p and q of the response to be modulated significantly, because: Secondly, considering the soft spring characteristics of the combination resonance tail, the equivalent resonant frequencies [ω x0 , ω y0 ] in Equation (23) [16][17][18][19][20], respectively. In addition, [ω x0 , ω y0 ] predicted from Equation (23) becomes smaller gradually, except ω y0 at Ω = 355 rad/s (see Table 5). Thirdly, in essence the dynamic stiffness (7) is time-varying and influenced by the response characteristics, so the equivalent resonant frequencies ω x0 and ω y0 might also fluctuate to some extent. Furthermore, due to the dangerous bifurcation [58] shown at points E 7 and E 8 , there are also other complex hysteresis and jumping behaviors in the case of ω x0 + ω y0 = Ω VC combination resonance.
Herein, the primary resonance R1 and the 1/2-order subharmonic resonance R4 are close to a slight superposition.

Discussion on the Verification of Global VC Parametric Resonances
Different from previously works on parametric VC primary resonances, in this paper we extend to study many other parametric VC resonance types such as 1/2-order subharmonic resonances, 2order superharmonic resonances, and various VC combination resonances. For verification of these new findings, we use our proposed method to reanalyze a classical work from Fukata et al. [14], where they presented that sub/super-harmonic, multi-period, quasi-periodic, and chaos-like motions emerge in System (5). As shown in Table 6, this ball bearing model has a larger bearing clearance δ0 = 20 μm and a lighter radical load W = 58.8 N compared to the above-studied model.  Figure 28 presents frequency-response curves of the VC parametric response of the model in Table 6. In this case, the equivalent resonant frequencies of the model are about ωx0 = 1463.75 rad/s

Discussion on the Verification of Global VC Parametric Resonances
Different from previously works on parametric VC primary resonances, in this paper we extend to study many other parametric VC resonance types such as 1/2-order subharmonic resonances, 2order superharmonic resonances, and various VC combination resonances. For verification of these new findings, we use our proposed method to reanalyze a classical work from Fukata et al. [14], where they presented that sub/super-harmonic, multi-period, quasi-periodic, and chaos-like motions emerge in System (5). As shown in Table 6, this ball bearing model has a larger bearing clearance δ0 = 20 μm and a lighter radical load W = 58.8 N compared to the above-studied model.  Figure 28 presents frequency-response curves of the VC parametric response of the model in Table 6. In this case, the equivalent resonant frequencies of the model are about ωx0 = 1463.75 rad/s

Discussion on the Verification of Global VC Parametric Resonances
Different from previously works on parametric VC primary resonances, in this paper we extend to study many other parametric VC resonance types such as 1/2-order subharmonic resonances, 2-order superharmonic resonances, and various VC combination resonances. For verification of these new findings, we use our proposed method to reanalyze a classical work from Fukata et al. [14], where they presented that sub/super-harmonic, multi-period, quasi-periodic, and chaos-like motions emerge in System (5). As shown in Table 6, this ball bearing model has a larger bearing clearance δ 0 = 20 µm and a lighter radical load W = 58.8 N compared to the above-studied model. Table 6. Specifications of the ball bearing in [14].

Item Value
Contact stiffness C b (N/m 3/2 ) 1.334 × 10 10 Ball diameter D b (mm) 11 Figure 28 presents frequency-response curves of the VC parametric response of the model in Table 6. In this case, the equivalent resonant frequencies of the model are about ω x0 = 1463.75 rad/s and ω y0 = 370.74 rad/s. Compared to the case of δ 0 = 4.0 µm in Section 4.2 (i.e., ω x0 = 1881.47 rad/s and ω y0 = 1047.21 rad/s), the resonant excitation frequencies fall a lot in Fukata's model. Herein, the VC primary resonances as Ω rotates around ω x0/ N b = 182.97 rad/s and ω y0/ N b = 46.34 rad/s are also excited (i.e., the R 1 and R 3 ranges). These two primary resonances exhibit soft spring and soft-hard spring characteristics similar to the case of δ 0 = 4.0 µm studied. Through supercritical and subcritical period-doubling bifurcations trigger at points F 1 and F 2 , the 1/2-order subharmonic resonance in the horizontal direction is strong excited in a wide range R 4 as Ω rotates around 2·ω y0/ N b = 92.68 rad/s, and VC period-2 resonant responses occur. It should be note that 2-order superharmonic resonance in the vertical direction dominates the resonant amplitude R 6 in the vertical direction as Ω rotates around 1/2·ω x0/ N b = 91.48 rad/s, because the main frequency component of x(τ) in this range is 2-order harmonics (i.e., the response characteristics at Ω = 88 rad/s presented in Figure 29). Moreover, the period-2 solution branch F 1 -F 2 undergoes period-doubling bifurcations at F 3 , F 4 , F 5 and F 6 , generating period-4 branch F 3 -F 4 and F 5 -F 6 . Especially, there are various VC motion types arising subsequently on branch F 5 -F 6 , including chaotic responses, and this branch contains the beat/chaos-like range b-b described by Fukata et al. [14] (see Figure 30).  [14] (see Figure 30).  [14] (see Figure 30).  Figure 30. The peak-to-peak frequency-response curves of literature [14], in (a) x and (b) y directions. Herein, the blue and red lines are the period-1 and period-2 solution branches of Figure 28, respectively, obtained by our work.
Furthermore, the stable period-1 solution branch undergoes supercritical Hopf bifurcations at points F 7 and F 8 , and the coming unstable solution branch F 7 -F 8 (i.e., R 5 range around ω x0 /N b + ω y0/ N b = 229.31 rad/s) contains various VC responses such as period-5, period-6 and beat/chaos-like motions (see Figure 30). The response characteristics of the quasi-period motion as Ω rotates at 242.2277 (i.e., in literature [14], ω s takes 6000 rpm) are presented in Figure 31. Herein, the values of the main incommensurable dimensionless VC frequencies p and q are about 0.7715 and 0.2285, respectively, and p + q = 1, where the value 1 is the dimensionless VC excitation frequency Ω VC = N b ·Ω. Moreover, for the original time scale, the corresponding main frequency components of the response are N b ·Ω ·p ≈ 1495.03 rad/s and N b ·Ω ·q ≈ 442.79 rad/s, which are close to the predicted equivalent resonant frequencies ω x0 = 1463.75 rad/s and ω y0 = 370.74 rad/s of the two degrees-of-freedom of System (5). As stated in Section 4.1, this indicates that the case of ω x0 + ω y0 = Ω VC combination resonance occurs in range R 5 . Different from the previous analysis (see Figures 7b and 16b), as shown in Figure 31b, in this case the response-frequency spectrum contains component of 2·q, so there may be the case of ω x0 + 1/2·ω y0 = Ω VC combination resonance also excited in R 5 , and actually the value ω x0 /N b + 1/2·ω y0/ N b = 206.14 rad/s (see Figure 28) is in range R 5 . This can explain why the instability range R 5 excited by the combination resonance is so large in Fukata's model. In other words, due to the more asymmetric stiffness characteristics (where ω x0 = 1463.75 rad/s and ω y0 = 370.74 rad/s), Fukata's model can contain more types of parametric resonances, and these resonances may overlap (literature [27] contains a detailed discussion of this overlapping phenomenon). Due to the vibration mechanism being more complicated, many researchers [44,60] have devoted themselves to investigating the dynamic characteristics of asymmetric rolling bearing-rotor systems. The response characteristics of the period-5 motion as Ω rotates at 211.9492 rad/s (i.e., in literature [14], ω s takes 5250 rpm) are presented in Figure 32, where its main frequency components p and q are 4/5 and 1/5, respectively; satisfying p + q = 1, and using the above analysis process, it can be found that this period-5 response is also excited by the case of ω x0 + ω y0 = Ω VC combination resonance. Therefore, the case of ω x0 + ω y0 = Ω VC combination resonance is the internal mechanism of excitation of various types of VC responses in R 5 range described by Fukata et al. [14] (see Figure 30). In addition, it can be easily verified by our analysis process that the resonant amplitude R 12 comes from 6-order superharmonic resonance in the vertical direction, because this resonant range emerge as Ω rotates around 1/6·ω x0 /N b = 30.50 rad/s (see Figures 28 and 30). As shown in Figure 33, our above analysis results agree well with the numerical bifurcation diagrams. Figure 30. The peak-to-peak frequency-response curves of literature [14], in (a) x and (b) y directions. Herein, the blue and red lines are the period-1 and period-2 solution branches of Figure 28, respectively, obtained by our work. Furthermore, the stable period-1 solution branch undergoes supercritical Hopf bifurcations at points F7 and F8, and the coming unstable solution branch F7-F8 (i.e., R5 range around ωx0/Nb + ωy0/Nb = 229.31 rad/s) contains various VC responses such as period-5, period-6 and beat/chaos-like motions (see Figure 30). The response characteristics of the quasi-period motion as Ω rotates at 242.2277 (i.e., in literature [14], ωs takes 6000 rpm) are presented in Figure 31. Herein, the values of the main incommensurable dimensionless VC frequencies p and q are about 0.7715 and 0.2285, respectively, and p + q = 1, where the value 1 is the dimensionless VC excitation frequency ΩVC = Nb·Ω. Moreover, for the original time scale, the corresponding main frequency components of the response are Nb·Ω ·p ≈ 1495.03 rad/s and Nb·Ω ·q ≈ 442.79 rad/s, which are close to the predicted equivalent resonant frequencies ωx0 = 1463.75 rad/s and ωy0 = 370.74 rad/s of the two degrees-of-freedom of System (5). As stated in Section 4.1, this indicates that the case of ωx0 + ωy0 = ΩVC combination resonance occurs in range R5. Different from the previous analysis (see Figures 7b and 16b), as shown in Figure 31b, in this case the response-frequency spectrum contains component of 2·q, so there may be the case of ωx0 + 1/2·ωy0 = ΩVC combination resonance also excited in R5, and actually the value ωx0/Nb + 1/2·ωy0/Nb = 206.14 rad/s (see Figure 28) is in range R5. This can explain why the instability range R5 excited by the combination resonance is so large in Fukata's model. In other words, due to the more asymmetric stiffness characteristics (where ωx0 = 1463.75 rad/s and ωy0 = 370.74 rad/s), Fukata's model can contain more types of parametric resonances, and these resonances may overlap (literature [27] contains a detailed discussion of this overlapping phenomenon). Due to the vibration mechanism being more complicated, many researchers [44,60] have devoted themselves to investigating the dynamic characteristics of asymmetric rolling bearing-rotor systems. The response characteristics of the period-5 motion as Ω rotates at 211.9492 rad/s (i.e., in literature [14], ωs takes 5250 rpm) are presented in Figure 32, where its main frequency components p and q are 4/5 and 1/5, respectively; satisfying p + q = 1, and using the above analysis process, it can be found that this period-5 response is also excited by the case of ωx0 + ωy0 = ΩVC combination resonance. Therefore, the case of ωx0 + ωy0 = ΩVC combination resonance is the internal mechanism of excitation of various types of VC responses in R5 range described by Fukata et al. [14] (see Figure 30). In addition, it can be easily verified by our analysis process that the resonant amplitude R12 comes from 6-order superharmonic resonance in the vertical direction, because this resonant range emerge as Ω rotates around 1/6·ωx0/Nb = 30.50 rad/s (see Figures   28 and 30). As shown in Figure 33, our above analysis results agree well with the numerical bifurcation diagrams.   For experimental verification, because the VC vibration response is a weak signal, the vibration is usually measured by non-contact displacement sensor in experiments [24,34,61]. The hysteresis characteristics of primary contact resonances have been validated by our previous work [24], and other parametric VC resonances types shown in this paper may be verified in the similar experimental process.

Conclusions
In this paper, we investigated the response characteristics of the global VC parametric resonances in a ball bearing system around the bearing clearance-free operations, considering the nonlinearities of bearing clearance and the Hertzian contact between the rolling elements and raceway. The VC period-M motions and their stabilities were traced using the HB-AFT method and Floquet theory. It was shown that the system can contain primary resonances, 1/2-order subharmonic resonances, 2-order superharmonic resonances, and various cases of combination resonances. The 1/2-order VC subharmonic resonances are excited by supercritical/subcritical period-doubling bifurcations, leading to the emergence of VC period-2 motions. Moreover, it was found that the cases of ωx0 + ωy0 = ΩVC, 1/2·(ωx0 + ωy0) = ΩVC, 1/2·(1/4·ωx0 + ωy0) = ΩVC, 1/2·(1/2·ωx0 + ωy0) = ΩVC, and 1/2·(ωx0 + 1/2·ωy0) = ΩVC combination resonances can be excited by supercritical/subcritical secondary Hopf bifurcations. Due to the frequency modulation from hysteresis and dynamic stiffness characteristics, the VC parametric combination resonance can result in various VC motion types, such as period-4, period-8, period-35, different cases of quasi-period, and even chaotic responses in the system. Furthermore, it was indicated that the parametric combination resonance is an internal physical  For experimental verification, because the VC vibration response is a weak signal, the vibration is usually measured by non-contact displacement sensor in experiments [24,34,61]. The hysteresis characteristics of primary contact resonances have been validated by our previous work [24], and other parametric VC resonances types shown in this paper may be verified in the similar experimental process.

Conclusions
In this paper, we investigated the response characteristics of the global VC parametric resonances in a ball bearing system around the bearing clearance-free operations, considering the nonlinearities of bearing clearance and the Hertzian contact between the rolling elements and raceway. The VC period-M motions and their stabilities were traced using the HB-AFT method and Floquet theory. It was shown that the system can contain primary resonances, 1/2-order subharmonic resonances, 2-order superharmonic resonances, and various cases of combination resonances. The 1/2-order VC subharmonic resonances are excited by supercritical/subcritical period-doubling bifurcations, leading to the emergence of VC period-2 motions. Moreover, it was found that the cases of ωx0 + ωy0 = ΩVC, 1/2·(ωx0 + ωy0) = ΩVC, 1/2·(1/4·ωx0 + ωy0) = ΩVC, 1/2·(1/2·ωx0 + ωy0) = ΩVC, and 1/2·(ωx0 + 1/2·ωy0) = ΩVC combination resonances can be excited by supercritical/subcritical secondary Hopf bifurcations. Due to the frequency modulation from hysteresis and dynamic stiffness characteristics, the VC parametric combination resonance can result in various VC motion types, such as period-4, period-8, period-35, different cases of quasi-period, and even chaotic responses in the system. Furthermore, it was indicated that the parametric combination resonance is an internal physical For experimental verification, because the VC vibration response is a weak signal, the vibration is usually measured by non-contact displacement sensor in experiments [24,34,61]. The hysteresis characteristics of primary contact resonances have been validated by our previous work [24], and other parametric VC resonances types shown in this paper may be verified in the similar experimental process.

Conclusions
In this paper, we investigated the response characteristics of the global VC parametric resonances in a ball bearing system around the bearing clearance-free operations, considering the nonlinearities of bearing clearance and the Hertzian contact between the rolling elements and raceway. The VC period-M motions and their stabilities were traced using the HB-AFT method and Floquet theory. It was shown that the system can contain primary resonances, 1/2-order subharmonic resonances, 2-order superharmonic resonances, and various cases of combination resonances. The 1/2-order VC subharmonic resonances are excited by supercritical/subcritical period-doubling bifurcations, leading to the emergence of VC period-2 motions. Moreover, it was found that the cases of ω x0 + ω y0 = Ω VC , 1/2·(ω x0 + ω y0 ) = Ω VC , 1/2·(1/4·ω x0 + ω y0 ) = Ω VC , 1/2·(1/2·ω x0 + ω y0 ) = Ω VC , and 1/2·(ω x0 + 1/2·ω y0 ) = Ω VC combination resonances can be excited by supercritical/subcritical secondary Hopf bifurcations. Due to the frequency modulation from hysteresis and dynamic stiffness characteristics, the VC parametric combination resonance can result in various VC motion types, such as period-4, period-8, period-35, different cases of quasi-period, and even chaotic responses in the system. Furthermore, it was indicated that the parametric combination resonance is an internal physical mechanism of VC quasi-period responses trigged in the rolling bearing system. The above results indicate that adjusting the bearing clearance can control the resonant amplitude and its location. Finally, the classical study from Fukata et al. [14] was reanalyzed by our proposed method and the experimental verification of our findings was also discussed. This paper helps clarify the mechanism of complex nonlinear vibration behaviors in rolling bearings.