A comprehensive nonlinear dynamic model for ball screw feed system with rolling joint characteristics

Modern tendency of machine tools design requires more accurate model to predict the system dynamics, in order to anticipate its interaction with machining process. In this paper, a comprehensive dynamic model of ball screw feed system (BSFS) considering nonlinear kinematic joints is introduced to investigate the varying dynamic characteristics when worktable is subjected to combined load from six directions. The load–deformation relationship of each kinematic joint is dealt with a set of translational and angular spring elements. The nonlinear restoring force function of each joint involving coupling displacement is calculated. Based on the lumped mass method, the analytical 18-DOF dynamic equation is formulated by the analysis of the interaction force between joints. Model verification tests are conducted. The worktable response exhibits the abundant and fascinating nonlinear phenomena arising in nonlinear joint and coupling effect. The nonlinear behavior behaves significant difference owing to the variations of excitation, platform position, screw length, preload and damping of joints. Thus, the model is promising for comprehension of machine dynamic behavior and for development of sophisticated control strategy.


Introduction
With the ever-increasing development of modern manufacturing industry, demands for more precise and high-speed CNC (computer numerical control) machine tools have been put forward. The quality and productivity of CNC machine tools are often dominated by the static and dynamic performances of feed drive system [1][2][3][4]. Further, vibration during processing is the essential factor that deteriorates surface finish and dimensional accuracy of components, such as the inferior product machining of thin-walled parts, integral impellers and high-precision crankshafts [5][6][7]. An accurate dynamic model of feed drive system is essential for the prediction of structural vibration and the enhancement of operating behavior of CNC machine tools. The BSFS accompanied with rolling linear guideway platform is the popular unit to realize linear motion of CNC machine tools owing to good kinematic accuracy, smooth motions, low wear and long service lifetime without a stick-slip effect [8]. The approximately 60% of the total dynamic stiffness and approximately 90% of the total damping of the entire machine tool structure originates in the kinematic joints [9,10], and then the complication and instability of vibration may be dramatically affected by the nonlinearity of joints among supported bearing, nut and carriage [5,[11][12][13][14]. On the other hand, more accurate model considering structural nonlinearity is required by the virtual prototype and virtual machining [15]. Therefore, inclusion of the ball point contact of joints allows the present nonlinear vibration analysis to ensure controllable dynamic behaviors of BSFS.
Numerous scientific works have dealt with the analytical mechanical model of BSFS, with the aim of replicating the dynamical behavior as detailed as possible and improve the translational and positioning accuracy. Vicente et al. [16] developed a dynamic model for the first three resonant mode shapes analysis of BSFS. The axial-torsional coupling effect was investigated when screw shaft served as a continuous subsystem using a field function. However, the model was proposed in the absence of experimental verification. Similarly, Dong and Tang [17] raised a mathematical model describing the continuous deformation of flexible screw using a field function with unknown bounds. Their model also contained timevarying parametric uncertainties and disturbances to investigate the tracking control for a high-speed BSFS. Based on the Timoshenko beam assumption, Wu et al. [18] modeled the lead screw as a continuous body when the other parts of BSFS were treated as lumped mass units. The high-order modes of system and flexible feature of screw were characterized by the proposed hybrid model. Instead, Frey et al. [19] described the screw as a simple lumped mass and identified the ability of proposed model on expressing the most relevant characteristics of feed drives with respect to their application in machine tools. Sato [20] also reported a simple and physically intuitive 2-DOF model for expressing the first-order mode of axial and rotational vibrations. It means that the axial stiffness of BSFS can be estimated as a serial connection of stiffness of the support bearing, screw shaft and nut. Ansoategui and Campa [21] created a detailed 7-DOF lumped parameters dynamic model and compared with a classical 2-DOF model. The advantage of proposed model was that the 7-DOF model modal analysis could identify the weaker components of transmission during a redesign of the drive or transmission chain. Guo et al. [22] performed a novel idea for the hybrid dynamic model of BSFS through combining classical dynamics theory and data driving. Nguyen et al. [23] constructed a discrete dynamic model to investigate the interaction between preload, the displacement of working table and the axial natural frequency of screw nut.
The nonlinearity of joints may result in substantial variation in natural frequencies and frequency responses of machine tool, which are rarely mentioned in antecedent literature. Considering the nonlinear stiffness of rolling guide joints owing to the point contact between ball and raceway, Wang et al. [24] developed a dynamic model of three-axis coupling machine tool to evaluate the axis coupling effect as well as position-dependent response of tool center point. On the basis of nonlinear contact characteristics of screw-nut and bearing joints, for vertical structure BSFS without counterweight [25], the slender BSFS [26] and BSFS with lead screw pre-stretching [27] and the position-dependent transmission stiffness models of system were established, respectively. Nevertheless, the abovementioned studies pay little attention to nonlinear vibration behavior of BSFS which is necessary to evaluate the operational behavior of CNC machine tools under periodic cutting force treatment. Li et al. [28] and Xu et al. [29] developed a timedependent dynamic model for clarifying the micro posture fluctuation mechanism of ball screw and carriage, respectively. Due to the ball's point contact property, the inherent frequency multiplication and demultiplication components were captured when system was excited by the number variation of inserted loaded ball. Gu and Zhang [30] and Xu et al. [31] investigated the resonance state of BSFS using maximum amplitude-frequency curves. The result manifested the remarkable nonlinear phenomena, like super-and subharmonic resonance and jumping discontinuous behavior. Kong et al. [32] demonstrated the existence of unstable chaotic motion in linear guide through bifurcation diagram and maximum Lyapunov exponent. It is worth mentioning that these models ignored the coupling effect when the BSFS suffers from multidirectional loads, and were directed to the unidirectional vibration characteristics. Liu et al. [33] formulated a nonlinear dynamic model considering the assembly error of deflection angle caused by the misalignment between front and rear bearings. But the model discarded the guide platform which acts as the dominant carrying parts and the effect of torque was also not included. Although Wang et al. [13] derived the coupling relevance of linear guideways, screw nut, screw shaft and supported bearings and studied the mechanical behavior of nonlinear kinematic joints for ball screw-driven CNC axis, the yawing vibration behavior of worktable failed to be estimated and only the axial vibration response of screw nut, shaft and bearings was highlighted. From the current research situation, most of the mechanical model for BSFS is limited to the axial dynamic characteristics, without involving the coupling effect of multidirectional loadings, which is incomplete for the accurate prediction of dynamic response of BSFS. This paper will try to establish a more comprehensive dynamic model to expand the research on structural dynamics of CNC machine tools.
In this paper, the combined load from vertical, lateral, feed, yawing, pitching and rolling directions is applied on the worktable of BSFS. In conventional machine tools, due to intensive FE optimization of machine structure in design, the structural parts are much stiffer than joints. The structural modal is at a relatively high frequency, and the predominant modes in the low-frequency and mid-frequency are entirely caused by limited joint stiffness. Thus, the dynamics of system is mainly dependent on the stiffness of joint rather than structures [24]. It is appropriate to model BSFS by lumped mass elements and a set of translational and angular spring elements for each kinematic joint. Aiming to describe the nonlinear contact of kinematic joint, the restoring force function of each supported part with respect to translation and angular displacements is calculated, involving two translations and three deflections of linear guideway platform, three translations and three deflections of screw nut, three translations and two deflections of matched angular contact ball bearings, as well as two radial translations of supported deep groove ball bearing. Additionally, the screw shaft is treated as a simple lumped mass element and its three translations and two deflections (excluding the torsion deformation) are considered; furthermore, the position-dependent restoring force expression of lumped screw shaft unit is derived by the classical deformation formula of loaded beam. The approach aims to identify the dominant effects of screw shaft and express the most relevant position-dependent characteristics of BSFS. With the aid of lumped mass method, the analytical 18-DOF equation of motion with coupling displacement and nonlinear restoring force term is introduced. To the best of authors' knowledge, the proposed detailed model is not shared by any former author. Moreover, experiment is also arranged for the validations of calculation methodology and vibration response through a BSFS test bench. Nonlinear vibration analysis based on numerical solution is carried out according to amplitude-frequency curve, bifurcation diagram, 3-D frequency spectrum, velocity-displacement phase diagram, Poincaré section and frequency domain. Effect of external harmonic excitation, position-dependent, length of screw shaft, preload and damping of joint on dynamic response is discussed comparatively. The proposed model can be used to identify the mechanism of time domain vibration response, understand the posture variation, and develop precise numerical model algorithms to determine the vibration characteristics of BSFS.

Dynamic model
The most common BSFS for precision positioning consists of screw shaft, which is assembled to the machine pedestal by rotary bearings, and is driven by a servomotor through a flexible coupling, as shown in Fig. 1. The screw-nut joint is attached to worktable that is constrained to move axially on carriage and rail. For the fixed-supported installation of Fig. 1, the fixed end close to the motor is performed by a pair of angular contact ball bearings, and a deep groove ball bearing achieves the radial support of the other end. The right-hand Cartesian coordinate XYZ-O is defined, and the origin O lies on the symmetry center of top surface of worktable. The external loads F x0 , F y0 , F z0 , M x0 , M y0 , M z0 from six directions are applied on the platform, and the derived translation and angular displacements of platform center are x, y, z, h x , h y , h z corresponding to vertical, lateral, feed, yawing, pitching and rolling directions, respectively.

Equation of motion
It is assumed that screw is suitably mounted in the servomechanism and non-concentric force produced by misalignments is discarded. The stiffness of servomotor and torsion deflection of flexible screw shaft are not considered. The rolling joints deformation is supposed to only occur in the contact area between ball and raceway. Figure 2 shows the equivalent dynamic model of BSFS using the lumped mass method. The multi-body dynamical system is formulated by four lumped mass points, and it refers to linear guideway platform which consists of four carriages, worktable and nut; flexible screw shaft; matched bearings; and radial bearing. The rolling joints and flexible screw shaft are regarded as the spring-damper parts. The displacement vector is defined as [x, y, z, h x , h y , h z , x ss , y ss , z ss , h ssx , h ssy , x mb , y mb , z mb , h mbx , h mby , x db , y db ], as shown in Fig. 2. Then, the 18-DOF differential equation of motion with time-varying restoring force term under harmonic excitation can be written as follows: For linear guideway platform: where m p and I px , I py , I pz denote the mass and moment of inertia of platform (including table, carriage blocks, nut bracket and nut); . and .. are the first and second derivatives with respect to time t; x is excitation angular frequency; x ss , y ss , z ss , h ssx , h ssy represent translation and angular displacements of lumped screw shaft unit; F lgx , F lgy , M lgx , M lgy , M lgz are the synthetic restoring force and moment of four carriages; F snx , F sny , F snz , M snx , M sny , M snz are restoring force and moment of nut joint; c lgx , c lgy , c lghx , c lghy , c lghz and c snx , c sny , c snz , c snhx , c snhy , c snhz are the corresponding viscous damping coefficients. The approximated same damping for all ball contacts within the same joint is presumed. In reality, damping characteristic of elastic contact parts is nonlinear as well. In order to focus on the study of nonlinear elastic restoring force, the simplified constant proportional viscous damping coefficients are adopted for analyzing dynamic response owing to the fact of relatively small damping by analogy with rolling element bearing [34][35][36][37]. Because screw shaft is fixed on the base by bearings at both ends, the relationship between absolute displacement and deformation of lumped screw shaft unit should be clearly given before deriving its equation of motion. Figure 3 describes the geometric relationship of screw displacement, deformation and bearing displacement. The translation and angular deformations of screw can be obtained: At the XOZ plane at the YOZ plane where L is the length of screw shaft (interval between matched and radial bearings) and L p is the distance between platform and matched bearings; x mb , y mb and x db , y db are the radial displacements of matched and radial bearings, respectively; h mbx and h mby are the angular displacements of matched bearing around vertical and lateral directions, respectively; w 1 and w 2 are the auxiliary parameter; x ss , y ss and h ssx , h ssy are the absolute translation and angular displacements of screw; w ssx , w ssy and # ssx , # ssy are the corresponding translation and angular deformations. Further, the equation of motion for screw can be formulated: where m ss and I ssx , I ssy denote the mass and moment of inertia of screw; F ssx , F ssy , F ssz , M ssx , M ssy are the restoring force and moment of screw; c ssx , c ssy , c ssz , c ss0x , c ss0y are the corresponding damping coefficients; z mb is the displacement of matched bearing joints in feed direction. For matched bearings, the equation of motion is expressed: Fig. 3 The relationship between absolute displacement and deformation of lumped screw shaft unit where m mb and I mbx , I mby denote the mass and moment of inertia of matched bearings (including inner and its shaft segment); x mb , y mb , z mb , h mbx , h mby are the displacements of matched bearings joint in five vibration directions; F mbx , F mby , F mbz , M mbx , M mby are the restoring force and moment; F ssmbx , F ssmby , M ssmbx , M ssmby are the interaction force that exerts on matched bearings by screw; the specific expression is given in Sect. 2.4; c mbx , c mby , c mbz , c mbhx , c mbhy , c ssmbx , c ssmby , c ssmbhx , c ssmbhy are the corresponding damping coefficients. For radial bearings, where m db denotes the mass of deep groove ball bearing (including inner and its shaft segment); x db , y db are the radial displacements; F dbx , F dby are restoring force; F ssdbx , F ssdby are the interaction force that exerts on radial bearing by screw, the specific expression is also given in Sect. 2.4; c dbx , c dby , c ssdbx , c ssdby are the corresponding damping coefficients.
2.2 Calculation for elastic restoring force of linear guideway platform supported four carriages Figure 4 shows the geometric information of platform and loaded ball arrangement between carriage and rail raceways. There are four rows of ball (i = 1, 2, 3, 4) with face-to-face configuration inside one carriage. Each ball groove has a contact profile of circular arc forming a two-point contact mode. The oversized ball within groove enables a negative clearance d lg0 between raceway and ball that is necessary for high stiffness and stability of system. The groove curvature radii of carriage and rail are r lg . The nominal diameter of ball is d lg . The initial distance O ci O ri between groove curvature centers of carriage and rail is A lg0 can be expressed: On the basis of Hertz contact theory, the normal preload Q lg0 of individual ball along contact angle a lg0 is [38,39]: where K lg is constant depending on material and geometric properties of contact area.

The relationship between ball deformation and platform displacement
As shown in Fig. 4b, the identical loaded ball number n lg and ball arrangement in each row are assumed and loaded balls are distributed with equal intervals. A straight raceway of carriage with defect-free is supposed and the crowning design which is used to reduce ball passage vibration is not considered. According to the dimensions of carriage and platform, ball location along feed direction can be given: where j(j = 1, 2, …n lg ) denotes the ball serial number; 2l p is the distance between carriage I and II or III and IV; l 0 is the length of carriage block. Then, the deformation of each ball can be obtained by the platform displacement and geometric relationship. For carriage I and II (j = I or II), the ball deformation along X, Y, Z-axes can be calculated: where d lgjijx , d lgjijy and d lgjijz represent the jth ball elastic deformation of the ith row along vertical, lateral and feed directions, respectively. As shown in Fig. 4a, 2wp is the distance between carriage I and IV or II and III; 2w 0 is the distance between 1 and 2 rows of ball; h is the distance between carriage center and top surface of platform and h = 0 denotes that the Xaxis passes through carriage center; h 0 is the distance between 2 and 3 rows of ball. Similarly, for carriage III and IV (j = III or IV), the ball deformation along X, Y, Z-axes can be calculated: Then, the new curvature center distance A lgjij of i = 1, 2, 3, 4 row of the jth ball under load can be obtained by: For carriage I and II (j = I or II), For carriage III and IV (j = III or IV), The normal deformation d lgjij of each ball under combined load can be described by: According to Hertz contact theory, the normal load Q lgjij of individual ball can be determined by: It is pointed out that when d lgjij \ 0, the ball will lose contact with raceway surface and d lgjij is set to Fig. 4 Ball arrangement within carriage; a ball contact geometry with face-to-face configuration; b ball location along feed direction equal 0. The spatial normal contact load Q lgjij can be decomposed into X, Y, Z-axes directions: where F lgjijx , F lgjijy and F lgjijz denote the contact load of individual ball along vertical, lateral and feed directions, respectively; a lgjij is the angle between Yaxis and the projection of new curvature center distance A lgjij on plane XOY; b lgjij is the angle between new curvature center distance A lgjij and plane XOY; and they can be formulated by: For carriage I and II (j = I or II), For carriage III and IV (j = III or IV),

The synthetic restoring force function by four carriages
In mathematical modeling, carriage and rail are connected with a series of spring elements with adequate stiffness. Based on a superposition of each ball load, the restoring force and moment of each carriage can be calculated by: For carriage I and II (j = I or II), For carriage III and IV (j = III or IV), where F lgjx , F lgjy , F lgjz , M lgjx , M lgjy and M lgjz (j = I, II, III, IV) denote the restoring force and moment of carriage j. The synthetic force and moment can be obtained by a superposition of each carriage load. The expression can be given by:

Calculation for elastic restoring force of nut joint
The offset-type preload screw-nut joint used in this paper has the benefit of a compact single nut with high stiffness via small preload force for light to medium loading applications in contrast with the traditional double nut preload method. The preload of ball is adjusted by a r value offset on the center pitch P sn of nut on the two ball tracks as depicted in Fig. 5a. The negative clearance d sn0 between raceway and ball can be achieved by the offset. The groove curvature radii of screw and nut are r sn . The nominal diameter of ball is d sn . The initial groove curvature center distance O sl O nl or O sr O nr of left or right-side screw nut is A sn0 and can be expressed: It is supposed that nut and screw in the nut length are rigid. On the basis of Hertz contact theory, the normal preload Q sn0 of individual ball along contact angle a sn0 is [40]: where K sn is constant depending on material and geometric properties of contact area.

The relationship between ball deformation within nut joint and platform displacement
The loaded ball number n sn and ball arrangement in left and right-side nut are identical. The ideal screw loop raceway is assumed. As shown in Fig. 5b, ball location along circumferential and axial direction can be given: The radial position angle / ksn of the k sn th ball (measured from ? X-axis): where n t is ball's number of turns in one-side nut (n t = 2.5 in this paper) and k sn (k sn = 1, 2, …n sn ) denotes the ball serial number, the axial position p k of the k sn th ball: Then, the deformation of each ball under load can be obtained by platform and screw superimposed displacements and geometric relationship. For left and right-side nut, the k sn th ball deformation along radial(r-axis) and axial(z-axis) direction can be calculated by analogy with matched bearings [41]: where d snlzk , d snlrk and d snrzk , d snrrk denote the k sn th ball deformation along z-axis and r-axis in the left and right nuts, respectively; the pitch circle radius is R psn ; H represents vertical distance between origin located at the center of platform and screw axis, as shown in Fig. 2.
Then, the new curvature center distance of the k sn th ball under load is: where A snlk and A snrk denote the k sn th ball curvature center distance in the left and right nuts, respectively; the lead angle is u and can be expressed: The normal deformation of each ball under combined load can be described by: where d snlk and d snrk denote the k sn th ball normal deformation in the left and right nuts. According to Hertz contact theory, the normal load of individual ball can be determined by: where Q snlk and Q snrk denote the k sn th ball normal load. It is pointed out that when d snlk or d snrk \ 0, the ball will lose contact with raceway surface and d snlk or d snrk is set to equal 0. The spatial normal contact load Q snlk and Q snrk can be decomposed into X, Y, Z-axes directions: where F snlxk , F snlyk , F snlzk and F snrxk , F snryk , F snrzk denote the load of k sn th ball in left and right nut along vertical, lateral and feed directions, respectively; a snlk and a snrk are the new contact angles and it can be formulated by: Fig. 5 Ball arrangement within screw-nut joint; a ball contact geometry; b ball location tan a snlk ¼ d snlzk cos u À A sn0 sin a sn0 A sn0 cos a sn0 þ d snlrk ; tan a snrk 2.3.2 The restoring force function of screw-nut joint Based on a superposition of each ball load, the restoring force and moment of screw-nut joint can be calculated by: 2.4 Calculation for position-dependent elastic restoring force of lumped screw shaft unit In this section, a simple approximation using lumped mass method is employed to model the screw shaft and this approach has been performed in Ref. [11,19,23,24]. The axial deformation d ssz of screw shaft may be expressed as: where E¯and A are the elastic modulus and crosssectional area of screw; the round rod with uniform cross section of diameter d is used to model the screw. Then, the axial restoring force function of screw can be written: Figure 6 shows the bending and deflection of screws under combined load and bearing constraints. The bending deformation at the radial bearing end is zero. At the XOZ plane, on the basis of deformation formula of classical loaded beam, the following equation can be obtained: where F ssdbx is the supporting force of XOZ plane at radial bearing end; w Fssx , w Mssy and w Fssdbx are deformation of radial bearing end along X-axis under F ssx , M ssy and F ssdbx , respectively; I is equivalent area moment of inertia. The constraint can be derived: where F ssmbx and M ssmby are the supporting force and moment of XOZ plane at fixed matched bearings end. Further, the translational w ssx and angular # ssy deformation of lumped mass point of screw can be calculated: Then, the time-dependent restoring force function F ssx and M ssy can be expressed: where Moreover, at the YOZ plane, the following equation can be obtained: where F ssdby is the supporting force of YOZ plane at radial bearing end; w Fssy , w Mssx and w Fssdby are deformation of radial bearing end along Y-axis under F ssy , M ssx and F ssdby , respectively. The constraint can be derived: where F ssmby and M ssmbx are the supporting force and moment of YOZ plane at fixed matched bearings end. Further, the translational w ssy and angular # ssx deformation can be calculated: Then, the restoring force function F ssy and M ssx can be described: where

The elastic restoring force of matched angular contact ball bearings
The matched bearings with back-to-back configuration achieve the support of fixed end of screw. The outer ring is fixed to the base. The preload with negative clearance d mb0 is applied to eliminate radial clearance, as depicted in Fig. 7a. The groove curvature radii of inner and out rings are r mb . The ball diameter is d mb . The initial groove curvature center distance O il O ol or O ir O or of left or right bearing is A mb0 and can be expressed: On the basis of Hertz contact theory, the normal preload Q mb0 of individual ball along contact angle a mb0 is where K mb is constant depending on material and geometric properties of contact area.
As shown in Fig. 7b, ball location along circumferential direction can be given: where / kmb is the radial position angle of the k sn th ball (measured from ? X-axis); n mb is ball's number of one bearing; k mb (k mb = 1, 2, …n mb ) denotes the ball serial number. Then, the deformation of each ball under load can be obtained by inner ring displacement (x mb , y mb , z mb , h mbx , h mby ) and geometric relationship. For left and right bearing, the k mb th ball deformation along radial(r-axis) and axial direction can be calculated: where d mblzk , d mblrk and d mbrzk , d mbrrk denote the k mb th ball deformation along z-axis and r-axis in the left and right bearings, respectively; the pitch circle radius is R pmb ; 2w b represents width of bearing. Then, the new curvature center distance of the k mb th ball under load is: where A mblk and A mbrk denote the k mb th ball curvature center distance in the left and right bearings, respectively. The normal deformation under combined load can be described by: where d mblk and d mbrk denote the k mb th ball normal deformation in the left and right bearings, respectively. According to Hertz contact theory, the normal load of individual ball can be determined by: It is pointed out that when d mblk or d mbrk \ 0, the ball will lose contact with raceway surface and d mblk or d mbrk is set to equal 0. The spatial normal load Q mblk and Q mbrk can be decomposed into X, Y, Z-axes directions: ; where F mblxk , F mblyk , F mblzk and F mbrxk , F mbryk , F mbrzk denote the contact load of k mb th ball in the left and right bearing along vertical, lateral and feed directions, respectively; a mblk and a mbrk are the new contact angles and it can be formulated by: Fig. 7 Ball arrangement within matched bearing joint; a ball contact geometry; b ball location tan a mblk ¼ Based on a superposition of each ball load, the restoring force and moment of matched bearings joint can be calculated by:

The elastic restoring force of deep groove ball bearing
For the deep groove ball bearing at the free end of screw shaft, deflection and axial deformations of ball are ignored. Similar to matched bearings, the ball location along circumferential direction can be given: where / kdb is the radial position angle of the k db th ball; n db is ball's number of one bearing, k db (k db = 1, 2, …n db ) denotes the ball serial number. Then, the deformation of each ball under load can be obtained by inner ring displacement (x db , y db ) and geometric relationship. The k db th ball deformation along radial direction can be calculated: where d dbrk denotes the k db th ball deformation along radial direction; c is radial clearance. On the basis of Hertz contact theory, the normal load Q dbrk of the k db th ball is: where K db is constant depending on material and geometric properties of contact area. When d dbrk \ 0, the ball will lose contact with raceway surface and d dbrk is set to equal 0. Q dbrk can be decomposed into X, Y-axes directions: where F dbxk and F dbyk denote the contact load of individual ball along vertical and lateral directions, respectively. Based on a superposition of each ball load, the restoring force of deep groove ball bearing joint can be calculated by: 3 Result and discussion

Model verification
The model parameters of studied ball screw feed system are given in Table 1. The above differential equations of proposed model with 18-DOF are solved by the fourth-order Runge-Kutta algorithm. The initial displacement and velocity of mechanical system are chosen as zero. The error for numerical approximation caused by time discretization matters the convergence of method. Thus, the termination criterion is specified to an error tolerance of less than 10 -7 to satisfy the requirements for convergent solution. A linear viscous damping is adopted in this study, and the damping of bearing is about 0.25-2.5 9 10 5 times of linearized stiffness of bearing [42], and the damping coefficient is determined within this range. Figure 8a shows the excitation test arrangement. The axial excitation position is selected for the response validation due to the coupled pitching and feed vibrations. The harmonic excitation signal is generated by shaker (type: Sinocera JZK-50) and  Hertz coefficient K lg (N/m 3/2 ) 9.14849 9 10 9 Hertz coefficient K sn (N/m 3/2 ) 8.07171 9 10 9 Hertz coefficient K mb and K db (N/m 3/2 ) 8.53657 9 10 9 frequency f shows a high conformity with measured one while there are some distinction between amplitudes of 2f and 3f. Some differences exist between the vibration amplitudes of the predicted and measured results, but the change law of predicted vibration is considerably similar to that of measured one when various excitation is applied. With the aid of the impact test arrangement in Fig. 8b, the FRF (frequency response function) in axial direction of system is extracted. The impact testing instruments include impact hammer (type: Sinocera LC-04A) and piezoelectric accelerometer (type : Brüel&Kjaer 4524-B-001). The force impacting in axial direction from the hammer is used as input, and acceleration is used as output to calculate the FRF. The raw impact signal is collected by data acquisition system (DH 5956), and the results of repeated five times are displayed in Fig. 11a, b when platform is located at 160 mm and 480 mm, respectively. The detected axial natural frequency varies with the position of screw-nut, but the results of multiple hammering tests at one position are the same. The nine positions in the 0-800 mm with the interval of 80 mm are selected to obtain various main peak frequency, and the result is given in Fig. 11c. It indicates the  position-dependent effect of axial natural frequency of system through experiment and proposed model, and the relative errors between measurement and simulation are 0.5%, 4.6%, 0.7%, 3.9%, 4.5%, 0.8%, 7.7%, 3.4%, 0.1%. The measured result agrees well with predicted axial natural frequency, and the maximum deviation does not exceed 8%. When the platform is far away from the fixed end of matched bearings, both  Fig. 11 The results of impact test; a-b FRF when platform is located at 160 mm and 480 mm; c axial natural frequency comparison between measured and simulated results the measured and predicted results indicate that the axial natural frequency shows a downward trend. When the platform is close to the radial bearing end, the axial natural frequency increases slightly. The comparison analysis indicates that the proposed model has an ability to predict the dynamic responses of ball screw feed system.
The axial sweep frequency experiment is carried out when platform is fixed at L p = 450 mm. Figure 12a shows the arrangement of equipment. The sine signal is produced by signal generator (Sinocera: YE1311) and delivered to the shaker through power amplifier (Sinocera: YE5874A). The scan frequency range is set at (50 Hz, 300 Hz), and the axial response in Fig. 12b is collected by the accelerometer. The eleven excitation amplitudes 0.2 g-1.2 g with the step of 0.1 g are employed. The experimental resonance frequencies are 145.6 Hz, 144.5 Hz, 142.9 Hz, 140.9 Hz, 139.2 Hz, 136.1 Hz, 134.7 Hz, 131.3 Hz, 129.7 Hz, 129.3 Hz, 128.4 Hz, respectively. The maximum deviation of resonance frequency is 11.85%, and the deviation cannot be found in the linear system. The phenomenon that resonance frequency point shifts with the change of excitation amplitude demonstrates that the nonlinear behavior exists in the ball screw feed system. The predicted Fig. 12 Sweep frequency experiment; a sweep frequency experiment site; b waves form when excitation amplitude is 0.6 g; c frequency response curve with various excitation amplitude; d comparison of resonance frequency points between theory and experiment result in Fig. 12d agrees well with the experimental value, and it proves the effectiveness of proposed nonlinear dynamic model.

Effect of external harmonic excitation on dynamic response
The response of worktable is output for simulation analysis. Figure 13 shows the frequency response  : 28kN), this value can be achieved in rare cases in engineering. For the purpose of theoretical analysis and scientific research, the simulation analysis under extreme conditions is carried out. The proposed general model and simulation within a wide range of parameters could provide the theoretical basis for dynamic structural optimization design and expand the application range of ball screw feed system, such as mining machinery, heavy-duty equipment, aerospace, high-speed and ultra-high-speed systems, and so on. The ordinate value of Fig. 13 is generated by peak-to-peak value of steady-state solution. The system undergoes resonances in a wide frequency range. As excitation frequency increases, the multiple resonance zones of nonlinear system are excited, and super-and subharmonic solutions and jumping phenomenon as the manifestation of nonlinear system occur. The harden effect is illustrated when jumping point bends to the side of increased frequency [43]. Figure 13c, e and f also captures the soften effect when jumping point bends to the side of decreased frequency corresponding to E = 1000. Taking E = 3000 for example, there is obvious distinction among the main peaks of frequency sweep plot in six directions. For Fig. 13a-f, the main peak frequencies are located at 1470 Hz, 1395 Hz, 135 Hz, 1320 Hz, 2160 Hz, 1860 Hz, respectively. The lowest and highest resonance frequencies are in feed and pitching directions. In Fig. 13c, the peak-to-peak value of main peak in feed direction is higher than vertical and lateral direction due to the fact that axial stiffness of BSFS is the weakest part. Besides, the frequency corresponding to main peak points of six vibration directions can be found in each frequency sweep plot. By the compared results of E = 1000, 2000 and 3000, when excitation amplitude is larger, vibration amplitude increases and jumping point corresponding to higher frequency is obtained. From Fig. 13c, e and f, the soften effect turns to harden effect at the certain frequency range when larger excitation amplitude is used. With the aid of the frequency sweep plot, the distribution of resonance frequency region can be observed. The corresponding frequency range of unstable motion is also captured when excitation frequency is used as the control parameter. In the design stage of ball screw feed system, in order to obtain a smaller vibration amplitude and dynamic stress, appropriate configuration and selection of structural parameters should be adopted to avoid the resonance frequency range and unstable motion zone related to severe vibration. The bifurcation diagram is employed to investigate the stability of BSFS. Although the result of frequency sweep plot exhibits obvious differences in resonance frequency of six vibration directions, bifurcation diagram in Fig. 14 captures the same evolution rule of dynamic behavior in six directions when the control parameter is excitation frequency and E = 3000, 5000 and 10,000, respectively. In Fig. 14c, the main peak frequency in the feed direction is also implied. The jumping discontinuous behavior is displayed, especially for the larger excitation amplitude of E = 10,000. The non-single-period motion, such as period-n and quasi-period, indicates that BSFS is in the unstable vibration state. Furthermore, the frequency range corresponding to unstable region expands and moves to the higher frequency when larger E is used.
The response of feed direction is chosen to conduct the next analysis. Figure 15 shows the bifurcation diagram and corresponding 3-D frequency spectrum with excitation amplitude E as the control parameter. The four typical frequency points (f = 1500 Hz, 2000 Hz, 2370 Hz, 2500 Hz) are selected for exhibition of rich evolution law of dynamic behavior. The subjectively selected frequency points are all near the resonance frequency to facilitate the observation of the nonlinear motion window. It can be concluded that the nonlinear behavior is sensitive to excitation frequency and amplitude. With the increase in excitation amplitude, the system tends to be unstable, manifested as jumping discontinuity, quasi-period and chaos. The specific motion state is marked in Fig. 15a1-d1. From the result of 3-D frequency spectrum, the primary frequency f almost dominates the spectrum and its frequency multiplication such as 2f and 3f is also captured. The jumping discontinuous behavior of spectrum is triggered, like non-smooth border of frequency f. Besides, the continuous frequency  Hz and E = 6600, phase trajectories repeat themselves completely every five periods and have no differences between adjacent trajectories. Five isolated points occupy Poincaré section. There are f/5, 2f/5, 3f/5… in spectrum, and non-integral multiples and continuous frequency disappear. It suggests that 5 T-period motion is found. At f = 2370 Hz and E = 12,000, chaos is captured and its phase trajectories almost fill the entire phase space. Poincaré section is consisted of 3000 irregular distributed Poincaré points. The typical continuous frequency occurs. Further, at f = 2370 Hz and E = 12,600, chaos evolves into quasi-period motion. The phase trajectories are sparser than that of chaos. The disordered Poincaré points transit to a closed track. The continuous frequency becomes clearer, but non-integer frequency still exists and leads the spectrum. Similarly, at f = 2500 Hz and E = 27,300, 28,800, the same developed phenomena of vibration response is found when chaos evolves into quasi-period motion.

