Investigation on Modeling and Formation Mechanism of Dynamic Rotational Error for Spindle-Rolling Bearing System

: In the ﬁeld of precision machining, the spindle-rolling bearing (SRB) system is widely used on the machine tool as one of the most fundamental and important components. The rotational error motions of the SRB system have signiﬁcant e ﬀ ects on the machining accuracy (contour accuracy and surface roughness). Over the past decades, much work has been focused on the measurement of spindle balancing and rotational error motions, the vibrations response induced by the nonlinear sti ﬀ ness and surface waviness of the bearing. However, the formative mechanism of the rotational error motions for the SRB system is not well understood. In this paper, the dynamic model of the SRB system considering the bearing nonlinearity is established. Seeking to reveal the e ﬀ ects of surface waviness of the bearing raceway, unbalance mass and disturbance force on the dynamic rotational error, the modeling method and formative mechanism of the dynamic rotational error for the SRB system is explored both theoretically and experimentally. Then, numerical simulation is performed to analyze the inﬂuence of the bearing raceway waviness, unbalance mass and disturbance force on the dynamic rotational error. An experimental setup is established based on a typical SRB system and a series of experiments are carried out. The experimental results are in good agreement with the theoretical and simulation results, which can demonstrate the feasibility and validity of the modeling method. Furthermore, this method can be e ﬀ ectively applied to the design and development phases of an SRB system to improve dynamic rotational accuracy.


