Analysis on Dynamic Contact Characteristics and Dynamic Stiffness Estimating Method of Single Nut Ball Screw Pair Based on the Whole Rolling Elements Model

: The purpose of this paper is investigating the characteristics of dynamic contact and dynamic stiffness of the single nut ball screw pair (SNBSP). Then a new sensorless method is proposed to extract the SNBSP dynamic contact stiffness of a mass production CNC machine tool feed system. First of all, the transformation relationship between each coordinate system of SNBSP is established. Secondly, the dynamic model of all ball–raceway contact pairs is established. Based on this, a dynamic contact stiffness model is established. The dynamic contact parameters are obtained by the numerical method. It is found that the influence of screw speed on screw and nut raceway normal force distribution are opposite. This will affect the variations of dynamic contact stiffness. It is also clear that the effect of axial load on dynamic stiffness is significant. Then, an effective method is proposed to estimate the dynamic contact stiffness of a mass production CNC machine tool feed system without any external sensors. The axial force of feed system is estimated by using torque current of servo motor. Current signals can be obtained through FANUC Open CNC API Specifications (FOCAS) library functions, and then dynamic contact stiffness can be calculated through the stiffness model without external sensors. Finally, a feed system dynamic model is built, and the contact model and sensorless stiffness estimating method are verified by experiments in this dynamic system.


Introduction
Machining accuracy is a key feature to measure the overall performance of CNC machine tools. Ball screws are widely used in CNC machine tool feed systems because they have great stiffness, small axial clearance, high bearing capacity and low frictional force [1]. Although the ball screws have used in CNC machine tool systems for many years, the contact parameters of the screw nut pair, especially dynamic contact stiffness for keeping the feed system's positioning accuracy and precision, are still not fully explained.
The contact stiffness of ball screw pair depends on ball and raceway contact state. However, the ball screw pair has complex surfaces, thus internal motion of rolling elements is difficult to directly analyze. Therefore, a lot of numerical analysis methods have been studied to analyze it. A half lead model of double nuts ball screws is built by Hirokazu Shimoda et al. [2], and load distribution properties were analyzed. It was assumed that the normal force of the contact pair in every half lead was the same, and the phenomenon could be found from the computational results that load distribution was uneven with the unit of half a lead. Bertolaso et al. [3] verified the ball screw nut pair's uneven distribution, through photoelasticity, tag tracking and the FEM numerical method based on the half lead model. However, the models proposed by the above scholars cannot describe the contact state of the whole ball-raceway pairs, and the research on the screw nut is still in the stage of static analysis. Huang et al. [4] introduced the method of medial axis transform (MAT), and analyzed normal force, contact angle and contact stress of contact surface between the ball and the raceway. Lin et al. [5] researched ball screw pair kinematic characteristics by using the Frenet-Serret coordinate system, and proved that a sliding state existed in the screw and nut raceway. At the same time, Lin et al. [6] also proposed three methods to obtain ball screws' mechanical efficiency. Based on his theory, the optimization design of the ball screw pair mechanical structure was carried out. According to the literature [5], Wei et al. [7] carried out further research, and analyzed the dynamic variation of the contact angle, frictional force and slip angle of the contact surface between the ball and the raceway. However, the above scholars' ball screw pair dynamic analysis was in view of the assumption that the loading force of all the ball-raceway contact pairs is evenly distributed, which is contrary to the research results of Hirokazu and Bertolaso. Xu [8] investigated friction characteristics. The effect of creep parameters on the ball screw friction and friction distribution were analyzed. Chen [9] researched rolling element bearing, creep in one lap of ball screw nut pair, but he did not analyze all laps of the pairs. The axial, torsional and bending dynamic models of the ball screw were built, and the stiffness, mass and damping matrix of the system were derived [9,10]. Wu et al. [11] established the mathematical model of the axial and torsional stiffness of the ball screw based on the fixed constraint of the ball screw. By using Simulink, a simulation model of CNC machine tool feed system considering the stiffness was established. However, the contact stiffness model of ball screw pair is not accurate because all contact states of ball and raceway are not considered. Moreover, the extraction of dynamic contact stiffness of ball screw pair has not been applied to mass production CNC machine tools.
In view of the above problems, a dynamic contact model of a single nut ball screw pair (SNBSP) was proposed. Dynamic contact characteristics of an SNBSP were analyzed. Then, the stiffness sensorless estimating method of SNBSP was carried out. This paper consists of the following parts: In Section 2, dynamic contact mechanics model of an SNBSP was proposed. In Section 3, the influence of dynamic factors on contact normal force, contact angle and contact stiffness distribution characteristics of contact pairs are analyzed. In Section 4, a stiffness sensorless estimating method of SNBSP was proposed. An experiment is carried out by which the model is verified. Finally, Section 5 contains the conclusions.