Effect of worktable position on dynamic response
In this section, the position-dependent property of BSFS is investigated when excitation amplitude E = 5000 is fixed and different worktable position L p is used. Figure 17 shows the bifurcation diagram with excitation frequency as the control parameter when L p = 100, 200…700 mm (the screw length of test bench is 800 mm) and f [ (0, 2200 Hz). The selected external amplitude and frequency are consistent with Fig. 14. In Fig. 17a1-g1, the main resonance frequency point in the feed direction reveals a significant difference. The lowest frequency is 135 Hz corresponding to the middle position of screw shaft. As worktable is close to two supported ends, the main resonance frequency increases symmetrically and they are 150 Hz for 300 and 500 mm, 150 Hz for 200 and 600 mm, 165 Hz for 100 and 700 mm, respectively. This phenomenon is different from the result in Fig. 11. It may be attributed to the effect of load-  Fig. 17a2-g2, the unstable frequency region is slightly widened as worktable is far from fixed end of matched bearings. It also demonstrates the existence of unstable region when f [ (0, 1400 Hz) and L p = 100 and 700 mm, as shown in Fig. 17a1 and g1. The detailed development process of nonlinear behavior is displayed in Fig. 18. Only single-period motion is captured when f [ (0, 1400 Hz) and L p = 200, 300… 600 mm. From Fig. 18a1, the six detailed responses are selected to observe the evolution of orbit and Poincaré section with the increased excitation frequency, and the result is depicted in Fig. 19. It is a complete process for a stable period-3 transiting into chaotic motion. The distinct process is formed through period-3 of f = 224 Hz, period-6 of f = 228 Hz, period-12 of f = 231 Hz, quasi-period of f = 231.4 Hz, chaotic motions of f = 231.8 and 232 Hz. It is a kind of common evolution path that changes a stable periodic motion into chaos. Chaos happened in BSFS means that system works in a state of unpredictability and non-repeatability. For the actual mechanical equipments, chaos is partly detrimental to the working performance and causes noises. According to the bifurcation diagram, the frequency range that probably causes the chaotic motion should be avoided in order to ensure the stability and reliability of system [44].
When platform position L p is chosen as the control parameter, bifurcation and 3-D spectrum are displayed in Fig. 20. The harmonic excitation is fixed at f = 206 Hz which is close to the resonance frequency of feed direction. It can be concluded that BSFS tends to lose the stability and period-n and quasi-period motions appear when platform is located at the both end of screw shaft. At L p [ (200 mm, 600 mm), Fig. 20a indicates that BSFS goes through singleperiod motion and only frequency f and its integer multiples occur in spectrum of Fig. 20b. Additionally, period-n motion is always accompanied by f/n frequency components, such as period-2 motion and f/2, period-3 motion and f/3, 2f/3, as shown in Fig. 20b and a1, a2. The effect of worktable position demonstrates that the dynamic behavior of system behaves differently when platform is located at various positions of screw shaft. During the cutting operation, for structural stability and controllable processing performance, the middle part should be ensured in the actual working stroke of screw to avoid the limit positions at the two ends prone to unstable motion.

Effect of screw shaft length on dynamic response
The screw length is chosen as the control parameter, and various external excitation and fixed platform position of L p = L/2 are employed. The short screw stroke L [ (0.5 m, 1.5 m) is used due to the installation that one end of ball screw is fixed and the other end is supported. In Fig. 21a1-a3, when higher excitation amplitude E is used, the screw length region corresponding to unstable motion is enlarged, and chaotic motion is even captured. From the bifurcation diagrams of Fig. 21b1, a1, b2 and b3, the jumping discontinuous, quasi-period, chaos and single-period are found when f is 1500 Hz, 2000 Hz, 2370 Hz and 2500 Hz respectively. In the corresponding 3-D frequency spectrum, fundamental frequency f dominates the frequency domain, and diverse frequencies are captured which demonstrate the rich nonlinear behavior with the variation of screw length. The parameter range of continuous frequency is consistent with chaotic motion area. The screw length has a substantial impact on the nonlinear behavior of BSFS. However, the regular evolution law is hardly generalized with the variation of screw length since dynamic characteristic is sensitive to external excitation. The four detailed responses from Fig. 21a3 and b2 are selected for the observation of quasi-period and chaotic motions, and the result is depicted in Fig. 22.
Based on the analysis of screw shaft length, it can be found that the stable and controllable motion can be achieved by adjusting the supporting span of screw when system is in an unstable state. In the optimization design stage of dynamic parameters of ball screw feed system, screw length should also be considered as a key factor for a high-performance feed system.
3.5 Effect of preload on dynamic response

The influence of carriage preload
Platform supported by four carriages bears the most of vertical and lateral loads and torque. The carriage preload is the main factor that affects the stiffness and dynamic characteristics of BSFS. Figure 23 exhibits the bifurcation diagram and 3-D spectrum when carriage preload is chosen as the control parameter, f = 2000 Hz, L = 800 mm, L p = L/2 and E = 10 k, 20 k, 30 k, respectively. The frequency near the resonance zone is selected for a clear observation of nonlinear behavior development. Larger excitation amplitude is to obtain richer nonlinear phenomena for the purpose of the theoretical investigation. For Fig. 23a1, a2, system develops into the stable singleperiod from quasi-period when larger carriage preload is employed. The frequency component except f and its multiples gradually disappears with the increased carriage preload, and basic frequency f leads the frequency domain. For Fig. 23b1, b2 and c1, c2, when larger amplitude E is employed, the unstable region is expanded and jumping discontinuous is observed from bifurcation diagram and frequency component f. The conclusion is similar to that in the above b Fig. 17 Bifurcation diagram with excitation frequency as the control parameter under different worktable position when excitation amplitude E = 5000 Fig. 18 Partial bifurcation diagram with excitation frequency as the control parameter at L p = 100 and 700 mm when E = 5000 section. Besides, in the spectrums of Fig. 23b2, c2, the highest frequency is still basic frequency f.

The influence of nut preload
The nut preload is taken as the control parameter, and other condition is the same as Sect. 3.5.1. When E = 10 k in Fig. 24a1, a2, system evolves from singleperiod through quasi-period to single-period motion.
In spectrum, jumping fails to be captured and noninteger multiple frequency occurs. When E = 20 k in Fig. 24b1, b2, system develops into single-period from quasi-period. At preload [ (0, 2016 N), several single-period and period-n motion windows are observed. When E = 30 k in Fig. 24c1, c2, as preload increases, system transforms to stable single-period from period-n and quasi-period motions. Besides, larger excitation amplitude triggers the jumping phenomenon in the bifurcation diagram and spectrum. The characteristic of period-31 motion is illustrated in Fig. 25 when preload is 1050 N and E = 30 k. The phase trajectories repeat themselves completely every 31 periods. Poincaré section is formed by 31 discrete points. The shape is similar to quasi-period motion, but a closed curve loop fails to appear. There are 2f/31, 29f/31… in the frequency domain, and non-integral multiples and continuous frequency disappear.
b Fig. 19 The evolution of orbit and Poincaré section under different excitation frequency when L p = 100 mm and E = 5000

The influence of matched bearing preload
The effect of matched bearing preload on stability of BSFS is investigated through bifurcation diagram and 3-D frequency spectrum, and other conditions are consistent with preceding section. From the result of Fig. 26, the system shows unstable period-n and quasiperiod responses at low values of preload. As preload increases, system transforms to stable single-period motion. Furthermore, at E = 30 k, the discontinuous behavior of bifurcation is captured. The results demonstrate the effectiveness of larger preload in suppressing non-period motions of BSFS, thereby improving its dynamic response. The main purpose of applying a preload is to eliminate the backlash between ball and raceway, and increase the repeatability and rigidity of ball screw feed system. From bifurcation diagram and 3-D frequency spectrum, it can be concluded that adding preload is a feasible method to obtain the good stability of feed system. However, too much preloading results in high friction and heat generation [23]. This also has a negative impact on the quality of the feed motion. Therefore, the preload at a reasonable level is very important.

Effect of damping on dynamic response
In this section, damping coefficient between rolling interfaces is used as the control parameter when f = 2 kHz, E = 5000, L = 800 mm, L p = 400 mm. The preloads of carriage, nut and bearing are 623 N, 1690 N and 290 N which are light and medium preload levels less than 5% of rated dynamic load. The system reveals the stable single-period motion at the large value of damping. When the damping coefficient of carriage and nut is at a low value, the period-n and quasi-period motions of system are found by Fig. 27a1, b1. The spectrum shows a rich component including non-integer and integer multiples, and the amplitude of primary frequency f changes nonlinearly. In Fig. 27c1, c2, at damping of matched bearings [ (0, 320), a non-periodic chaotic motion of system is captured and continuous frequency is found in spectrum. With the increased damping of matched bearings, period-n and quasi-period motions appear alternately at (320, 8240). Then, system enters single-period motion and non-integer multiples disappear in spectrum. The corresponding detailed response windows are exhibited by Fig. 28. The common evolution path from chaos through quasiperiod to period-19 motion is displayed. The stability analysis shows that as the damping coefficient of rolling joint increases, dynamic behavior of system transits to an increasingly stable periodic motion. In other words, the results demonstrate the effectiveness of a higher damping coefficient in suppressing nonperiodic motions of system, thereby improving its dynamic response. By on-purpose designing, the proper damping coefficient could be acquired through adjusting the related structure and the physical parameters of BSFS, such as materials with high damping characteristics. As a result of that, the system is expected to have narrow interval of chaotic motion, extended life time, enhanced reliability and lower noise.

Conclusion
This paper contributes a comprehensive dynamic model of BSFS to investigate the nonlinear behavior induced by nonlinear joint characteristics and coupling effect. The combined loads from six directions are applied on the worktable, and the derived translational and angular deformations of each joint are calculated. The restoring force expression comprising coupling displacement of each joint is obtained. According to mass-spring-damping model, the analytical 18-DOF dynamic equation is established. Model verification tests are carried out, and the results show a good match between simulation and measurement. Effect of external excitation, worktable position, screw length, preload and damping of joint on b Fig. 22 Vibration response under certain excitation and screw shaft length  (1) The jumping phenomenon, super-and subharmonic solutions, soften and harden effects are captured by frequency response curve. The multiple resonance zones of nonlinear system are excited. The different main resonance frequencies among six directions are revealed, and the lowest frequency is in feed direction and the highest one is in pitching direction. Larger excitation amplitudes lead to higher peak-topeak value of response and postpone the resonance frequency corresponding to jumping point. (2) The identical evolution path of nonlinear behavior in six directions is discovered by bifurcation diagram. As the excitation amplitude increases, the unstable region is widened and jumping discontinuous behavior is triggered. The integral and non-integral multiples of primary frequency, as well as continuous frequency, are captured in spectrum corresponding to quasi-period and chaotic motions. (3) The lowest main resonance frequency point in feed direction appears in the middle position of screw. As worktable is close to two supported ends, the main resonance frequency increases symmetrically. The unstable frequency area is slightly expanded as the worktable is far from fixed end of matched bearings. Additionally, BSFS tends to lose the stability when worktable is located at the both end of screw. (4) Larger preload and damping coefficient of rolling joints are beneficial for the stability of system since the unstable quasi-period and chaotic motions develop into single-period motion with the increased preload and damping.
Possible practical applications of the developed model can be concerned with optimizing structure, reducing the unstable motion state and vibration  amplitude, improving accuracy of virtual prototype, etc.
Author contribution Mengtao Xu contributed to methodology; investigation; experimental; writing-original draft; writing-review and editing. Changyou Li contributed to methodology; resources and supervision. Hongzhuang Zhang contributed to supervision; writing-review and editing. Zhendong Liu contributed to supervision; writing-review and editing. Yimin Zhang contributed to resources and supervision.

Declarations
Conflicts of interest There is no conflict of interest in the submission of this manuscript, and the manuscript is approved by all the authors for publication.
Availability of data and material All the experimental data used in this paper have been deposited into the Northeastern University Library and are publicly available for download.
Code availability Not applicable.
Ethical approval Not applicable.
Consent to participate Not applicable. The article involves no studies on humans.
Consent for publication All authors have read and agreed to the published version of the manuscript.