Multi-Position Identification of Joint Parameters in Ball Screw Feed System Based on Response Coupling

Existing methods of parameters identification do not consider the torsion characteristics of a ball screw and the worktable position simultaneously. Therefore, this paper proposes a multi-position identification method based on receptance coupling. Firstly, the mathematical model of the feed drive system is established by the improved receptance coupling, and this model considers both axial and torsional vibration of the ball screw. Secondly, the identification equation is established by minimum error of the modal parameters of multiple worktable position, and differential evolution algorithm is used to calculate the stiffness and damping of the joint. Finally, the self-developed ball screw feed drive system is used for experimental study. The maximum error of the first four orders of natural frequencies predicted through multi-position identification results is 2.95%, and the multi-position method is more robust than the common method identification at one position. The experiment study showed that the proposed method is accuracy and necessity.


INTRODUCTION
Ball screw feed drive system is extensively used in CNC machine tools [1,2]. Reliable dynamic models, dynamic characteristics analysis, and servo controller design for ball screw feed drive system are based on accurate joint parameters [3,4]. The accuracy of joint parameter identifications becomes increasingly important with the development of machine tools to high speed and high precision.
The joints of the ball screw feed drive system have two types, namely, fixed and rolling. The stiffness of the fixed joint is sufficiently large that the fixed joint is typically regarded as rigid coupling [5]. Moreover, parameter identification of the rolling joints is complex and challenging. The parameter identification of the rolling joint of the ball screw feed drive system can be classified as separate and simultaneous. In separate identification, the modeling and identification of rolling joints are considered separately [6][7][8][9]. Results inevitably contain errors to a certain degree given the coupling of components. Thus, simultaneous identification methods have been proposed in the literature to overcome this limitation. At present, there are few researches on this aspect. Deng et al. [10] developed a finite element (FE) model for ball screw feed drive system, where an objective function is established with the minimum error between the measured and simulated values of natural frequency and damping ratio, and the joint parameters of the feed drive system have been identified using an optimization algorithm. Due to the complex structure of the ball screw pair, whose FE model usually only considers the axial and radial stiffness of the joint, and this induces that the torsional characteristics of the ball screw are not considered in the model. Therefore, the FE model is often used in the whole machine model (the influence of ball screw torsional characteristics is small and can be ignored) while it is seldom used for the modeling of the feed drive system because the influence of ball screw torsional characteristics is large and cannot be ignored. Wang et al. [11] proposed a stiffness calculation method that uses Hertz contact theory to calculate the stiffness of the joints of the ball screw feed drive system. Some parameters in the method of the Hertz contact theory are design parameters of enterprises, which are not obtained easily, and the method cannot afford continuous analysis according to the position change of the worktable. Zhu et al. [12] established a mechanical model for ball screw feed drive system and identified the joint parameters through a genetic algorithm. In this model, the position of the worktable, frequency of the excitation force, and vibration amplitude of the ball screw are considered. However, the model does not consider the torsional characteristics of ball screw. In addition to the axial characteristic of the ball screw, its torsional characteristics and the position of the worktable affect the dynamic characteristics of the feed drive system [13,14] and must be considered during the identification of the joint parameters. However, existing research does not consider the torsional characteristics and worktable position simultaneously. The present study aims to achieve an accurate identification method for joint parameters that consider the position change of a worktable, the torsion characteristics and axial characteristics of a ball screw.
Receptance coupling method is a common method for parameter identification of joint. This method is rarely used in ball screw feed drive system and is mainly used to identify the joint parameters of a bolting system [15,16] and the milling machine spindle system [17,18]. Liu et al. [19] established the mathematical model of the two-axis feed drive system by using the receptance coupling, and analyzed the influence of worktable positions on the dynamic characteristics of the ball screw feed drive system. However, this study did not introduce how to obtain the joint parameters. Thus, the work of Liu et al. inspired us to utilize the receptance coupling method to identify the joint parameters of the ball screw feed drive system, while considering the position of a worktable and the torsion characteristics of a ball screw. Although Liu established the mathematical model of the two-axis feed drive system by the receptance coupling, there are still the following problems. At present, the commonly used model of ball screw feed drive system is a rigid-flexible model. The worktable is regarded as rigid body, without FRF. While modeling by the receptance coupling is coupling the FRF of substructure. Liu regards the carriage, the worktable and the joint as one substructure, and establishes the model of feed drive system by coupling the FRF of this substructure. However, the single-axis feed drive system has a worktable, which cannot form a substructure and cannot be modeled by this method. However, the modeling and the dynamic characteristic analysis of single-axis feed drive system are equally important [20]. Therefore, this paper proposed a new modeling method according to the idea of lumped parameter model of feed drive system. The ball screw is divided into several segments, which are regarded as flexible body and rigid body respectively. The flexible body is simulated by a shaft with axial and torsional vibration. The rigid body and the worktable are regarded asa substructure. This method can be applied not only to the modeling of single-axis feed drive system, but also to the modeling of two-axis feed drive system.
The objective of this paper is to consider the influence of the axial and torsional vibration of the ball screw and the worktable position on the feed drive system comprehensively. Firstly, this paper improves the modeling method of receptance coupling, and solves the problem that the existing method cannot be applied to the modeling of single-axis feed drive system. Secondly, the multi-position identification is proposed for the joints of feed drive system, and the joint parameters are identified. At last, a self-developed ball screw system is used for experimental study. The results showed that the multiposition identification method has higher precision and better robustness.