Setting and Transformation of the Coordinate System
For researching the movement and contact characteristics of a screw and nut conveniently, three coordinate systems are introduced in this paper.
Rotating coordinate system O xyz : The direction of z axis is consistent with that of screw, and its rotating speed will be the same as the screw speed ( Figure 1). : In order to describe the contact pair during movement, the Frenet frame is selected and the Cartesian coordinate system is carried out to establish CS O nb τ ′ . The center of the coordinate is the ball center. The direction of the n axis is consistent with the screw radial direction at the O′ point. The tangent direction of the helix at point O′ is the τ axis. The direction of the b axis is determined by the τ axis and n axis through the right-hand rule ( Figure   1).
where λ denotes the lead angle of the ball screw, m r denotes the screw pitch radius, and θ denotes the position angle of the ball.
where, when Here, set the tiny position and the orientation of SNBSP in CS O xyz is can be express as follows [13], The following formula can be obtained through differential transformation theory Suppose the differential transform operators, respectively 0 0 Δ =Δ − P I (10) 2 2 ς ς Δ =Δ − P I (11) Put Equation (10) and Equation (11) into Equation (9), and transform it; 2ς Δ can be represented by The right side of the Equation (12) (14) where ( ) where, when

Dynamic Contact Mechanics Model of SNBSP
In the quasi-static model in which the screw speed is not considered, it is considered that the normal forces and contact angles of the ith ball and the screw and nut raceway are equal, that is, Then the relationship between the axial load and the normal force of contact pair of balls and raceway can be represented as [12], where Fa denotes the axial force; i Q ς is the normal contact force of the ith contact pair;  But the rotating screw will drive the balls to generate the angular velocity of the revolution around screw axis and the angular velocity of the ball rotation, respectively, which will make the balls produce outward centrifugal force C F , and produce the transverse frictional force denotes the transverse friction of the ith contact pair, and Ci F denotes the centrifugal force generated by the rotation of the ith ball. This can be represented as: where B m is the ball mass; m d is the screw helix diameter; and mi ω is the angular velocity of ball revolution along the screw helix of the ith ball. Then, the relationship between the axial force of the screw nut pair and the loads of the ball and raceway contact pairs can be represented as follows: 1 1 1 sin cos cos cos sin It is known from the deformation relationship between the ith and (i-1)th contact pairs that the deformation increment of the contact pairs of the adjacent balls can be expressed respectively as: where L Δ denotes the axial distance between adjacent balls, as shown in Figure 5. S E and N E denote the Young modulus of the screw and nut materials, respectively. S A and N A denote effective cross-sectional area of the screw and nut, respectively. Substituting Equation (23) into Equation (24a), it can be obtained as follows: Substituting Equation (23) into Equation (24b), it can be obtained as follows: According to the deformation coordination relationship, the deformation increment of the ith and the (i-1)th ball screw contact pair and ball nut contact pair can be represented as: Then, according to the Hertz contact theory, there is a relationship between the deformation and the contact force [14,15]: where B υ , S υ and N υ denote the Poisson ratio of the ball, screw and nut, respectively. , , can be obtained through the table in reference [16]. Simultaneous Equations (26a) and (27a), and substituting Equation (28) into the equation, it can be obtained as follows: In the same way, simultaneous Equations (26b) and (27b), and substituting Equation (28) into the equation, it can be obtained as follow: where: As a high-speed ball screw nut pair, on the one hand, it mainly bears axial force during operation. On the other hand, the influence of centrifugal force caused by the high-speed rotation of the ball should also be considered. The effect of contact angle, contact force and centrifugal force is considered comprehensively, and the deformation is shown in Figure 6a [17,18]. For analyzing the motion state and deformation of the ball screw nut pair in the motion process, the center of the nut raceway N O is fixed. As shown in Figure 6b, when the preload is not considered, there is no force between the ball and the screw and nut raceway before the SNBSP under load. So, there is no contact deformation, and the ball center is aligned with the screw and nut raceway curvature center. The internal and external contact angles on both sides of the ball and the raceway are the initial design angle 0 α . When the ball screw is bearing a load, due to internal (screw) and external (nut) contact force and the centrifugal force and the influence of internal and external contact deformation, the ball center moves from the point B O to B O′ , the screw raceway curvature center moves from the point S O to S O′ , and the contact angles of the ball and screw raceway and nut raceway change from 0 α to S α and N α , respectively. It can be seen from the geometric relation that, where S t and N t are the curvature radius coefficient of the inner and outer raceways, respectively, The geometric analysis of the ith ball after bearing load is as follows: Then, the dynamic contact mechanics model of SNBSP is established which contains Equations (20), (29), (31) and (34a)-(34f).

