Dynamic modeling and experimental research on position-dependent behavior of twin ball screw feed system

To establish the dynamic model of machine tool structure is an important means to assess the performance of the machine tool structure during the cutting process. It is necessary to study the dynamics of the machine tools in different configurations for the sake of analyzing the dynamic behavior of the machine tools in the entire workspace. In this paper, a robust approach is presented to build an efficient and reliable dynamic model to evaluate the position-dependent dynamics of the twin ball screw (TBS) feed system. First, the TBS feed system is divided into several components, and a finite element (FE) model is built for each component. Second, the Craig-Bampton method is proposed to reduce the order of the substructures. Third, a multipoint constraints (MPCs) method was introduced to model the mechanical joint substructures of the TBS system, and the spring-damper element (SDE) is employed to connect the condensation nodes. Finally, a series experimental tests and full-order FE analysis are conducted on the self-designed TBS worktable in the four positions to validate the effectiveness of the proposed dynamic model. The results show that the proposed approach evaluates accurately the position-dependent behavior of the TBS system.


Introduction
TBS feed system is widely employed in gantry precision CNC machine tools with high rigidity, high stability, and high control bandwidth. Its dynamic properties make a big influence on the control performance machining accuracy and stability [1]. Therefore, the study of dynamic characteristics of TBS system has a significant effect on improving the feed accuracy and machining accuracy [2][3][4][5][6]. The machine tool structure exhibits a position-dependent dynamic behavior which affects the stability limit of the process lying on the worktable position. When considering a ball screw feed system, it should not be overlooked that the dynamics of the system depends on the position of the nut on the ball screw. Modeling the coupling points between the flexible substructures for different axial positions is challenging and can be dealt with different ways. Zaeh et al. [7] utilized the node model of ball screw and guide rail together with the function of parabolic load, to connect the two components. The other way round, Spescha et al. [8,9] employed triangular interpolation based on Fourier series expansion, to introduce loads between moving interface of the feed axes of the machine tools. Brecher et al. [10,11] proposed rigid-body elements (RBE) that move with a component and search for the nearest node at the desired position. Law [12,13] demonstrated Lagrange-multipliers method and penalty method to calculate the dynamic response of the reduced order model. Garitaonandia et al. [14] proposed the linear spring model to represent the flexible connections between two substructures of the machine tools.
In order to evaluate the position-dependent dynamic performance in the machine tool design stage, the FE analysis method for machine tool design and analyses has been certified effective for subcomponent level design and structure optimization. Nevertheless, dynamic response analyses of a full-order machine tool FE model which are usually millions or more degrees of freedom (DOFs) that are expensive and time-consuming to calculate and accounts for a large part of the total computational effort are required to design and analyze a machine tool [15]. In addition, the establishment of position-dependent in such large-order FE models requires complex adaptive strategies, which takes a lot of time in practice. In order to simulate the dynamic performance of a machine tool efficiently, the machine tool model needs to be partitioned into substructures, and the dynamic substructuring approach is also applied [16]. This allows each axis to be calculated separately and coupled to reduce the model in different axis positions without reducing the entire machine model in each axis position. There are many different methods to partitioning mechanical structures into substructures, of which the most popular are component mode synthesis (CMS) and frequency-based substructuring (FBS) approaches [17]. References [18][19][20] all proposed FBS method, and applied to machine tools, using a simple multipoint constraint coupling method to represent the connection between substructures with simulated models. A system with ball screws and linear guides has a large number of coupling points along an axis, and each coupling point requires a transfer function. In addition to FBS approach, subspace techniques, the Krylov subspace method based on subspace technology can also be used to reduce the order of machine tool structures [21]. However, the CMS approach, and in particular the Craig-Bampton reduction [22], is one of the most common techniques for dividing mechanical structures into substructures and reducing their orders and has usually been implemented in commercial FE analysis software. The order of the dynamic system matrix obtained is relatively low and can be handled pretty well in further simulation. Since the coupling points are part of the state vector, the system can be easily adjusted by leading in new stiffness or damping between the two coupling points [23,24]. Applied the Craig-Bampton reduction method to build a reduced state-space model of a machine tool, which was contained in a framework for optimally designed electromechanical systems [25]. Adopted the Craig-Bampton reduction method to minimize the DOFs of a machine tool.
In the assembly process of substructures, the mechanical joints between substructures also need to be taken into account. Due to the difficulty of modeling joints in a FE environment, it is still a challenge to accurately predict the dynamic performance of machine tools during the design stage. Joint performance lie on critical factors such as contact interface conditions, friction, stiffness, and damping, and since approximately 60% of the total dynamic stiffness and 90% of the total damping in the machine tool structure come from the joints [26]. Therefore, the contact joints must be considered in the model, and the detailed and accurate modeling of machine tool engagement characteristics is an independent and important research topic and is the subject of many hot research surveys [27,28]. For coupling two substructures, multipoint constraint (MPCs) can be applied to connect to the spring-damper element (SDE) system. In order to allow coupling of substructures in different axis positions, MPCs needs to be used after model order reduction. References [29][30][31] showed several methods to use the MPCs to a reduced order model of machine tools. Nevertheless, these substructuring methods do not consider the TBS feed system, which limits the accurate prediction of the position-dependent dynamic behavior of the machine tool equipped with the TBS system.
In this paper, a dynamic substructuring method with position flexible coupling of the substructures is presented, which considers mechanical joint models for each direction. Emphasis is placed on comparing different methods to predict the computational time and modeling accuracy of positiondependent dynamic behavior. First, the CMS method is presented to reduce the internal DoFs of the TBS system and greatly decrease the calculation time. Next, considering the mechanical joints among the substructures, MPCs model approach-based spring-damper element (SDE) is used to couple the flexible substructures. Third, for the sake of simulating the position-dependent behavior of the TBS system, Lagrange multipliers method is introduced. Furthermore, a full-order FE analysis in the ANSYS software is also added to compare with the proposed method. Finally, a series of modal experimental tests are conducted on the TBS system as depicted in Fig 1a. The results showed that the proposed dynamic substructuring method based on MPCs joint model can accurately and quickly predict the dynamic characteristics of TBS feed system.