MODELING OF A BALL SCREW FEED DRIVE SYSTEM
The typical structure of a ball screw feed drive system is illustrated in Fig. 1. The servomotor drives the ball screw rotation through the coupling. The ball screw is supported on the bed through the bearing group at both ends. The screw nut is rigidly coupled with the worktable. The worktable is supported on the bed through guides. The screw nut turns the rotation movement of the ball screw to the axial movement of the worktable. Its basic principle is to transform torsion force into axial force by the ball screw pair. Given that the ball screw is a slender shaft, the flexible performance in axial and torsional directions under the action of torsional and axial forces cannot be ignored [13]. In the present study, the ball screw feed drive system is simplified into an axis-torsion model (Fig. 2) while considering the axial and torsional deformations of the ball screw. In Fig. 2, x s is the axial displacement of a ball screw. θ is the torsional displacement of the ball screw. k al and c al are the axial stiffness and damping of the left bearing group, respectively. k ar and c ar are the axial stiffness and damping of the right bearing group, correspondingly. k m and c m are the torsional stiffness and damping of coupling, respectively. k n and c n are the axial stiffness and damping of ball screw pairs, correspondingly.

Substructure Division of the Ball Screw Feed Drive System
The ball screw feed drive system must be divided into several substructures to apply the receptance coupling method. Its principle and aim are the conveniences for solving the FRF of each substructure. In this study, the ball screw feed drive system is divided into three substructures, namely, axial, torsional, and worktable substructures, as depicted in Fig. 3.
The ball screw is divided into four segments. Segment F is the part contacted with the screw nut. Segment H is the part contacted with coupling. The segment of the ball screw between Segments F and H is Segment A. The rest of the ball screw is Segment B. Segments A and B are regarded as flexible bodies with axial and torsional vibration. Segments F and H are regarded as rigid bodies. C and D are the beds (Two joints exist between the bed and ball screw). Meanwhile, E is the worktable. The axial substructure includes the axial vibration model of Segments A and B and related joints. Segments A and B of the ball screw are simulated as the elastic shaft with axial deformation. The connection between the axial vibration model of Segments A and B is the rigid coupling.
The torsional substructure includes the torsional vibration model of Segments A, B, and H, motor, and coupling. Segments A and B are simulated as the elastic shaft with torsional deformation. Segment H and the motor are simulated as the lumped parameter. Moreover, the connection between the torsional vibration model of Segments A and B and that between the two torsional vibration models of Segments A and H is the rigid coupling.
The worktable substructure consists of Segment F and worktable. The modeling for Segment F contains the axial and torsional vibrations. The axial vibration of Segment F, torsional vibration of Segment F, and the worktable are modeled by lumped parameter. The connection of the axial vibration model between Segment F and axial substructure is the rigid coupling. The connection of the torsional vibration model between Segment F and torsional substructure is the rigid coupling.

Axial Substructure Modeling
According to the analysis in Subsection 2.1, the axial substructure model of the feed drive system can be obtained, as demonstrated in Fig. 4. The segments of the ball screw are marked with 1 -4, which are the excitation and response points of the FRF, to represent the FRF of the substructure.