Calculation of the Frictional Force of SNBSP and Other Related Parameters
where a ς and b ς denote the semi-major and semi-minor axis of contact ellipse, respectively, the detailed calculation method can be detailed in the reference [16,19]. q′ ′ denotes the dimensionless coefficient in the direction of x ς ′′ axis, = / q x a ς ′′ ′′ . t ′′ denotes the dimensionless coefficient in the direction of y ς ′′ axis, = / t y b ς ′′ ′′ . X ς τ ′′ and Y ς τ ′′ denote the component of the shear stress of the ball and screw or nut in the direction of the x ς ′′ axis and y ς ′′ axis, respectively.
= cos( ) where ς τ denotes the shear stresses of the ball and screw or nut. ς ψ denotes the slip angle, which are described in Section 2.3.3.

Calculation of Shear Stress of the Contact Pairs
The shear stress ς τ generated in the contact area between the ball and the screw or ball and nut can be expressed as [4]: where f ς denotes the frictional coefficient of the contact pair. P ς denotes the Hertz contact pressure on the contact area. The relationship between P ς and Q ς can be represented as:

Calculation of the Contact Pairs Slip Angle
The slip angle of the contact pairs reflects the direction of the slip velocity at the contact point in the CS O x y z ς ς ς ς ′′ ′′ ′′ ′′ , as shown in Figure 7; it can be expressed as [4]: ′′ denote the component of the relative slip speed between the contact pair in the direction of the x ς ′′ axis and y ς ′′ axis, respectively. The relationship between the two parameters is as follows [5]: Relative slip velocity component of the ball screw contact pair: Relative slip velocity component of the ball nut contact pair: ω ω ω and m ω can be obtained in reference [7], so there is no repeat on it. The direction of each velocity is shown in Figure 2.

Dynamic Contact Stiffness Model of SNBSP
In the dynamic contact stiffness model, the rolling element is usually regarded as a spring without mass, which is connected along the contact point with the screw and nut raceway ( Figure 8). Stiffness refers to the ability of the material to resist the deformation of external forces. When the ball is in contact with the internal and external raceway, the screw raceway and the nut raceway all have contact deformation. For the ith ball, its normal contact stiffness with the raceway is represented as [20].
The contact stiffness matrix of the SNBSP can be derived though the transformation relations of , the stiffness of the ith ball with the raceway can be represented as: where i k ς can be obtained through Equation (41).
Then the contact stiffness between the ith ball and raceway in CS O xyz can be represented as The leading diagonal elements reflect the local variation of contact pairs. Therefore, the concept of local contact stiffness is proposed in this paper. The expression of local contact stiffness in all directions is shown in Equation (45).
In the ball screw pair, the relationship between each ball and the raceway contact pair is parallel. Therefore, the contact stiffness in all directions in CS O xyz is represented as:

Flow Chart of Numerical Calculation Method
The flow chart in Figure 9 describes the process of calculating the contact parameters of each contact pair by using the numerical solution based on the dynamic contact mechanics model of SNBSP in this paper.  Table 1 is the structure and size parameters of SNBSP. The parameters are taken into each equation, and the results of each dynamic contact parameter can be solved by iterative operation.