Introduction
The spindle-rolling bearing (SRB) system is a key component of the machine tool or other rotating machines, which plays a key role in working. The rotational accuracy is a critical technical indicator for the SRB system when evaluating the dynamic performance [1,2]. With the increasing demands on machining accuracy of the machine tool, much effort has been made to research the measuring methods on the rotational accuracy of the SRB system. Metrologists began to measure rotational accuracy as early as the 1900s. Scheslinger [3] used a mechanical indicator to measure the total indicator reading (TIR), which is recognized to be a static measuring method. Tlusty [4] and Bryan [5] pioneered modern measuring methods by using more accurate data acquisition (DAQ) system and evaluating methods. Their methods lay the foundation for ANSI(American National Standards setup is designed to verify the proposed model. A series of measuring experiments for the dynamic rotational error corresponding to different influenced factors are performed. The conclusions and prospects of this study are given in Section 5.

Modeling of Dynamic Rotational Error for Spindle-Rolling Bearing System
Angular contact ball bearings are widely used in the SRB system. The nonlinear and time-varying characteristics of the bearing, under the influence of geometric error, bearing structure and dynamic load, is the root cause of the nonlinearity and has great important effects on the dynamic rotational error of the SRB system. In this section, the dynamic model of the SRB system considering the nonlinear and time-varying behaviors of the bearing are established, based on which the dynamic rotational error is evaluated.

Nonlinear Bearing Model
The dynamic supporting stiffness and time-varying bearing force are important parameters in modeling the nonlinear behaviors of rolling bearings. During operation, the dynamic supporting stiffness of the SRB system is nonlinear according to the Hertz contact theory (the Hertz contact theory illustrates the distribution rule of local stress and strain generated by the contact of two objects under compression), while the bearing force is time-variable due to the machining error of the bearing rings and rolling elements.

Modeling of Dynamic Supporting Stiffness
The geometry of an angular contact ball bearing and the coordinate systems is shown in Figure 1. The global coordinate system (O-XYZ) is fixed at the rotating axis center of the outer ring. The inner ring rotates at a constant speed ω 0 with Z axis while the outer ring is fixed. α o is the initial contact angle of the bearing. The outer ring is fixed and assume that δ x , δ y , δ z , θ x and θ y are the three translational and two angular displacements of the inner ring. The displacement vector of the inner ring is expressed as Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 24 the dynamic rotational error corresponding to different influenced factors are performed. The conclusions and prospects of this study are given in Section 5.

Modeling of Dynamic Rotational Error for Spindle-Rolling Bearing System
Angular contact ball bearings are widely used in the SRB system. The nonlinear and timevarying characteristics of the bearing, under the influence of geometric error, bearing structure and dynamic load, is the root cause of the nonlinearity and has great important effects on the dynamic rotational error of the SRB system. In this section, the dynamic model of the SRB system considering the nonlinear and time-varying behaviors of the bearing are established, based on which the dynamic rotational error is evaluated.

Nonlinear Bearing Model
The dynamic supporting stiffness and time-varying bearing force are important parameters in modeling the nonlinear behaviors of rolling bearings. During operation, the dynamic supporting stiffness of the SRB system is nonlinear according to the Hertz contact theory (the Hertz contact theory illustrates the distribution rule of local stress and strain generated by the contact of two objects under compression), while the bearing force is time-variable due to the machining error of the bearing rings and rolling elements.

Modeling of Dynamic Supporting Stiffness
The geometry of an angular contact ball bearing and the coordinate systems is shown in Figure  1. The global coordinate system (O-XYZ) is fixed at the rotating axis center of the outer ring. The inner ring rotates at a constant speed ω0 with Z axis while the outer ring is fixed. αo is the initial contact angle of the bearing. The outer ring is fixed and assume that δx, δy, δz, θx and θy are the three translational and two angular displacements of the inner ring. The displacement vector of the inner ring is expressed as There are Nb rolling elements in the bearing and the position angle of the kth rolling element, ϕbk, could be obtained from the following formula: The local coordinate system (o-xyz), whose center is at the mass center of the rolling element, rotates at an angular speed ωbk with the Z axis. ωbk is the orbital angular speed of the rolling element and can be expressed as There are N b rolling elements in the bearing and the position angle of the kth rolling element, φ bk , could be obtained from the following formula: Appl. Sci. 2020, 10, 5753 4 of 23 The local coordinate system (o-xyz), whose center is at the mass center of the rolling element, rotates at an angular speed ω bk with the Z axis. ω bk is the orbital angular speed of the rolling element and can be expressed as where d m is the pitch diameter of the bearing, D is the diameter of the rolling element. The rolling elements also rotate at an angular speed ω sk around the axis of itself, which is called the spinning axis. The self-rotation angular speed ω sk can be easily obtained by [31].
The angle between the spinning axis of the rolling elements and Z axis is β k , which is called as rolling element attitude angle. β k can be expressed as [31].
When the rotational speed ω 0 and the axial preload are zero, both the inner contact angle and the outer contact angle are the same and equal to the initial contact angle α o . The inner raceway groove curvature center, outer raceway groove curvature center and the rolling element center are collinear (red line), as shown in Figure 2. The distance between the inner raceway groove curvature center and the rolling element center, as well as the distance between the outer raceway groove curvature center and the rolling element center before deformation are where r i and r o are the inner and outer raceway groove curvature radius, respectively, f i and f o are the inner and outer raceway groove curvature radius coefficient, respectively.
where dm is the pitch diameter of the bearing, D is the diameter of the rolling element. The rolling elements also rotate at an angular speed ωsk around the axis of itself, which is called the spinning axis. The self-rotation angular speed ωsk can be easily obtained by [31].
The angle between the spinning axis of the rolling elements and Z axis is βk, which is called as rolling element attitude angle. βk can be expressed as [31].
When the rotational speed ω0 and the axial preload are zero, both the inner contact angle and the outer contact angle are the same and equal to the initial contact angle α o . The inner raceway groove curvature center, outer raceway groove curvature center and the rolling element center are collinear (red line), as shown in Figure 2. The distance between the inner raceway groove curvature center and the rolling element center, as well as the distance between the outer raceway groove curvature center and the rolling element center before deformation are where ri and ro are the inner and outer raceway groove curvature radius, respectively, fi and fo are the inner and outer raceway groove curvature radius coefficient, respectively. With the rotational speed ω0 and axial preload increase, the inner raceway groove curvature center, outer raceway groove curvature center and the rolling element center are not collinear (blue line). The inner contact angle and the outer contact angle are not equal to each other. The inner contact angle increase while the outer contact angle decreases due to the centrifugal force. Besides, the distance between the inner/outer raceway groove curvature center and rolling element center changes in the Y-Z plane. The distance between the inner raceway groove curvature center and the With the rotational speed ω 0 and axial preload increase, the inner raceway groove curvature center, outer raceway groove curvature center and the rolling element center are not collinear (blue line). The inner contact angle and the outer contact angle are not equal to each other. The inner contact angle increase while the outer contact angle decreases due to the centrifugal force. Besides, the distance between the inner/outer raceway groove curvature center and rolling element center changes in the Y-Z plane. The distance between the inner raceway groove curvature center and the rolling element center, as well as the distance between the outer raceway groove curvature center and the rolling element center, after deformation are where δ ik /δ ok is the contact deformation between the inner/outer bearing ring and the kth rolling element.
In accordance with the displacement coordinate relations, the axial distance between the inner and outer raceway groove curvature centers is where η i is the distance between the rotating axis Z and inner raceway groove curvature center and can be calculated by It convenient to introduce new variables X 1k and X 2k , as shown in Figure 2. It can be seen from Figure 2 that at any rolling element location Using the Pythagorean Theorem, it can be seen from Figure 2 that [31].
As shown in Figure 3, the load equilibrium equations for the kth rolling element can be expressed as [31].
where Q ik /Q ok is the contact force between the inner/outer ring and the rolling element. Q ik and Q ok can be expressed as Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 24 λik / λok represent the friction coefficients between the inner/outer raceway and rolling element. As for the outer raceway control, λik = 0 and λok = 2 The gyroscopic moment Mgk and centrifugal force Fck acting on the kth rolling element can be described as follows: where mb, Ib is the mass and mass moment inertia of the kth rolling element, respectively. Equations (11) and (12) are highly coupled and nonlinear and can be resolved by Newton-Raphson method. Then the contact forces Qi:, Qo:, contact angles αi:, αo:, and gyroscopic moments Mi: Mo:, between each rolling element and inner/outer rings can be obtained. The resultant force vector of the inner ring due to all of the contacting rolling elements is According to the relationship between the resultant forces and displacements of the inner ring, the dynamic supporting stiffness matrix of the bearing can be obtained λ ik /λ ok represent the friction coefficients between the inner/outer raceway and rolling element. As for the outer raceway control, λ ik = 0 and λ ok = 2 The gyroscopic moment M gk and centrifugal force F ck acting on the kth rolling element can be described as follows: (14) where m b , I b is the mass and mass moment inertia of the kth rolling element, respectively.
Equations (11) and (12) are highly coupled and nonlinear and can be resolved by Newton-Raphson method. Then the contact forces Q i: , Q o: , contact angles α i: , α o: , and gyroscopic moments M i: M o: , between each rolling element and inner/outer rings can be obtained. The resultant force vector of the inner ring due to all of the contacting rolling elements is According to the relationship between the resultant forces and displacements of the inner ring, the dynamic supporting stiffness matrix of the bearing can be obtained The machining errors of the bearing parts, which has great important effects on the bearing force and dynamic rotational error, are mainly include shape errors, surface waviness and roughness. The surface roughness is not considered in modeling the bearing force for the wavelength of the surface roughness is smaller than the width of the Hertz contact [31].
The outer ring is fixed and the displacements for the inner raceway groove curvature center corresponding to the kth rolling element, including two translations along the y-, z-axes, and one rotation with respect to the x-axis, are expressed as δ ik = u yk , u zk , θ xk T (18) Under the assumption of small deflection, δ ik could be obtained through coordinate transformations where T k is the transformation matrix, and can be expressed as ∆w k denotes the additional displacements caused by the waviness of the bearing rings and rolling elements. Similarly, ∆D k represents the additional displacements due to roundness error of the bearing rings and rolling elements.
As illustrated in Figure 4, the waviness exists on each bearing components. p ik , p ok and q ik , q ok represent radial and axial waviness of inner and outer raceways, w ik and w ok denote the waviness of rolling elements. The harmonic functions are utilized to express the surface waviness due to the waviness vary with time periodically.
A out cos n out ω bk t + C ba cos(n ba ω sk t + ϕ ba ) C ba cos(n ba ω sk t + π + ϕ ba ) (21) where A in , A out , B in , B in , C ba , C ba are the amplitudes of the surface waviness; n in , n out , n ba , are the harmonic orders of the surface waviness; φ in , ψ in , φ out , ψ out , ϕ ba are the initial phases angle of the surface waviness.
Appl. Sci. 2020, 10, 5753 where Ain, Aout, Bin, Bin, Cba, Cba are the amplitudes of the surface waviness; nin, nout, nba, are the harmonic orders of the surface waviness; ϕin, ψin, ϕout, ψout, φba are the initial phases angle of the surface waviness. In accordance with the geometric relations, one can obtain The effects of the surface waviness on the dynamic supporting stiffness are very small; however, the effects of the surface waviness on the dynamic rotational error can be reflected in the form of time-varying bearing force. The distances between the inner/outer raceway groove curvature center and the rolling element center, considering the effect of surface waviness, are The geometric relationships can be written as follows: In accordance with the geometric relations, one can obtain The effects of the surface waviness on the dynamic supporting stiffness are very small; however, the effects of the surface waviness on the dynamic rotational error can be reflected in the form of time-varying bearing force. The distances between the inner/outer raceway groove curvature center and the rolling element center, considering the effect of surface waviness, are The geometric relationships can be written as follows: where v yk and v zk are the displacement variation of the k th rolling element along the y-axis, zaxis, respectively. The time-varying forces from the kth rolling element to the inner ring considering the surface waviness can be stated as: According to Equations (11) and (12), the time-varying force vector Q ik can be expressed as follows: Appl. Sci. 2020, 10, 5753 9 of 23 Then the sum of the time-varying forces from all of the contacting rolling elements can be obtained: Figure 5 shows the layout of the SRB system. The spindle is supported by two pairs of angular contact bearings. The front and rear bearings are arranged back to back. The types of the front and rear bearing are 7017 and 7016 C, respectively. The light preload is exerted on the SRB system in assembling process. The system coordinate O-U s V s W s is fixed at the mass center of the spindle. Table 1 shows the parameters of the SRB system. In this section, the dynamic model of the SRB system considering the dynamic supporting stiffness, time-varying bearing force caused by the surface waviness, dynamic unbalance force and the disturbance force, will be derived, which lays a foundation for evaluating the dynamic rotational error.  The simplified dynamic model of the SRB system is shown in Figure 6. The spindle can be regarded as a rigid body and has five degrees of freedom, i.e., three translational freedom along the Us-, Vs-and Ws-axes, and two rotational freedom with respect to the Us-and Vs-axes. The dynamic displacement vector of the rotor is Figure 5. Layout of the spindle-rolling bearing (SRB) system. Table 1. Calculation parameters of the SRB system.

Items Parameters Items Parameters
The simplified dynamic model of the SRB system is shown in Figure 6. The spindle can be regarded as a rigid body and has five degrees of freedom, i.e., three translational freedom along the U s -, V s -and Ws-axes, and two rotational freedom with respect to the U s -and V s -axes. The dynamic displacement vector of the rotor is q = [u s , v s , w s , θ u s , θ v s ] T (28) regarded as a rigid body and has five degrees of freedom, i.e., three translational freedom along the Us-, Vs-and Ws-axes, and two rotational freedom with respect to the Us-and Vs-axes. The dynamic displacement vector of the rotor is where [M] = [ms, ms, ms, Is, Is] T , is the mass matrix of the system, ms is the mass of the spindle and Is is the inertia moment of the spindle about the Us-axis and Vs-axis. K b is the stiffness matrix which can be calculated in Section 2. C is the structural damping matrix constructed from modal damping ratios identified experimentally. G is the gyroscopic matrix and can be expressed as ..
where [M] = [m s , m s , m s , I s , I s ] T , is the mass matrix of the system, m s is the mass of the spindle and I s is the inertia moment of the spindle about the U s -axis and V s -axis. K b is the stiffness matrix which can be calculated in Section 2. C is the structural damping matrix constructed from modal damping ratios identified experimentally. G is the gyroscopic matrix and can be expressed as J s is the inertia moment of the rotor about the W s -axis. F e F dr and F b are the dynamic unbalance force caused by unbalanced mass, disturbance force and the time-varying bearing force caused by surface waviness, respectively. F e can be states as {F e } = m u e s ω 2 0 [cos(ω 0 t), sin(ω 0 t), 0, l 3 cos(ω 0 t), l 3 sin(ω 0 t)] T (31) where m u is the eccentric mass and e s is the eccentric distance. The cutting forces have not been considered, since they involve more complicated issues not covered in this paper. The effect of cutting force on the dynamic rotational error will be studied in the future.

Evaluation of Dynamic Rotational Error
It is difficult to measure the dynamic rotational error of the SRB system directly. Precision artifacts are adopted to aid the metrology to identify the dynamic features of the SRB system. Figure 7 shows the diagram of radial error motion measurement and discrete time samples of displacement data measured in radial direction.

Evaluation of Dynamic Rotational Error
It is difficult to measure the dynamic rotational error of the SRB system directly. Precision artifacts are adopted to aid the metrology to identify the dynamic features of the SRB system. Figure  7 shows the diagram of radial error motion measurement and discrete time samples of displacement data measured in radial direction. The displacement data express sinusoidal trends, with a deviation from zero level. In accordance with the ISO or ANSI/ASME test standard, different kinds of error motion values can be evaluated by using the displacement signals of the artifact [6,7]. The least-square method is used for evaluating the synchronous error motion value. The solution of the synchronous error motion value is transformed into a nonlinear least squares problem. The distances xLSC and yLSC of the least squares center from the polar chart center can be calculated as follows: where xi and yi are the values of the ith data points in the X and Y direction, respectively.
RLSC is the radius the least squares circle and can be calculated as follows:  The displacement data express sinusoidal trends, with a deviation from zero level. In accordance with the ISO or ANSI/ASME test standard, different kinds of error motion values can be evaluated by using the displacement signals of the artifact [6,7]. The least-square method is used for evaluating the synchronous error motion value. The solution of the synchronous error motion value is transformed into a nonlinear least squares problem. The distances x LSC and y LSC of the least squares center from the polar chart center can be calculated as follows: where x i and y i are the values of the ith data points in the X and Y direction, respectively. R LSC is the radius the least squares circle and can be calculated as follows: Error motion polar plots and corresponding error motion values are depicted in the Figure 7.

Numerical Simulation
Based on the dynamic model established in Section 2 and the ISO 230-7 or ASME B89.3.4 test standard, the influence of different error factors on the dynamic rotational error can be analyzed. Polar plots of the dynamic rotational error can be obtained by using the spindle displacement response, which are measured by the displacement sensor at the end of the SRB system. The least-square method can be used for evaluating the value of synchronous error motion. In this section, the dynamic rotational error of the SRB system has been analyzed quantitatively.

Influence of Bearing Raceway Waviness on Dynamic Rotational Error
Time-varying restoring bearing force can be generated due to the surface waviness of the bearing parts [25]. The excitation caused by bearing waviness is applied on the SRB system. The spindle displacement response is represented as dynamic rotational error when the spindle rotates at high speed. In the simulation, only the effect of a rear bearing waviness is taken into account. Assuming that the bearing inner ring raceway is ideally circular, the rolling elements do not have any shape or dimensional errors. The outer raceway surface has a single order of radial and axial sinusoidal waviness. The amplitude of the bearing waviness is 1-2 µm, and the harmonic order is 11-20 UPR (Undulations Per Revolution). To show the relationships more distinctly between radial errors and the number of undulations and rolling elements, values that are larger than actual errors are selected. The number of rear bearing rolling elements is 20. The length of numerical simulation is set to 20 revolutions. The discrete time response signals and the polar plots of the radial error motion are shown in Figure 8. speed. In the simulation, only the effect of a rear bearing waviness is taken into account. Assuming that the bearing inner ring raceway is ideally circular, the rolling elements do not have any shape or dimensional errors. The outer raceway surface has a single order of radial and axial sinusoidal waviness. The amplitude of the bearing waviness is 1-2 μm, and the harmonic order is 11-20 UPR (Undulations Per Revolution). To show the relationships more distinctly between radial errors and the number of undulations and rolling elements, values that are larger than actual errors are selected. The number of rear bearing rolling elements is 20. The length of numerical simulation is set to 20 revolutions. The discrete time response signals and the polar plots of the radial error motion are shown in Figure 8.  As can be seen from Figures 8 and 9, when the bearing waviness has the same amplitude, the radial displacement response in the V s -axis is related to the order of the bearing waviness. When the order of the waviness and the number of rolling elements are the same, the waviness has almost no effect on the rotational error motion. However, when the difference between the order and the number of rolling elements is 1, a very obvious asynchronous error motion can be generated in the SRB system. Meanwhile, waviness hardly causes synchronous error motion. The same trends can be obtained by comparing the influence of the waviness on the rotational error of a single bearing in the literatures [32,33]. In the comparing literature, an angular contact ball bearing (SKF, 7012/CD) is chosen for the calculation. The comparison of the effects on radial error motion is shown in Figure 10.
When the order of surface waviness is the same, the magnitude of the bearing waviness determines the magnitude of the rotational error. The rotational error motion value, especially for asynchronous error, increases with the increase of the waviness amplitude. effect on the rotational error motion. However, when the difference between the order and the number of rolling elements is 1, a very obvious asynchronous error motion can be generated in the SRB system. Meanwhile, waviness hardly causes synchronous error motion. The same trends can be obtained by comparing the influence of the waviness on the rotational error of a single bearing in the literatures [32,33]. In the comparing literature, an angular contact ball bearing (SKF, 7012/CD) is chosen for the calculation. The comparison of the effects on radial error motion is shown in Figure 10.   number of rolling elements is 1, a very obvious asynchronous error motion can be generated in the SRB system. Meanwhile, waviness hardly causes synchronous error motion. The same trends can be obtained by comparing the influence of the waviness on the rotational error of a single bearing in the literatures [32,33]. In the comparing literature, an angular contact ball bearing (SKF, 7012/CD) is chosen for the calculation. The comparison of the effects on radial error motion is shown in Figure 10.

Influence of Dynamic Unbalance on Dynamic Rotational Error
As a typical rotary machine, dynamic unbalance force is a critical factor attributed to the error motions of the SRB system when its center of mass deviates from the axis of rotation. As a kind of dynamic excitation, the frequency of the dynamic unbalance force is the same as the rotation frequency of the spindle. In the traditional test of the spindle rotational error motion, the fundamental frequency components, including installation eccentricity error of the target, are filtered. The dynamic unbalance force is a function of the rotational speed. From dynamic point of view, the vibration caused by dynamic unbalance force directly affects the total indicator reading (TIR) of the SRB system. Assume that the spindle has an eccentric mass m u = 0.005 kg on the pulley side, and the eccentricity e c = 0.054 m. The spindle rotational speed ranges from 1000 to 10,000 rpm. The response displacements at the front end of the SRB system are shown in Figure 11.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 24 Figure 11. Effects of dynamic unbalanced force on total indicator reading (TIR) at different speeds.
The trends of the overall effect of spindle unbalance on TIR, synchronous and asynchronous errors are shown in Figure 12. As can be seen in Figure 12, the TIR caused by the spindle unbalance increases from 0.07 μm in 1000 rpm to 5.61 μm in 10,000 rpm. Spindle unbalance has little effect on the synchronous error, and the overall synchronous error is less than 0.4 μm. The asynchronous error increases from 0.03 to 2.27 μm, as the speed increasing from 1000 to 10,000 rpm. The trends of the overall effect of spindle unbalance on TIR, synchronous and asynchronous errors are shown in Figure 12. As can be seen in Figure 12, the TIR caused by the spindle unbalance increases from 0.07 µm in 1000 rpm to 5.61 µm in 10,000 rpm. Spindle unbalance has little effect on the synchronous error, and the overall synchronous error is less than 0.4 µm. The asynchronous error increases from 0.03 to 2.27 µm, as the speed increasing from 1000 to 10,000 rpm. The frequency domain analysis of the simulated displacement response at 10,000 rpm is shown in Figure 13. The frequency domain analysis shows that it mainly includes two peaks. The frequency of the first peak corresponds to the excitation frequency of the spindle dynamic unbalance force, that is, the fundamental frequency (spindle rotational frequency). The frequency of the second peak corresponds to the spindle resonance frequency (natural frequency of the SRB system). The asynchronous error motion presents large value due to the resonance frequency is non-integer multiples of the fundamental frequency.

Influence of Disturbance Force on Dynamic Rotational Error
The spindle is rotated by the drive system. The mechanical spindle often uses a wedge belt or a timing belt drive, while the electric spindle is driven by built-in motor. During the rotation of the spindle, the transmission system and the electromagnetic force will generate a disturbance force on the SRB system. The dynamic response of the SRB system to disturbance forces is also reflected in the dynamic rotational error. Based on the established model, the influence of disturbance force on dynamic rotational error can be analyzed numerically. Assuming that the frequency of the disturbance force is an integer multiple of the fundamental frequency, the frequency order ranges from 2 to 10 UPR. Polar plots of the dynamic rotational error, which are caused by the disturbance force with the amplitude of 200 N, are shown in Figure 14. The frequency domain analysis of the simulated displacement response at 10,000 rpm is shown in Figure 13. The frequency domain analysis of the simulated displacement response at 10,000 rpm is shown in Figure 13. The frequency domain analysis shows that it mainly includes two peaks. The frequency of the first peak corresponds to the excitation frequency of the spindle dynamic unbalance force, that is, the fundamental frequency (spindle rotational frequency). The frequency of the second peak corresponds to the spindle resonance frequency (natural frequency of the SRB system). The asynchronous error motion presents large value due to the resonance frequency is non-integer multiples of the fundamental frequency.

Influence of Disturbance Force on Dynamic Rotational Error
The spindle is rotated by the drive system. The mechanical spindle often uses a wedge belt or a timing belt drive, while the electric spindle is driven by built-in motor. During the rotation of the spindle, the transmission system and the electromagnetic force will generate a disturbance force on the SRB system. The dynamic response of the SRB system to disturbance forces is also reflected in the dynamic rotational error. Based on the established model, the influence of disturbance force on dynamic rotational error can be analyzed numerically. Assuming that the frequency of the disturbance force is an integer multiple of the fundamental frequency, the frequency order ranges from 2 to 10 UPR. Polar plots of the dynamic rotational error, which are caused by the disturbance force with the amplitude of 200 N, are shown in Figure 14. The frequency domain analysis shows that it mainly includes two peaks. The frequency of the first peak corresponds to the excitation frequency of the spindle dynamic unbalance force, that is, the fundamental frequency (spindle rotational frequency). The frequency of the second peak corresponds to the spindle resonance frequency (natural frequency of the SRB system). The asynchronous error motion presents large value due to the resonance frequency is non-integer multiples of the fundamental frequency.

Influence of Disturbance Force on Dynamic Rotational Error
The spindle is rotated by the drive system. The mechanical spindle often uses a wedge belt or a timing belt drive, while the electric spindle is driven by built-in motor. During the rotation of the spindle, the transmission system and the electromagnetic force will generate a disturbance force on the SRB system. The dynamic response of the SRB system to disturbance forces is also reflected in the dynamic rotational error. Based on the established model, the influence of disturbance force on dynamic rotational error can be analyzed numerically. Assuming that the frequency of the disturbance force is an integer multiple of the fundamental frequency, the frequency order ranges from 2 to 10 UPR. Polar plots of the dynamic rotational error, which are caused by the disturbance force with the amplitude of 200 N, are shown in Figure 14. As shown in the polar plots, the influence of external disturbance force, such as the drive system, on the dynamic rotational error is mainly forced vibration, and the frequency of the radial synchronous error is the same as the frequency of the disturbance force. By increasing the amplitude of the disturbance force, the comparison results of the dynamic rotational errors of different orders and amplitudes can be obtained, as shown in Figure 15. The simulation result shows that both synchronous and asynchronous error motion values increase slowly with the increasing of disturbance force order. When increasing the amplitude of the disturbance force from 200 to 400 N, the range of synchronous error motion values increase from 2.13-2.99 to 4.26-6.03 μm, and the range of asynchronous error motion values increase from 0.28-1.63 to 0.47-3.26 μm. Both the synchronous error and the asynchronous error increase as the amplitude of the disturbance force increases.

Experimental Validation
The experimental platform is built up to further verify the validity of the established model. Figure 16 shows the photo of the experimental setup. A mechanical machine tool spindle, which is driven by belt driving system, is applied to explore these issues. A high-precision dual ball artifact mounted in the rotor is used as the testing target, through which eccentricity can be adjusted manually. The NI (National Instruments) data acquisition (DAQ) card (USB 6356) and capacitive displacement sensor with nanometer-level resolution are used to acquire the displacement signals. Table 2 shows the technical parameters of the capacitive sensors used in the experiment. Data As shown in the polar plots, the influence of external disturbance force, such as the drive system, on the dynamic rotational error is mainly forced vibration, and the frequency of the radial synchronous error is the same as the frequency of the disturbance force. By increasing the amplitude of the disturbance force, the comparison results of the dynamic rotational errors of different orders and amplitudes can be obtained, as shown in Figure 15. As shown in the polar plots, the influence of external disturbance force, such as the drive system, on the dynamic rotational error is mainly forced vibration, and the frequency of the radial synchronous error is the same as the frequency of the disturbance force. By increasing the amplitude of the disturbance force, the comparison results of the dynamic rotational errors of different orders and amplitudes can be obtained, as shown in Figure 15. The simulation result shows that both synchronous and asynchronous error motion values increase slowly with the increasing of disturbance force order. When increasing the amplitude of the disturbance force from 200 to 400 N, the range of synchronous error motion values increase from 2.13-2.99 to 4.26-6.03 μm, and the range of asynchronous error motion values increase from 0.28-1.63 to 0.47-3.26 μm. Both the synchronous error and the asynchronous error increase as the amplitude of the disturbance force increases.

Experimental Validation
The experimental platform is built up to further verify the validity of the established model. Figure 16 shows the photo of the experimental setup. A mechanical machine tool spindle, which is driven by belt driving system, is applied to explore these issues. A high-precision dual ball artifact mounted in the rotor is used as the testing target, through which eccentricity can be adjusted manually. The NI (National Instruments) data acquisition (DAQ) card (USB 6356) and capacitive displacement sensor with nanometer-level resolution are used to acquire the displacement signals. Table 2 shows the technical parameters of the capacitive sensors used in the experiment. Data The simulation result shows that both synchronous and asynchronous error motion values increase slowly with the increasing of disturbance force order. When increasing the amplitude of the disturbance force from 200 to 400 N, the range of synchronous error motion values increase from 2.13-2.99 to 4.26-6.03 µm, and the range of asynchronous error motion values increase from 0.28-1.63 to 0.47-3.26 µm. Both the synchronous error and the asynchronous error increase as the amplitude of the disturbance force increases.

Experimental Validation
The experimental platform is built up to further verify the validity of the established model. Figure 16 shows the photo of the experimental setup. A mechanical machine tool spindle, which is driven by belt driving system, is applied to explore these issues. A high-precision dual ball artifact mounted in the rotor is used as the testing target, through which eccentricity can be adjusted manually. The NI (National Instruments) data acquisition (DAQ) card (USB 6356) and capacitive displacement sensor with nanometer-level resolution are used to acquire the displacement signals. Table 2 shows the technical parameters of the capacitive sensors used in the experiment. Data acquisition and processing software is developed by LabVIEW software and can be used for online analysis. The displacement data can be also saved for off-line processing.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 24 acquisition and processing software is developed by LabVIEW software and can be used for online analysis. The displacement data can be also saved for off-line processing.  In order to investigate the dynamic rotational error quantitatively, two experiments are designed. The schematic diagram of unbalanced mass and disturbance force verification experiment is illustrated in Figure 17. In the actual operation of the SRB system, there are many factors that affect its dynamic rotational error. It is difficult to eliminate the impact of these factors effectively. To study the influence of a certain factor, the proportion of this factor can be artificially enlarged. Compared with other factors, it will produce more obvious changes. When the SRB system is not subjected to any artificial excitation force, rotational error motion can also occur due to the influence of some unknown factors. The spindle rotational error measured in this state is used as the initial value for comparison. The polar plot is shown in Figure 18. The whirling speed is 300 rpm. The corresponding TIR, synchronous and asynchronous errors are 14.41, 0.55 and 2.86 μm, respectively.  In order to investigate the dynamic rotational error quantitatively, two experiments are designed. The schematic diagram of unbalanced mass and disturbance force verification experiment is illustrated in Figure 17. In the actual operation of the SRB system, there are many factors that affect its dynamic rotational error. It is difficult to eliminate the impact of these factors effectively. To study the influence of a certain factor, the proportion of this factor can be artificially enlarged. Compared with other factors, it will produce more obvious changes.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 24 acquisition and processing software is developed by LabVIEW software and can be used for online analysis. The displacement data can be also saved for off-line processing. Figure 16. Test setup of disturbance force verification experiment. In order to investigate the dynamic rotational error quantitatively, two experiments are designed. The schematic diagram of unbalanced mass and disturbance force verification experiment is illustrated in Figure 17. In the actual operation of the SRB system, there are many factors that affect its dynamic rotational error. It is difficult to eliminate the impact of these factors effectively. To study the influence of a certain factor, the proportion of this factor can be artificially enlarged. Compared with other factors, it will produce more obvious changes. When the SRB system is not subjected to any artificial excitation force, rotational error motion can also occur due to the influence of some unknown factors. The spindle rotational error measured in this state is used as the initial value for comparison. The polar plot is shown in Figure 18. The whirling speed is 300 rpm. The corresponding TIR, synchronous and asynchronous errors are 14.41, 0.55 and 2.86 μm, respectively. When the SRB system is not subjected to any artificial excitation force, rotational error motion can also occur due to the influence of some unknown factors. The spindle rotational error measured in this state is used as the initial value for comparison. The polar plot is shown in Figure 18. The whirling speed is 300 rpm. The corresponding TIR, synchronous and asynchronous errors are 14.41, 0.55 and 2.86 µm, respectively. Appl. Sci. 2020, 10, x FOR PEER REVIEW 19 of 24 Figure 18. Spindle rotational error of the initial state.

Validation of Dynamic Unbalance Force
The spindle rotor used in the experiment is handled by dynamic balancing and it produces very little centrifugal force at high speeds. An unbalanced mass is mounted at the location of pulley to generate dynamic unbalance force when it rotates, as shown in Figure 17. The dynamic unbalance force generated by unbalanced mass increases as the increasing of whirling speed. The dynamic rotational error can be measured by the DAQ system.
In the experiment, the whirling speed of SRB system ranges from 300 to 6000 rpm. The length of sampling is set at 60 revolutions, with the sampling frequency of 15 kHz. The dynamic unbalance force generated at 300 rpm is about 0.54 N. The TIR, synchronous and asynchronous errors are selected as the initial value due to some unknown factors except unbalance force. The rotational errors generated in the process of increasing the speed are all considered to be caused by the dynamic unbalance force. The measured TIR, synchronous and asynchronous errors under different speeds, as well as the comparisons with the simulated values, are shown in Figure 19. It can be seen from Figure 19a that experimental results agree well with the simulated values. The simulated value of TIR varies from 14.53 μm at 300 rpm to 16.78 μm at 6000 rpm, while the experimental value varies from 14.41 μm to 17.91 μm. The difference between the simulated and experimental TIR values become large from 3000 to 6000 rpm, due to the vibration of transmission system. Figure 19b indicates that dynamic unbalance force basically has no influence on synchronous error, which is consistent with the simulated results. However, dynamic unbalance force has a

Validation of Dynamic Unbalance Force
The spindle rotor used in the experiment is handled by dynamic balancing and it produces very little centrifugal force at high speeds. An unbalanced mass is mounted at the location of pulley to generate dynamic unbalance force when it rotates, as shown in Figure 17. The dynamic unbalance force generated by unbalanced mass increases as the increasing of whirling speed. The dynamic rotational error can be measured by the DAQ system.
In the experiment, the whirling speed of SRB system ranges from 300 to 6000 rpm. The length of sampling is set at 60 revolutions, with the sampling frequency of 15 kHz. The dynamic unbalance force generated at 300 rpm is about 0.54 N. The TIR, synchronous and asynchronous errors are selected as the initial value due to some unknown factors except unbalance force. The rotational errors generated in the process of increasing the speed are all considered to be caused by the dynamic unbalance force. The measured TIR, synchronous and asynchronous errors under different speeds, as well as the comparisons with the simulated values, are shown in Figure 19.

Validation of Dynamic Unbalance Force
The spindle rotor used in the experiment is handled by dynamic balancing and it produces very little centrifugal force at high speeds. An unbalanced mass is mounted at the location of pulley to generate dynamic unbalance force when it rotates, as shown in Figure 17. The dynamic unbalance force generated by unbalanced mass increases as the increasing of whirling speed. The dynamic rotational error can be measured by the DAQ system.
In the experiment, the whirling speed of SRB system ranges from 300 to 6000 rpm. The length of sampling is set at 60 revolutions, with the sampling frequency of 15 kHz. The dynamic unbalance force generated at 300 rpm is about 0.54 N. The TIR, synchronous and asynchronous errors are selected as the initial value due to some unknown factors except unbalance force. The rotational errors generated in the process of increasing the speed are all considered to be caused by the dynamic unbalance force. The measured TIR, synchronous and asynchronous errors under different speeds, as well as the comparisons with the simulated values, are shown in Figure 19. It can be seen from Figure 19a that experimental results agree well with the simulated values. The simulated value of TIR varies from 14.53 μm at 300 rpm to 16.78 μm at 6000 rpm, while the experimental value varies from 14.41 μm to 17.91 μm. The difference between the simulated and experimental TIR values become large from 3000 to 6000 rpm, due to the vibration of transmission system. Figure 19b indicates that dynamic unbalance force basically has no influence on synchronous error, which is consistent with the simulated results. However, dynamic unbalance force has a It can be seen from Figure 19a that experimental results agree well with the simulated values. The simulated value of TIR varies from 14.53 µm at 300 rpm to 16.78 µm at 6000 rpm, while the experimental value varies from 14.41 µm to 17.91 µm. The difference between the simulated and experimental TIR values become large from 3000 to 6000 rpm, due to the vibration of transmission system. Figure 19b indicates that dynamic unbalance force basically has no influence on synchronous error, which is consistent with the simulated results. However, dynamic unbalance force has a significant effect on asynchronous errors. Asynchronous errors due to unbalanced mass increase with the speed increasing. The experimental value of asynchronous error varies from 3.21 µm at 300 rpm to 4.69 µm at 6000 rpm, while the simulated value varies from 3.22 to 4.12 µm.

Validation of Disturbance Force
The same idea is adopted in this section that the disturbance force is enlarged artificially. External radial dynamic loading is applied by using a signal generator, power amplifier and electromagnetic vibrator [34]. A signal with known frequency is generated by the signal generator and amplified by a power amplifier to drive the electromagnetic vibrator to generate a dynamic disturbance force. A deep groove ball bearing, which is utilized to apply dynamic radial force to the SRB system, is mounted on the end of the rotor [34]. The error inside the load bearing can be omitted. A force sensor is used to measure the disturbance force in time domain. The schematic diagram and the experimental setup are shown in Figure 17.
In order to eliminate the influence of whirling speed change, the speed is set at 300 rpm. The amplitude and frequency of disturbance force are changed manually. The sinusoidal signal generated by the signal generator ranges from 10 to 40 Hz. The amplitude of the disturbance force can be adjusted by using the voltage gain of the power amplifier, approximately between 50 to 100 N. The theoretical values of the dynamic rotational error caused by the disturbance force are calculated by using the actual values measured by the force sensor. The spindle synchronous errors of different exciting orders are shown in Figure 20.

Validation of Disturbance Force
The same idea is adopted in this section that the disturbance force is enlarged artificially. External radial dynamic loading is applied by using a signal generator, power amplifier and electromagnetic vibrator [34]. A signal with known frequency is generated by the signal generator and amplified by a power amplifier to drive the electromagnetic vibrator to generate a dynamic disturbance force. A deep groove ball bearing, which is utilized to apply dynamic radial force to the SRB system, is mounted on the end of the rotor [34]. The error inside the load bearing can be omitted. A force sensor is used to measure the disturbance force in time domain. The schematic diagram and the experimental setup are shown in Figure 17.
In order to eliminate the influence of whirling speed change, the speed is set at 300 rpm. The amplitude and frequency of disturbance force are changed manually. The sinusoidal signal generated by the signal generator ranges from 10 to 40 Hz. The amplitude of the disturbance force can be adjusted by using the voltage gain of the power amplifier, approximately between 50 to 100 N. The theoretical values of the dynamic rotational error caused by the disturbance force are calculated by using the actual values measured by the force sensor. The spindle synchronous errors of different exciting orders are shown in Figure 20. By changing the amplitude of the disturbance force, the synchronous and asynchronous errors caused by the disturbance force with different amplitudes and frequencies can be obtained, as shown in Figure 21. It can be seen that the amplitude of the disturbance force decreases as the increasing of frequency. So, the corresponding synchronous error also shows an overall downward trend. However, the asynchronous error does not take obvious change. This can be attributed to the fact that the range of disturbance force is not large enough.  By changing the amplitude of the disturbance force, the synchronous and asynchronous errors caused by the disturbance force with different amplitudes and frequencies can be obtained, as shown in Figure 21. It can be seen that the amplitude of the disturbance force decreases as the increasing of frequency. So, the corresponding synchronous error also shows an overall downward trend. However, the asynchronous error does not take obvious change. This can be attributed to the fact that the range of disturbance force is not large enough.

Conclusions
An analytical modeling method, based on dynamic of the SRB system which considers the bearing nonlinearity, is presented in this paper, so as to study the formative mechanism of dynamic rotational error and predict the rotational performance for the SRB system. Core conclusions of the study as a whole are as follows.
The presented modeling method for dynamic rational error, based on the dynamic model of the SRB system, is feasible and valid to analyze the formative mechanism of the rotational error motion. Surface waviness, dynamic unbalance force and disturbance force are involved in the established model. The surface waviness hardly causes synchronous error motion. However, when the difference between the order and the number of rolling elements is one, the SRB system would produce a very obvious asynchronous error motion. The effects of dynamic unbalance and disturbance force are verified by the comparison experiments. Dynamic unbalance force has little effect on the synchronous error; however, the TIR and asynchronous error increases with the increasing of rotational speed and unbalance mass. Both synchronous and asynchronous error motion values increase when the amplitude of disturbance force and its frequency increases. Besides, the frequency of the radial synchronous error is the same as that of the disturbance force. The merit of this modeling method is that it can provide designers and/or field engineers with an informative guideline to improve the dynamic rotational accuracy.

Study prospects:
Based on the modeling method in this paper, investigations about the effect of temperature and lubrication on the dynamic rotational error during the spindle operation will be emphasized in future studies.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature δi
The displacement vector of the inner ring ϕbk The position angle of the kth ball Nb Number of the rolling element of the bearing ω0 The rotation angular speed of the inner ring (rad/s) ωbk The orbital angular speed of the rolling element (rad/s) ωsk The self-rotation angular speed of the rolling element (rad/s)

Conclusions
An analytical modeling method, based on dynamic of the SRB system which considers the bearing nonlinearity, is presented in this paper, so as to study the formative mechanism of dynamic rotational error and predict the rotational performance for the SRB system. Core conclusions of the study as a whole are as follows.
The presented modeling method for dynamic rational error, based on the dynamic model of the SRB system, is feasible and valid to analyze the formative mechanism of the rotational error motion. Surface waviness, dynamic unbalance force and disturbance force are involved in the established model. The surface waviness hardly causes synchronous error motion. However, when the difference between the order and the number of rolling elements is one, the SRB system would produce a very obvious asynchronous error motion. The effects of dynamic unbalance and disturbance force are verified by the comparison experiments. Dynamic unbalance force has little effect on the synchronous error; however, the TIR and asynchronous error increases with the increasing of rotational speed and unbalance mass. Both synchronous and asynchronous error motion values increase when the amplitude of disturbance force and its frequency increases. Besides, the frequency of the radial synchronous error is the same as that of the disturbance force. The merit of this modeling method is that it can provide designers and/or field engineers with an informative guideline to improve the dynamic rotational accuracy.

Study prospects:
Based on the modeling method in this paper, investigations about the effect of temperature and lubrication on the dynamic rotational error during the spindle operation will be emphasized in future studies.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature δ i
The displacement vector of the inner ring φ bk The position angle of the kth ball N b Number of the rolling element of the bearing ω 0 The rotation angular speed of the inner ring (rad/s) ω bk The orbital angular speed of the rolling element (rad/s) ω sk The self-rotation angular speed of the rolling element (rad/s) α o Initial contact angle of the bearing ( • ) α ik Inner ring contact angle of the bearing of the kth rolling element ( • ) α ik Outer ring contact angle of the bearing of the kth rolling element ( • ) β k The angle between the spinning axis of the rolling elements and Z axis ( • ) D The diameter of the rolling element (mm) d m Pitch Diameter (mm) δ ik /δ ok The contact deformation between the inner/outer bearing ring and kth rolling element f i /f o Inner/Outer groove curvature radius coefficient r i /r o Inner/Outer groove curvature radius (mm) η i The distance between the rotating axis Z and inner raceway groove curvature center (mm) L ik /L ok The distance between the inner/outer raceway groove curvature center and the rolling element center before deformation l ik /l ok The distance between the inner/outer raceway groove curvature center and the rolling element center after deformation Q ik /Q ok The contact force between the inner/outer ring and the rolling element (N) M gk Gyroscopic moment of the kth rolling element (N.mm) K ok Coefficient of the load-displacement between inner ring and kth rolling element (N/mm1.5) K ik Coefficient of the load-displacement between outer ring and kth rolling element (N/mm1.5) λ ik , λ ok , The friction coefficients between the inner/outer raceway and kth rolling element F ck The centrifugal force of kth rolling element (N) F i The resultant force vector of the inner ring due to all of the contacting rolling elements (N) K b The dynamic supporting stiffness matrix of the bearing (N/µm) δ ik The displacements for inner raceway groove curvature center corresponding to the kth rolling element (mm) T k The transformation matrix ∆w k The additional displacements caused by the waviness of the bearing rings and rolling elements (mm) ∆D k The additional displacements due to roundness error of the bearing rings and rolling elements (mm) p ik , p ok The radial waviness of inner and outer raceways p ik , p ok The axial waviness of inner and outer raceways w ik , w ok The waviness of kth rolling elements A in , A out , B in , B in , C ba , C ba The amplitudes of the surface waviness n in , n out , n ba , The harmonic orders of the surface waviness φ in , ψ in , φ out , ψ out , ϕ ba The initial phases angle of the surface waviness

F b
The sum of the time-varying forces from all of the contacting rolling elements (N) q The dynamic displacement vector of the rotor α i Outer ring contact angle of the bearing ( • ) C The structural damping matrix of the system G The gyroscopic matrix of the system J s The inertia moment of the rotor about the W s -axis R LSC The radius the least squares circle