Coupling the Axial Vibration Model of Segments A and B
By using Segments A and B as circular shafts with free ends, the axial FRF can be calculated directly by Eq. (1) to Eq. (2) [21].
where A and B indicate the substructures, and a is the axial direction of excitation and response. Moreover, H A−a is the axial FRF matrix of Segment A; h A−a12 is the axial FRF with response at Point 1 and excitation at Point 2; E is the elasticity modulus of the ball screw; S is the cross-section area of the ball screw; ω is the angular frequency; and ρ is the density of the ball screw. Then, λ can be expressed as Eq. (3).
where H AB−a is the axial FRF matrix of the entire structure composed of Segments A and B by rigid coupling. The connection between Segments A and B is the rigid coupling. The FRF coupling of Segments A and B can be obtained through the principle of receptance coupling [15][16].

Coupling the Axial Vibration Model of Segment ABand
Structure C Bed C is regarded as the rigid body. Thus, the axial FRF of Structure C can be expressed as The joint between Segments AB and Structure C is simulated as a spring damper. k al and c al are the stiffness and damping, respectively.
H AC−a can be obtained by coupling the FRF of Segments AB and Structure C.

Coupling the Axial Vibration Model of Structures ABC and D
The coupling method in this subsection is similar to that in Subsection 2.2.2. The only difference is as follows. In Subsection 2.2.2, Structures C and AB are coupled at Point 1, whereas Structures D and ABC are coupled at Point 4 in this section. The stiffness and damping of the joint between Structure ABC and D, correspondingly, are represented by k ar and c ar . Then, the FRF of its joint can be expressed as H AD−a is the FRF matrix of the axial substructure. It can be obtained by coupling the FRF of all four sections, as expressed in Eq. (9).

Torsional Substructure Modeling
According to the analysis in Subsection 2.1, the torsional substructure model of the feed drive system can be obtained, as exhibited in Fig. 5. Numbers 1 -4 represent the excitation and response points of the FRF.

Coupling Torsional Vibration Model of Segments A and B
Segments A and B of the ball screw are circular shafts with free ends. The torsional FRF can be calculated using Eq. (10) and Eq. (11) [21].
In Eq. (10) and Eq. (11), t indicates the torsion direction of excitation and response; H A−t and H B−t are the torsion FRF matrices of Segments A and B, respectively; G is the shear elasticity; and d is the diameter of the ball screw. J can be defined as: The connection between Segments A and B is rigid coupling. The FRF coupling of Segments A and B can be obtained as:

Coupling Rotor
The rotor of the motor and Segment H of the ball screw are connected by coupling, as displayed in Fig. 5. This part is modeled by lumped parameter. According to finite element theory, the mass matrix M H , stiffness matrix K H , and damping matrix C H can be calculated as follows: where J m is the inertia of the rotor; J sH is the inertia of Segment H. The FRF in this part is calculated and represented as h H−t00 , which means the torsional response at Point 0 of the torsional excite at Point 0. The connection between Segments A and H is the rigid coupling. The FRF of the torsional substructure can be obtained by coupling the torsional FRF of the motor with the torsional FRF of Structure AB, denoted by H AH-t .

Worktable Substructure Modeling
The model of the worktable substructure is presented in Fig. 6, which consists of the axial model of Ball screw segment F, the torsional model of Ball screw segment F, and the axial model of the worktable. This substructure is modeled by lumped parameter.  Fig. 6, 5a indicates the axial model of Segment F in the ball screw, 5t represents the torsional model of Segment F in the ball screw, 6a means the model of the worktable. Moreover, m s is the mass of segment F, J s is the inertia of segment F, m b is the mass of the worktable, x s is the axial displacement of segment F in the ball screw, θ s is the torsion displacement of Segment F, and x b is the axial displacement of the worktable.
By applying the Lagrange energy method to the dynamic equation of the worktable substructure [22][23], the mass, stiffness, and damping matrix of the substructure are expressed as: The coefficient i b can be expressed as Eq. (21).
where h is the lead of the ball screw.
where W indicates the worktable substructure.

Modeling of the Ball Screw Feed Drive System
The mathematical model of the assembly system can be obtained by coupling the FRFs of the three substructures in Subsections 2.2 -2.4. Point 2 of the axial substructure is rigidly coupling with Point 5a of the worktable substructure. We obtain: Point 2 of the torsional substructure is rigidly coupling with Point 5t of the worktable substructure. The connection of the torsional substructure through rigid coupling yields: Eq. (24) is the mathematical model of the ball screw feed drive system illustrated in Fig. 2.