Variation of Dynamic Contact Normal Force
When the screw speed is 0, friction is not considered, Si In a previous study, it was shown that the normal forces decrease with the contact pair departing from the load end of the SNBSP. The normal forces of the ball increase simultaneously with axial force increase [12].
However, in the dynamic state, Si Ni Q Q ≠ . Figure 10 shows the variation of the dynamic contact force with the contact pair position and the screw speed when the axial force is 1000 N. Firstly, it is obvious from Figure 10 that the normal contact force of the nut is greater than that of the screw in the same position, and the higher the screw speed is, the greater the difference is. In different contact position, the contact force difference between the nut and screw is not the same; when the screw speed is below 1500 rpm, with the contact position moving backwards, the contact force difference between the nut and screw gradually decreases. When the screw speed is below 500 rpm, though still conforming to this rule, the relative deviation is less than 0.05%, therefore, it is considered that the screw and nut contact force are approximately the same at the speed below that speed. However, when the screw speed of the screw is higher than 1500 rpm, the contact force difference between the nut and the screw increases gradually as the contact position moves behind. The higher the screw speed, the more obvious this phenomenon. Secondly, it is consistent with the uneven distribution static contact force that the dynamic contact force decreases from the first contact pair to the last one. When screw speed increases, the level of uneven load distribution of screw increases slightly but that of the nut decreases. When the axial force is 1000 N, the screw speed is 500 rpm, the first and the last contact force ratio of the screw is 1.1422, and that of the nut is 1.1427. When the screw speed reaches 5000 rpm, the ratio of the screw is 1.1427, and that of the nut reduces to 1.1387. Figure 11a,b exhibit variation of contact normal force with the contact position and the axial force when the screw speed is 2000 rpm. It is clear that, when the axial force increases, the uneven load distribution of both screw and nut are gradually increased.  Figure 12 illustrates the variation of the dynamic contact angle of the first contact pair with axial force and screw speed. It is clear that under the same axial force, the screw dynamic contact angle increases with the screw speed increase. However, the variation of dynamic contact angle of the nut is just the opposite. The contact angle of the nut increases with the axial force increase. When the screw speed is lower than 2000 rpm, the contact angle of the screw increases with the increase in axial force. But when the speed of the screw is higher than 2500 rpm, it decreases first and then increases with the increase in axial force. This is because the factor of screw speed has a coupling effect on the dynamic contact angle. At the same time, it can be seen that the higher the screw speed is, and the smaller the axial force is, the greater the difference in the contact angle between the screw and the nut. However, the lower the screw speed is, and the smaller the axial force is, the closer the two dynamic contact angles. Figure 13 shows the variation of dynamic contact angle with the contact position and the screw speed when the axial force is 1000 N. At the same speed of screw, the dynamic contact angles of screw and nut both decrease from the first contact pair to the last one. However, the characteristics of the two are not the same. The uneven distribution of the dynamic contact angle of the screw is not obvious with the increase in the screw speed, while that of the nut is opposite. The reason is that the ball's spinning angle velocity ( R ω ) and revolutional angle velocity ( m ω ) are the two important factors which affect the dynamic contact angle. The bigger the two velocities are, the greater the screw contact angle is and the smaller the nut contact angle is. Figure 14a illustrates the variations of R ω and m ω with screw speed under different axial force. It is obvious that both two angle velocities are almost linear with screw speed. The axial force has little effect on the two angular velocities. Figure 14b shows variations of R ω with the contact position. It is clear that when the screw speed is higher than 2000 rpm, R ω increases from the first contact pair to the last one. When the screw speed is higher than 4000 rpm, m ω increases from the first contact pair to the last one. Normal contact force decreases with variation of the contact position, which is also the positive factor that makes the contact angle change with the screw speed. All of these factors make the distribution of the dynamic contact angle of the screw less obvious but that of the nut more obvious.   For the contact pairs in the same position, both the screw and nut dynamic contact angles increase with the axial force increase. For contact pairs in different position, dynamic contact angle decreases with the contact pair position changing, which presents the phenomenon of uneven distribution. Under the same axial force, the uneven distribution of the screw dynamic contact angle is more significant. However, as the axial force increases, the uneven distribution of the dynamic contact angles of the lead screw will be more and more significant, while that of the nut is contrary.