Model order reduction of the TBS system
The TBS feed system analyzed in this study is the dual-drive CNC grinding wheel dressing machine tool as shown in Fig.  1a. Fig.1b shows the 3D model of the TBS feed system. As shown in the figure, the configuration mode of "dual motors + dual screw" is adopted in the X/Y feed direction to jointly drive the worktable, which has the characteristics of large thrust and good rigidity but is essentially a redundant drive. The ball screw adopts a double nut structure, which can effectively eliminate the axial clearance. Dual linear guide rails are symmetrically arranged inside the dual ball screws to realize the guiding and bearing function. A pair of angular contact ball bearings are used near the end of the motor, which plays a role in supporting and preventing the axial movement of the screw, and a deep groove ball used far away from the end of the motor only plays a supporting role. Linear grating ruler is installed on both sides of the worktable to detect and track the displacement of the worktable in real time and realize the closed-loop control of the feed system. The whole 2D worktable is fixed on the base by bolts. This structure is widely used in large gantry machine tools. Before developing the FE model, the CAD model of the TBS worktable should be simplified and then import to ANSYS software, as shown in Fig.2a. The TBS worktable is partitioned into 4 substructural components, as indicated in Fig.2b. FE modeling is carried out for each of the main substructural parts of the machine after partitioning. However, each substructure can still have over one million DOFs. In order to improve the calculation time of dynamic response analysis, each substructure was reduced using Craig-Bampton method [22], which has proved to be a practical and widely used CMS method implemented in commercial FE software such as FEM Tools.
Consider the substructure A, whose equation of motion is: whereM A , C A , and K A are the mass matrices, damping matrices, and stiffness matrices of the substructure A, respectively. u A ,u A , and :: u A are the vector matrices of displacements, velocities, and accelerations, respectively. The vector f A indicates as the forces acting on the substructure at the point of connection between the substructure and the adjacent substructure. External forces are not considered in this article, as they are independent of the dynamic analysis to be performed later. Though the different damping reasons can be researched and integrated in the FE model [16,19], the constant damping is considered in this paper with damping C = 0.02. The main source of damping is related to energy dissipation at the contact interface, and the modeling of the interfaces will be implemented in a later stage (Section 3).
Applying the Craig-Bampton method requires partitioning the system matrices into boundary (subscript b) and internal T is the transpose notation. Therefore, Eq. (1) becomes: Notice that since all the forces between the substructures are transferred through the boundary, the internal force f A i is zero. The Craig-Bampton transformation replaces the physical coordinates with a mixed set of coordinates, including physical and modal coordinates. If the full DOF is n, so n = n b + n i (n i is internal DoF, n b is boundary DoF). The boundary set u A b stays the same, while the internal set u A i is converted to a set of truncated modal coordinates q A m (m < n i ). The relationship between the physical coordinates u b u i f g AT and the mixed Craig-Bampton coordinates u b q m f g AT is as follows: The vectors in B A are called constraint modes, which are defined as structural displacements of physical coordinates caused by continuous unit displacements on the boundary DOFs, while keeping the other boundary constraints and int e r n a l D O F s f r e e , w h e r e B The vectors in Φ A F are usually specified as fixed-interface modes, which can be obtained from the constrained boundary set (u A b ¼ :: where Φ A i is the mode shape vectors, which are obtained by solving the eigenvalue problem corresponding to the internal DOFs of the equation In actual engineering, the mode is truncated into a simplified set, ignoring the contribution of the high frequency modes. In general, it is recommended that the substructure retain the modal content corresponding to the mode shape, with a natural frequency of 1.5 times the highest frequency of interest. Substituting the vectors B A and Φ A F into Eq. (3), the Craig-Bampton transformation is obtained: Applying this transformation to Equation (2), the resulting equation is pre-multiplied by the transpose of the matrix A, resulting in Suppose the mass of matrix Φ A i is normalized, without considering the internal forces (f A i ¼ 0 ). The matrices included in Equation (5) are those that define the Craig-Bampton formula and must be computed, where is the mass matrix reduced to the DOFs of the boundary, just the same as the Guyan transformation.
is the reduced stiffness matrix of Guyan, and ω 2 0 À Á A is the diagonal matrix including the square of the natural frequencies of the fixed-interface modes. Considering the substructures A and B and assembling the reduced order dynamic equation, the following can be obtained: By reducing the internal DOF of the substructure, CMS method reduces the higher order modes of the structure, while retaining the lower order modes which are more valuable to the machine tool design, thus greatly reducing the structural DOFs and greatly reducing the calculation time and improving the calculation efficiency on the premise of ensuring the accuracy. This formula is very common in structural analysis and can be implemented in most FE software, while previous matrices are extracted from ANSYS software. As shown in the Fig. 1c, all the internal DOFs of the TBS system is reduced, and the interface DOFs of the joints is left after assembly.

Interface modeling between the substructured TBS components
A standard machine tool contains many different kinds of contact interfaces, each with different performance that has a different impact on the dynamic response of the entire machine tool. The fixed interface composed of the bolted connections between structural members, the tool-holder interface, and the spindle interface. Here, the fixed joints of bolted connections are treated as groups. The moving interface include screw-nut interface and guide-slider interface. In this paper, a multipoint constraints (MPCs) approach based on SDE is proposed to model the interface.

Multipoint constraint method
Take the two substructures S A and S B as the research objects, shown in Fig.3a, which can represent the relative motion between any TBS components. Using the Craig-Bampton method discussed in Section 2, these substructures have been reduced to their interface DOFs. For the sake of realizing high efficiency and high precision error compensation, the dynamic model between substructures was established by introducing the MPCs model. In the TBS system, each linear axis contains a moving and a stationary substructure, which are usually connected by 2 ball screws and the 2 guide rails. This leads to 4 coupling points for each feed direction, 2 represent the joints of the screw nut, and 2 represent the joints of the guide rail slider, as shown in Fig.3b. To overcome the drawback of the incompatible meshes between the substructures, an approximate model of interface interactions is gained by defining a virtual condensation node at the center of each contact interface, as shown in Fig. 3a.
The contact nodes are coupled by the condensation nodes, which are connected by the SDE. Most of these elements are modeled as three-dimensional (3D) SDEs with nodes connected to contiguous substructures by MPCs, as shown in Fig.3b. The SDE consists of the stiffness and damping parameters of the coupling elements and is taken as part of the moving part.ΔL x and ΔL y as shown in Fig. 2a represent the distance from the right end of the X and Y direction, respectively, and change their values to change the position of the worktable.
Each of the four pairs of contact interfaces between the four ball screw groups and the guide group is represented by a pair of condensation nodes connected to the interface DOFs using MPCs formulations.
According to Eq. (7), the coordinate transformation consists of two items, namely the condensation node translation u A and the cross product of θ A and d i , where θ A is the rotating angular of the condensation node, and d i is the distance vector from the condensation node to the coupling nodes. Thus, the coordinate transformation matrix of each SDE connecting the coupling node and point to the condensation node can be gained by  If the coupling nodes do not have rotational DOFs, delete the last three rows from u A and T i . The rotational DoFs of coupling nodes are not considered in this paper. This modeling method can overcome the shortcoming of non-conforming mesh at joints, and for the sake of enhancing the approximate geometric coordination between substructures, Lagrange multiplier method is usually introduced.

Lagrange multiplier method
A famous approach to integrating MPCs into ordinary differential equations (ODEs) is to expand matrices through constraints. Using the discrete set of m Lagrange multipliers λ corresponding to m constraint equations, the modified potential functional form after the first variation is Ref. [32][33][34] The system matrices are extended by new DOFs λ, which can be interpreted as the resulting coupling forces when the MPCs are satisfied. Now the whole ODEs look like this Substitute Eq. (9) into Eq. (10), we get In Fig. 2a, the assembled motion equation guaranteeing the compatibility of the two substructures at the interface is þ The M and K matrices of substructures can be extracted from ANSYS software to obtain TXT file and then imported into MATLAB for modal comprehensive calculation. In this paper, the calculation method and value of the stiffness of the moving contact interfaces are referred to literature [35]. The computer CPU used in this paper is I7-3630QM, and the memory is 16G.
When the TBS worktable moves to a new position, by appropriately updating the displacement operator L in Eq. (10), a new set of constraint equations is developed to express the relationship between the new nodes in contact joints. Then, a new set of reduced synthesis constraint equations Eq. (12) is solved to gain the position-dependent dynamic response. By changing the position of substructure B set by ΔL x and ΔL y , as shown in Fig.3a, the dynamic response at any position is obtained. The proposed formulation presents a complete description of modeling the position-dependent dynamics based on defining constraint equations to synthesize reduced substructures. Based on the defined constraint equation to synthesize the simplified substructures, a complete modeling method of positiondependent dynamics is proposed, Fig. 4 shows the calculation flow of position-dependent dynamic characteristics of the TBS system.

FE analysis method
In the ANSYS environment, the FE analysis of the full DoF TBS feed system was carried out by referring to the joint stiffness in reference [35]. According to the real material of the parts, the material parameters are shown in Table 1. The first four natural frequencies at different positions can be gained by solving and analyzing the FE model. As shown in Table 2, because the stiffness of fixed interface was ignored in the FE analysis, the natural frequency obtained by FE analysis was greater than the natural frequency obtained by experiment when the worktable is at the same position. However, these will be improved in subsequent studies. In ANSYS, when the worktable is located at four positions of the stroke, Block Lanczos method is used for modal solution. Therefore, the CAD model of TBS workbench needs to be adjusted for four times and re-imported into ANSYS for simulation calculation. At the same time, the vibration modes of the two-dimensional TBS worktable are obtained.
The first four modes of the TBS table considering joint stiffness are shown in Fig.5. Fig.5 a and b show the firstorder and second-order translational vibration modes of the worktable in the Y and X directions, respectively. Fig.5c shows that the third vibration mode is high order yawing vibration of the worktable in the Z direction. Fig.5d shows that the fourth vibration mode is the pitching vibration of the worktable around Y direction.

Experimental test
To further validate the proposed method, a series of dynamic experiments were conducted on a self-designed TBS worktable depicted in Fig. 1a; the worktable is bolted to the base. X and Y feed directions are respectively driven by "dual servo motors + twin ball screws" to achieve the linkage of two feed axes, and the movements of the X-and Y-axis of the machine all cover 200 mm. On the basis of the transmission stiffness model, the transmission stiffness of the TBS feed system relies on the position of the worktable. Modal tests were carried out on the worktable at four typical positions (see Section 4.2) because it is important to determine the dynamic characteristics of the worktable at different positions for the actual machining process. The position of the worktable is controlled by the control system, as shown in Fig.6a. The TBS system is equipped with the Beckhoff motion control system and equipped with a grating ruler to control the position of the worktable through the human-computer interaction interface.
In this paper, the dynamical experiment of TBS system is carried out by hammering method. The dynamic characteristics of the worktable at different positions were tested by Econ dynamic analyzer. As shown in Fig.6b, four three-axis acceleration sensors (model: BK4507B) are set at the four corners where the vibration amplitude of the worktable is large. The method of single point excitation and multipoint vibration pickup is adopted, and the experiment is carried out by using the impact hammer (model: DYTRAN 5800B5). Econ MI7008 acquisition card is used to achieve high-precision data acquisition. Data processing was completed by Modal Genius software. The frequency response curve and mode shapes were obtained by data fitting.

Results and discussion
The motion position is defined by the coordinate (x, y) composed of the displacement of the working table in the X and Y axis. Four motion positions are selected in this paper, which are (0,100), (100,100), (200,100), and (200,200), respectively, as shown in Fig.7. The prediction in this paper concentrates on the frequency range of 0-500 Hz, including the relevant dynamic modes of the TBS system. These motion positions contain the minimum and maximum range working volume of the TBS worktable. The numerical simulations were made for four different motion positions of the worktable, and the comparison with experiments as well as the full-order FE analysis  indicated that the CMS method is able to identify the most representative modes in the input-output dynamic behavior and to evaluate the contribution of these modes in the response, which is shown in the Fig.7. By comparing the FRFs at different positions, it is not difficult to see that the dynamic performance of TBS system shows obvious position dependent; the peaks of the resonant frequencies are distorted differently. The established TBS system model has more than one million DOFs, and if a full-order FE model is used, it will take 5.38 h of calculation time for each desired axis position. By employing the proposed CMS method, the number of DoFs in the TBS system is reduced to 240, and the axis of motion can be coupled to different axis positions within a calculation time 78 s. It can be seen that using the CMS method proposed in this paper, the calculation time is shortened 248 times, and the calculation efficiency is very high.
The position-varying dynamic behavior is also tabulated in Table 2. With the change of position, the natural frequency of the TBS worktable changes significantly. Taking position 1 and position 2 as an example, when the worktable moves from position 1 to position 2, that is, it moves along the X direction to 100mm, the first-order natural frequency decreases from 53.4 to 51.7Hz and continues to decrease to 49.2Hz when it moves to position 3. This is because the ball screw of TBS worktable is installed with one end fixed and the other end free. When the worktable moves along the free end in the X direction, the overall equivalent stiffness of the system decreases, resulting in a decrease in the natural frequency. Because the TBS worktable travel is too short, the overall size is small, resulting in the dynamic parameters change that is not obvious. The natural frequencies obtained by experiments are all lower than those obtained by the other two mathematical methods, which is caused by the fact that the stiffness of a large number of fixed joints is neglected in the calculation process, thus leading to the increase of the overall equivalent stiffness. As shown in Fig. 8, the position-dependent nature frequency response comparisons between the reduced order CMS models and full order FE models. The first three nature frequencies of the CMS method are almost the same as those of the full-order FE method. With the increase of order, the truncated error is caused by order reduction. The normalized relative frequency difference (NRFD) is defined to determine the accuracy of the proposed reduced order model over the full-order FE method. NRFD = | 1 − f Full /f Rom | × 100%. As is shown in Table 2, the range of NRFD is 1-9.2%, indicating that the order reduction method has good accuracy, and the error generated by the CMS method for order reduction is all within 10%. Moreover, the accuracy of CMS method with MPCs constraint modeling is higher than that of FE method, which proves the effectiveness of the proposed method for the modeling process and validates the obtained model. In addition, nonlinearities in the contact joints are not currently considered but may significantly affect global mechanical dynamics. Fig. 6 Modal tests of the TBS feed system. a The control system of the TBS worktable. b The dynamic experiment of the TBS system. In this paper, a robust and computationally efficient method is proposed to model and evaluate the position-dependent dynamic characteristics of a self-designed TBS system considering contact nodes between substructures. In order to realize the effective simulation, a dynamic reduced model of TBS system is established by using the Craig-Bampton method. The contact joints effects are taken into account within the reduced model, and the MPCs model with SDE is introduced to describe the connections, which are convenient to model ball screw and linear guide joints during the worktable movements. The machine substructure components are reduced and then assembled using MPCs equations based on condensation nodes, overcoming mesh incompatibles on the contact interfaces and allowing a modular design.
The effectiveness of the proposed model reduction method is verified by comparing the simulation and measurement of the dynamic characteristics of the TBS workbench in four typical motion positions. On the one hand, the full-order FE analysis was conducted in the ANSYS software. On the other hand, a series of experimental tests were operated on the TBS worktable. By comparison, the proposed substructure method has higher computational speed and precision. The proposed method allows machine tool manufacturers to use high-precision detailed models when developing new machine tools, while reducing computation time when design changes or axisposition changes occur. However, some obstacles such as the modeling of position-dependent coupling elements and nonlinear effects need to be addressed in the future to enable an effective optimization approach.