PROCESS OF JOINT PARAMETER IDENTIFICATION
According to the simplified model of ball screw feed drive system, the parameters to be identified are: k al , c al , k ar , c ar , k n , c n , k m , c m . In this paper, the differential evolution algorithm is used for identification [24], and the optimization variable is expressed as Eq. (25).
The optimization variable is substituted into the mathematical model of the feed drive system in Section 2, and the theoretical natural frequencies and corresponding amplitude at different positions are calculated and expressed as nf_c(p)(q) and amp_c(p)(q), respectively. Where p is the position of the worktable, and q is the order of the natural frequency. In this paper, the multi-position identification is used to improve the accuracy of the joint parameters of the feed drive system. Therefore, the experimental FRFs of different worktable positions are tested by modal experiment. Test the FRF (the axial response on the worktable, the axial excitation worktable) of system, which is hall−6a6a in Eq. (24). The natural frequency and the corresponding FRF value are extracted respectively from the tested data. The objective function was established as Eq. (26). It is important to note that the used FRF in the test and modeling is axial FRF of the feed drive system. Therefore, the natural frequency is also the axial natural frequency.
Where, m is total number of worktable positions, n is the total number of orders. The specific identification process as shown in Fig. 7.

EXPERIMENTAL VERIFICATION
In this section, the FRFs of the feed drive system are tested by modal experiment. Firstly, the natural frequencies of the feed drive system are determined by analyzing the experimental results, and the joint parameters are identified by multi-position identification. Secondly, the accuracy of the multi-position identification is verified by the consistency between the predicted FRF and the tested FRF. Finally, the robustness and necessity of the multi-position identification are verified by comparing the overall error between the single-and multi-position identification results.

Experimental and Result Analyses
The test setup is demonstrated in Fig. 8. Modal test system (LMS SADAS and LMS Test Lab), impact hammer (Kistler 9724A2000), and miniature sensor (B&K 4525b) are used in the experiment. The experiment is carried out on the test bench, as shown in Fig. 9.  Due to the influence of many factors, such as the joint or the nonlinearity of the system components, it is difficult to determine the axial natural frequency of the feed drive system by an experimental FRF. The typical method to determine the axial natural frequency is by comparing the FRFs with different excitation and response, and analyzing the modal shape of the system. However, the feed drive system has its own characteristics. Firstly, the modal shapes of the feed drive system are expressed as the axial and torsional modal shapes of ball screw, and it is difficult to determine the vibration mode of the system. Secondly, with the worktable position change, the natural frequencies of the system is roughly the same, or changes in a small range in a certain rule. Therefore, this paper proposes a new method that the axial natural frequencies are determined by comparing the same FRFs of different worktable position and according to the change rule of the FRF peaks. The measured FRFs in six positions are exhibited in Fig. 10. Fig. 10a to Fig. 10b display the details of the dashed box. In Fig. 10a, the peaks of all the FRFs are approximately 24 Hz and that of amplitude is nearly the same. Fig. 10b, the peaks of all the FRFs are approximately 115 Hz, and the amplitude and frequency of the peak increases and decreases, respectively, with the increase in the distance l. In Fig. 10c, the peaks of all the FRFs are approximately 540 Hz and the amplitude of the peak is nearly the same, whereas the frequency initially increases and then decreases with the increase in l. That is, the frequency of the peak is high when the worktable is close to the middle of the ball screw. In Fig. 10d, the peaks of all the FRFs are approximately 720 Hz. The amplitude and frequency of the peak are nearly the same. The peaks beyond the four frequencies of the FRF has no regularity; thus, the ball screw feed drive system has four natural frequencies within 0 -1000 Hz, which are near the four frequencies. The experimental data of the three positions (i.e., l 1 , l 3 , and l 5 ) are selected to identify the joint parameters. The natural frequencies and corresponding amplitudes of the FRF for the first four orders at these positions are summarized in Tab. 1 and Tab. 2.