Variation of Dynamic Contact Stiffness
Substituting i Q ς into Equation (41), normal stiffness i k ς will be calculated. Figure 16 reveals the variation of the normal stiffness of the first contact pair with the axial force and the screw speed, respectively. It can be seen from the figure that the general variation rules of the screw and nut's normal contact stiffness are consistent. The normal contact stiffness increases as the axial force increases. The normal contact stiffness of the screw is less than that of the nut. With the screw speed increase, the normal contact stiffness decreases gradually, but the influence of the screw speed on the contact stiffness is less than that of the axial force.  Figure 17 exhibits variation of the normal contact stiffness with the screw speed and the contact position when the axial force is 1000 N. It is obvious that as the contact position is moving backward, the normal contact stiffness of both the screw and nut decreases gradually. In the same position, the normal contact stiffness is gradually reduced with the screw speed increase. The nut normal contact stiffness is always greater than screw. In the static state, the ball screw local contact stiffness in five directions can be calculated through Equations (42), (43) and (45). Figure 18 exhibits the variation of the local stiffness with contact position; the local lateral stiffness and oscillating stiffness fluctuate with the period of π and the phase difference of π/2. The amplitudes of lateral stiffness decrease slightly with the contact pair position far away from the SNBSP load side, but that of oscillating stiffness increases, as illustrated in Figure  18a,c. The load axial stiffness decreases with the contact pair position far away from the SNBSP load side, as shown in Figure 18b.     Figure 20a, it can be seen that the lateral stiffness increases gradually with the increase in the screw speed. The main reason for this phenomenon is that with the increase in screw rotational speed, the spinning and revolution angular velocities of the ball increase, and the ball moves to the outside gradually, which means that the nut contact angle decreases and the screw contact angle increases slightly. A closer contact state between the ball and the screw and nut raceway is generated, which leads to stronger support ability in the lateral direction. Therefore, lateral stiffness becomes bigger.  Figure 20b illustrates the variation of the local axial stiffness with the screw speed and the contact position. As the contact position moves backwards, the local axial stiffness decreases gradually. However, as the screw speed increases, the axial local stiffness gradually decreases. The reason for this is mainly that the screw and nut's normal contact stiffness decreases with the screw speed' increase. On the other hand, because the ball and raceway form a more close contact state which is caused by the changing of the contact angle, it leads to a stronger supporting ability along the lateral direction, and the axial support ability will decrease. Figure 20c exhibits the variation of the local oscillating stiffness with the screw speed and the position of the contact pair. It is different from lateral and axial stiffness that when the screw speed increases, the oscillating stiffness of the first three turns decreases, but the amplitude of the decrease gradually decreases with the phase fluctuation. The fourth half turn is a transition stage. From the fifth half turn, the contact pair oscillating stiffness increases gradually, and the amplitude of the increase gradually increases with the phase fluctuation. Figure 21 shows the difference value between the local oscillating stiffness of each contact pair at screw speeds of 5000 and 500 rpm when the axial force is 1000 N. Substituting local stiffness into Equation (46), and the screw-nut pair global stiffness will be calculated. As shown in Figure 22, the global contact stiffness in all directions increased with the increase in the axial force. Compared with the lateral stiffness x K and y K , the axial stiffness z K and the oscillating stiffness x K φ , y K φ increase more significantly, which indicates that the effect of the axial force on the latter two is more significant.  It is clear that with the rotating speed of the screw increase, the lateral stiffness x K and y K increase monotonically, while the axial stiffness z K decreases monotonically, which is in accordance with variation rule of the local contact stiffness. The oscillating stiffness x K φ and y K φ increase with screw speed increasing. Although the local oscillating stiffness of the first 1.5 turns decreases with screw speed increasing, it increases with the screw speed, increasing two cycles, and the increment is significant. So the total value of the oscillating stiffness increases monotonically, but it is clear that the amplitude of its increase is not very obvious.

Sensorless Stiffness Estimating Method Based on Contact Stiffness Model
Through the numerical analysis of the dynamic contact model, it is clear that the axial force and the speed of the screw have major influence on the dynamic contact stiffness of the SNBSP. At the same time, they are also important parameters to solve the model. If the values of them can be obtained when the machine tool is running, we can solve the dynamic contact stiffness of SNBSP according to the dynamic model proposed in this paper. It is significant to evaluate the dynamic characteristics of mass production CNC machine tools. The value of screw speed can be obtained through the CNC program. But the capture of axial force needs the help of an external sensor generally. That is why the stiffness analysis can only be carried out on components or a special test bench in references.
The relationship between the axial force of the SNBSP and the screw torque is as follows [21,22] where e k is the servo motor torque constant; i is the TCSM.
In this way, the problem of capturing the axial force is transformed into the problem of capturing the TCSM.
In a previous study, we developed a collection program by using VC + +. The working data of CNC system can be collected in real time through FOCAS1 database [23]. If the same type of data is collected, the sampling frequency can reach 28 Hz. For a more detailed processing, refer to reference [23]. In this way, an Ethernet cable is used to connect the CNC system with the PC, and then the TCSM signals at any time can be captured without the help of any external sensor. (Notes: A notebook PC was used in the earlier stage of study. However, it was replaced by a PC104 BUS computer for other functions.).
The process of estimating the dynamic contact stiffness of ball screw without a sensor can be seen in Figure 24.

Experiment
In order to verify the numeric results of the dynamic contact model and sensorless stiffness estimating method, an experiment is conducted to measure the axial stiffness of a SNBSP in the Zaxis feed system of a CNC lathe (Shenyang Machine Tool Group, Shenyang, China) with Fanuc Series (Figure 25). The experimental lathe parameters are exhibited in Table 2. The parameters of SNBSP(Best Precision, Jining, China) are exhibited in Table 1. In order to verify this, an external sensor is needed to measure the response signals during the operation of the machine tool. Fast Fourier transform (FFT) is used to process the signals, and the frequency components obtained are compared with the theoretical values to verification.
Hence an acceleration sensor, B&K4520-004 (Brüel & Kjaer, Naerum, Denmark), of which the maximum of sampling frequency is 25.6 kHz, was mounted on the worktable and moves together with it. The acceleration sensor was connected to NI-DAQ (National Instruments, Austin, State of Texas, US) which was connected to Notebook PC. The max travel of the worktable is 300 mm. However, considering the safety of high-speed operation, the travel was set to 270 mm with two end point coordinates 15 mm (denoted by Z0) and 285 mm (denoted by Z1), respectively. 285 mm The experimental process is as follows.
Step 1: Locate the worktable at the given position Z1 by using a NC program.
Step 2: Run acceleration sensor data acquisition program. Run servo axis data acquisition program to capture the TCSM.
Step 3: Run the NC program. The worktable moves from Z1 to Z0 at the specified speed.
Step 4: All the application programs are stopped and the data are saved.
Because the maximum feed speed of the Z feed shaft of the NC lathe is locked to 10,000 mm/min, and the lead of the Z axis screw shaft is 10 mm, that means that the screw speed is 1000 rpm. Hence, the high speed of the screw is selected as 1000 rpm, and 50 rpm is selected as the low speed. Under each rotating speed, the worktable is run three times according to the above test steps, in the forward (Z1-Z0) and reverse (Z0-Z1) direction. Then, the acceleration and TCSM signals are obtained. Figure  26 shows the acceleration and TCSM signals of the first positive feed at the screw speed 1000 rpm. It can be seen from Figure 26a that there are three states in the process: a state of acceleration, stable running and deceleration. The acceleration signals in the stability state can be used for FFT analysis.
It can be seen from Figure 26b that there are also three operation states in TCSM signals. Calculate the average value of the TCSM in the stability state. By substituting the average value of the current into Equations (47) and (48), the axial force of the Z axis screw of the feed system can be obtained as 637 N under the screw speed of 1000 rpm. According to the model presented in this paper, the dynamic axial contact stiffness of the SNBSP can be calculated. When the screw speed is 50 rpm, the axial force is about 684.8 N. Because the speed is very low, the axial stiffness of the SNBSP can be solved by static model, and the result is 3.89 × 10 5 N/mm.  Figure 27 shows the configuration of the z-axis feed system and the dynamic model. The mounting method of ball screw is fixed-fixed, as shown in Figure 27a. The stiffness values of the frame and worktable of the experimental machine tool are all large. Compared with the worktable mass, the mass of the screw can be neglected. Especially, low screw lead reduces the influence of the torsional stiffness and moments of inertia of the screw and other rotating parts. Therefore, the dynamic model of the worktable is simplified as in Figure 27b. The stiffness of the system is expressed as follows: where L χ denotes the distance between the nut and the two fixed ends, χ = 0 or 1. Then the feed system natural frequency is expressed as [24]: The axial stiffness of the feed system can be determined by numeric iteration of Equations (49) and (50).
It is clear from Equation (50) that the screw shaft axial stiffness varies with the position of the worktable, so the axial stiffness and natural frequency of the feed system are also time-varying. The theoretical value of the axial stiffness and natural frequency of the feed system variation with the position of the worktable under the experimental conditions are shown in Figure 28. It can be seen that although the axial stiffness and natural frequency vary with the position of the worktable, the variation is not obvious. For example, when the screw speed is 1000 rpm, the axial stiffness variation ratio is 0.087%, while the natural frequency variation ratio is only 0.046%. Therefore, the axial stiffness and natural frequency of the Z axis of the feed system can be considered to be approximately the same in different positions. The feed system natural frequency will be calculated by Equation (51). Herein, the contact stiffness and the natural frequency determined above are considered to be theoretical values of them, respectively.
On the other hand, the FFT algorithm is carried out to investigate the spectrum characteristics of the worktable acceleration. In Figure 29, the FFT results of the worktable acceleration are represented in the stable running state in Figure 26a. It is obvious that the 16.66 Hz and its two and three times harmonic frequency components are the frequency components corresponding to the screw speed of 1000 rpm, and there are several frequency components at the low frequency section. In fact, the feed system has shapes of many modes [10], and so the FFT results of its vibration contain multiple frequencies. However, the first mode of the feed system is generally axial [25]. Herein, the frequency in Figure 29, 108.6 Hz, should be the natural frequency, whose corresponding mode should be the axial vibration of the worktable. Therefore, 108.6 Hz corresponds to the natural frequency n f in Equation (51). The theoretical and experiment values of the feed system dynamic parameters under the two screw speeds in Section 4.1 are exhibited in Table 3. The results show that the theoretical values of natural frequency are consistent with the measured values. This fact verifies the effectiveness of the whole rolling elements dynamic contact model and sensorless stiffness estimating method proposed in this paper. This estimating method can be widely used in the mass production CNC machine tools.