Multi-Position Identification of Joint Parameters and Verification
According to the method discussed in Section 3, the theoretical and measured data of Positions l 1 , l 3 , and l 5 are substituted into Eq. (26). The model parameters used to calculate the theoretical data are listed in the Tab. A of the Appendix. The joint parameters of the ball screw feed drive system are identified using a differential evolution algorithm, as listed in Tab. 3.  The measured natural frequencies of the three other positions (i.e., l 2 , l 4 , and l 6 ) are used to verify the accuracy of the identification results. The comparison of the predicted and measured results are presented in Tab. 4. The overall error of a position is calculated using Eq. (27).
where q is the order of the natural frequencies. In Tab. 4, the maximum error of the natural frequency for the first four orders of the ball screw feed drive system is 2.95%. Due to l 4 being in the middle of the ball screw and being the vibration node of the third order (Fig. 12), the natural frequency of the third order cannot be determined from the FRF by using the measured and predicted FRFs. Thus, the result is excluded from Tab. 4. The accuracy of multi-position is confirmed by comparing the predicted and measured natural frequencies. Then, the comparison between the predicted FRFs at Positions l 2 , l 4 , and l 6 and the measured FRFs is exhibited in Fig. 11 to Fig.  13. The predicted FRFs of the three positions match well with the measured FRFs.

Results Analysis of Single-and Multi-Position Identification Methods
In this subsection, the advantages and disadvantages of the two methods are judged by comparing the prediction data of single-and multi-position identification results. Since the middle position is a node usually, the test position of modal experiment is not here usually, and the usual test position is set near the middle, such as position l 3 and l 5 . The experimental data of these two positions were used to identify the parameters by single-position identification, and make a comparison with multi-position identification.
There is a big difference of the results of one singleposition identification with others.Therefore, firstly, it is determined that the result of the single-position identification is the optimal solution of the identification equation by comparing the data of position l3 (l 5 ). Secondly, the accuracy and robustness of multi-position identification are determined to be higher by comparing the data of the other three positions (l 2 , l 4 , l 6 ).

Single-Position Identification and Result Analysis
The joint parameters are identified by the l 3 and l 5 position identifications. The results of the single-and multiposition identification are summarized in Tab. 5. In Tab. 5, the deviation between the l 3 -and multiposition identifications is small, whereas that of the l 5 -and multi-position identifications is large. Therefore, firstly, through comparing the theoretical and experimental values of positon l 3 (l 5 ), the result of single-position identification is determined to be the optimal solution of the identification equation. All the data are listed in Tab. 6 and Tab. 7, correspondingly. The FRFs are displayed in Fig.14 and Fig. 15, respectively. As the single-position identification considers the minimum error of one position, and the multi-position identification considers the minimum error of three positions, for the natural frequencies of position l 3 (l 5 ), the accuracy of single-position identification should be higher, and for the other three positions (l 2 , l 4 , l 6 ), the accuracy of single-position identification should be lower. As Tab. 6 and Tab. 7. In Fig. 14 and Fig. 15, the predicted FRF matches the measured FRF well, and the largest error between the identification result and measured data is only 2.12% and 2.14%, correspondingly. Thus, the results summarized in Tab. 5 are the optimal solutions of the identification equation.

Necessity of the Multi-Position Identification Method
The dynamic characteristics at Positions l 2 , l 4 , and l 6 are also predicted using two single-position identifications (i.e., l 3 and l 5 ) to confirm that the multi-position identification method reflects the dynamic characteristics of the feed drive system well. The results are presented in Tabl. 8 and Tab. 9, respectively.
In Tab. 8 and Tab. 4, the accuracy is higher in the multi-position identification than the l 3 -position identification.
In Tab. 9 and Tab. 4, the overall error is much larger in the l 5 -position identification than in the multi-position identification, and the maximum error of the first four orders is 7.69%. Although this error is acceptable in [25], the second order of natural frequencies increases with l, thus contradicting the regularity of the measured data presented in Fig. 10. Therefore, the parameters identification by l 5 -position cannot truly reflect the dynamic characteristics of the system.
The multi-position identification method considers the effect of the worktable position on the dynamic characteristics and utilizes the minimum overall error as the objective function. However, the single-position identification method disregards the effect of the worktable position, and the error of the natural frequency between the predicted and measured is large. Furthermore, the predicted value of the single-position identification may contradict the regularity of the measured data. The experiment shows that the multi-position identification method is robust and accurate for the joint parameter identification of the ball screw feed drive system, whose dynamic characteristic is affected by the worktable position. Table 8 Comparison between the predicted natural frequency using the l3-position identification and the measured natural frequency at Positions l2, l4, and l6