Conclusions
In this paper, the dynamic contact mechanics model and dynamic contact stiffness model of all balls of screw-nut pairs are established. Based on this, a sensorless stiffness estimating method is proposed. The following conclusions are obtained: (1) In the dynamic state, Si Ni Q Q ≠ . In the same position, the normal contact force of the nut is greater than that of the screw. The higher the speed of the screw, and the larger the axial force, the greater the difference between the two normal forces. When the screw speed is below 500 rpm, it is considered that the contact normal force of the screw and nut are approximately the same. Screw speed has different effects on the uneven distribution of the load of the screw and nut.
(2) In the dynamic state, s n i i α α ≠ . In the same position, the contact angle of the screw is greater than that of the nut. The higher the speed of the screw, and the smaller the axial force, the greater the difference between the two contact angles. Because the spinning and revolutional angle velocities of the ball are different at different contact positions, the coupling influence of the axial force and screw speed have different effects on the uneven distribution of the dynamic contact angle of the screw and nut.
(3) The effect of the screw speed on the lateral stiffness and axial stiffness is more significant, but the influence on the oscillating stiffness is little. However, the influence of the axial force on all the stiffnesses is significant.
(4) A stiffness sensorless estimating method of SNBSP was proposed. The axial force of SNBSP in feed system can be captured without a sensor when the machine tool is running. Then, the stiffness of SNBSP can be obtained. This method is verified by an experiment. Hence it is an effective method which can be used in mass production CNC machine tools.

Calculation Method of Transmission Efficiency
Transmission efficiency η is affected by axial force. The calculation is based on the theory of references [6,21]. In Ref. [21], Wei proposed the theoretical formula of transmission efficiency η of ball screw pair which was verified to predict ball screw performance. He also used Equation (A4) to calculate the torque Tca, which was proposed in Ref. [6]. Here, these core equations are listed. (Here, we unify the symbols in the equations in Refs. [6,21] with those in this paper).
Therefore, combined with the model and method proposed in this paper, the above equations are used to calculate the transmission efficiency η of ball screw pair. The algorithm flow chart is shown in the Figure A1.
Firstly, the initial value η 0 is given according to the classical equation of transmission efficiency. According to Equation (47) in this paper, Fa 0 is calculated. The parameters such as Fa 0 and screw speed are substituted into the dynamic contact model, and results such as