An improved algorithm to predict the pose-dependent cutting stability in robot milling

Chattering is one of the most important factors affecting productivity of robot machining. This paper investigates the pose-dependent cutting stability of a 5-DOF hybrid robot. By merging the complete robot structural dynamics with the cutting force at TCP, an effective approach for stability analysis of the robot milling process is proposed using the full-discretization technique. The proposed method enables the computational efficiency to be significantly improved because the system transition matrix can be simply generated using a sparse matrix multiplication. Both simulation and experimental results on a full size prototype machine show that the stability lobes are highly pose-dependent and primarily dominated by the lower-order structural modes.


Introduction
Thanks to many advantages such as lower cost, acceptable accuracy, flexibility, and reconfigurability, industrial robots have become popular for the in situ large parts manufacturing, light machining, milling, and drilling for example [1]. However, due to the pose-dependent static and dynamic behaviors [2,3], it is an essential requirement to develop an effective algorithm that enables the cutting stability to be rapidly and accurately predicted for the optimization of cutting parameters and the selection of suitable configurations.
Chatter is a type of self-excited vibration growing quickly in metal cutting [4]. Based on the dynamic model of milling process described by delayed differential equations (DDEs), extensive research has been conducted on the chatter stability analysis method including experimental, analytical, and semianalytical ones. Experimental methods [5,6] are the most accurate ones to get the stability boundaries in milling, but their setups are too expensive and time consuming. The most favorite analytical method is the zero-order approximation (ZOA) method in frequency domain presented by Altintas and Budak [7] due to its high efficiency. Though the ZOA method could satisfy most of milling processes, it cannot predict additional stability regions in the case of low radial depth of cut. Multi-frequency method [8] was proposed to overcome this problem with the sacrifice of computational efficiency. Semianalytical methods including temporal finite element analysis (TFEA) [9], semi-discretization method (SDM) [10,11], and full-discretization method (FDM) [12,13] usually solve the DDEs of milling process in time domain. The TFEA method has sufficient accuracy for small times in the cut besides the cases of full immersion. The SDM is another widely used method in milling stability analysis due to its ability for complex working conditions [14,15]. It formulates the response function of the robot milling dynamics as the delayed differential equations (DDEs) using direct integration scheme. By means of linear approximation of the delay items and zeroorder approximation of time-period directional matrix, the DDEs are transformed into ordinary differential equations (ODEs). After the construction of the transition matrix over one time period, the stability can then be determined by the Floquet theory. The FDM simultaneously approximates the state, delay, and time-periodic items in response function with linear interpolation. This treatment dramatically reduces the computational time with nearly the same accuracy compared with the SDM [16]. However, there still exist a huge number of high-dimensional matrix multiplications in the construction of transition matrix by FDM.
Cutting stability analysis in robot machining has been intensively investigated in recent years. Similar to conventional machine tools, two types of chatters can be encountered in practices, i.e., the regenerative chatter and the mode coupling chatter. The former results from the modulation of the chip thickness caused by the repeated cutting on the same surface [4], and the latter occurs if the structural modes are closely matched in the principal axes [17]. Pan et al. [18] seem to be the first to investigate the chatter mechanism in robot milling, and concluded that the cutting instability is mainly caused by the lower-order mode coupling effect. Cordes et al. [19] carried out an experimental investigation on an articulated robot in various cutting conditions. The results show that the regenerative chatter is primarily dominated by the local tool-spindle modes at high milling speed, while the mode coupling chatter is mainly dominated by the lower-order structural modes at low milling speed. Since the structural modes of robots are highly configuration-dependent, investigations of the variations of TCP dynamics of serial robots and its consequences on the milling stability were conducted by Mejri et al. [20] based on experimental modal analysis, and Mousavi et al. [21] based on multi-body dynamic model. Both of them found that the stability lobe in robotic milling is highly configuration-dependent. Mohammadi and Ahmadi [22] used a single degree-of-freedom system with the nonlinear restoring force to model the vibration response of a serial machining robot at its TCP. The proposed model can accurately capture the nonlinear behavior of the system, which can be used to improve the accuracy of predicted stability boundaries. Compared to serial robots, the pose-dependent structural modes of PKMs are more complex. Tunc and Shaw [23] experimentally investigated the configuration-dependent dynamics and stability of a Stewart platform-type robot. Influences of structural mode coupling and feed directions on stability lobes were observed in their study. However, experimental identifications are very expensive and time consuming over the entire robot's workspace. Wilck et al. [24] presented a method for optimizing process stability by taking advantage of the axis positiondependent dynamic behavior of machine tools and the flexible usability of spherical milling tools in a simultaneous 5-axis milling process. Law et al. [25] studied the reduced dynamic modeling method of a parallel robot. Based on the proposed model, the position-dependent limit depth of cut was simulated. However, the stability analysis of the robot milling with consideration of multiple order modes over entire workspace is still time consuming.
Taking a novel 5 DOF hybrid robot named TriMule as an example, an improved algorithm is proposed to realize rapid and accurate prediction of pose-dependent stability over entire workspace. Combined with dynamic model built in Sect. 2, a dynamic milling process model is presented in Sect. 3 to illustrate the pose-dependent characteristic of hybrid robot. In Sect. 4, the improved algorithm considering both efficiency and accuracy is depicted by introducing sparse matrix multiplication in the course of constructing transition matrix. Experimental verifications are conducted in Sect. 5 to show the effectiveness of the proposed method and some useful conclusions are given in the last section.
2 Dynamic modeling Figure 1 shows a 3D view of the TriMule robot under consideration. The robot essentially consists of a 1T2R (T-Translation, R-Rotation) spatial parallel mechanism for positioning and an A/C wrist for orientating. A spindle is attached to the rotary part of the A-axis, on which a tool-tool holder system is installed. In order to formulate dynamic equations of both robot and too-tool holder in the same coordinate frame, a body-fixed frame K C is attached to the tool center point (TCP) C with the w axis coincident with the tool axis, and the u axis parallel to the A-axis as clearly depicted in Fig. 1. Hereafter, all vectors and matrices are assumed to be expressed in K C unless indicated otherwise. For the detailed description of the robot, please refer to [26,27]. In order to investigate the pose-dependent cutting stability of the robot, it is necessary to develop a complete model of the system by considering the following: (1) dynamics of the robot itself, (2) dynamics of the tool-tool holder, and (3) dynamics of the interface between the spindle and tool holder. It should be noted that the spindle is treated as a rigid body rigidly connected with the robot because of its much higher rigidity compared with the tool-tool holder.
First, dynamics of the robot are modeled using the semianalytical approach proposed in Ref. [28] such that the full set of lower-order modes affecting dynamic compliances at TCP can be predicted effectively and accurately using only fourteen generalized coordinates. Consequently, the governing equations of the robot can be formulated by where A ∈ ℝ 6 denotes the nodal displacement vector at A ; ∈ ℝ 8 denotes the collection of the first-order bending modes of three actuated limbs in conjunction with the tor- where B ∈ ℝ 6 and C ∈ ℝ 6 denote the nodal displacement vectors at B and C , respectively, with ∈ ℝ 3 and ∈ ℝ 3 being the linear and angular displacement vectors at C . Similarly, f B and f C denote the nodal force vectors applied at the corresponding nodes with f ∈ ℝ 3 within f C being the cutting force vector applied at TCP. In addition, M B B denotes the spatial mass matrix of the tool holder about B.
Thirdly, the interface between the spindle and tool holder is treated as a virtual joint [30] having three linear/torsional springs along/about the axes parallel to those of K C as shown in Fig. 1. The contact stiffness coefficients of the virtual joint can be easily identified by the IRCSA method [31,32] using the measured FRFs of TCP by means of impact test [33]. These considerations lead to the following compatibility conditions where K c ∈ ℝ 6×6 denotes the stiffness matrix of the virtual joint.
Finally, merging three threads given in Eqs.
(1)-(3) results in the complete dynamic model of the entire system.

Dynamic milling process modeling
Equipped with the equations of motion of the robot developed in Sect. 2 at hand, the robot milling process is formulated by considering the cutting force applied at TCP. The cutting force in milling process is closely related to both chip thickness and tool geometry. As shown in Fig. 2, for a cylindrical milling cutter with diameter D and teeth number N , a rotation frame K i is firstly established fixed to tooth i ( i = 1, … , N ) with its w i axis coincident with the w axis. When the tool rotates at a constant rate ( rev/min ), the dynamic chip thickness of tooth i induced by self-excited vibration can be formulated as. (3) T sional deflections of the A/C wrist; and f A ∈ ℝ 6 denotes the interface forces applied at A. Then, dynamics of the tool-tool holder are established by treating the tool as a 2-node spatial beam element [29] while the tool holder as a rigid body such that Based on the general formulations presented in Ref. [34], the dynamic cutting force applied on tooth i evaluated in K i can be expressed by where R i denotes the orientation matrix of K i with respect to K C ; T(t) denotes the time-periodic directional matrix of the dynamic cutting force. Substituting Eq. (9) into Eq. (4) leads to the dynamic model of the robot milling process It can be seen from Eq. (10) that there exists not only the regenerative effect, but also structural mode coupling and load mode coupling effects in the robot milling process. Thus, the proposed model can be used to predict the stability lobes influenced by various effects.

Stability analysis
In this section, an improved stability analysis method with high efficiency is proposed at first. Then, the pose-dependent stability lobes of the TriMule robot are investigated. Finally, the distributions of limiting stable depth of cut over entire workspace, and with feed directions are predicted to guide optimizations of cutting parameters and configurations.

Improved full-discretization method
The improved stability analysis method proposed here is based on the well-known full-discretization method (FDM) [12]. To investigate the influence of specific modes on stability lobe, the dynamic model of the robot milling process in Eq. (10) is first transformed into modal space. The relationship between modal coordinate p(t) and physic coordinate is Fig. 2 Dynamic model of the robot milling process Tool Workpiece where K t , K r , and K a denote the tangential, radial, and axial cutting force coefficients; a p is the axial depth of cut; and w i (t) is the unit windows function determining whether tooth i is in the immersion zone where en and ex represent the entry and exit angles; en = 0 and ex = arccos(1 − 2a e ∕D) for up milling, en = π − arccos(1 − 2a e ∕D) and ex = π for down milling; and a e is the radial depth of cut.
Transforming Eq. (7) into K C and summing the cutting forces contributed by all teeth, the total cutting force can be formulated as where denotes the mode shape matrix. The column vectors of matrix denote the specific mode shape vectors of the system considered in the stability analysis.
Substituting Eq. (11) into Eq. (10), the dynamic model expressed in modal space is derived as where C represents the proportional damping matrix.
Based on the state space method, Eq. (12) can be transformed into the following form where A denotes the constant coefficient matrix; B(t) denotes the time-periodic coefficient matrix with the principal period at tooth passing interval T , i.e., B(t + T) = B(t).
To solve Eq. (13) numerically, T is divided into m − 1 time intervals (length = T∕(m − 1) ) firstly. Based on the direct integration scheme [12], the response function of Eq. (13) at any time t(t k ≤ t ≤ t k+1 , t k = (k − 1) , k = 1, … , m ) with initial condition y k = y(t k ) can be obtained as According to Eq. (14), y k+1 can be obtained as The easiest way to handle the integral items B(t k+1 − ) , y(t k+1 − ) , and y(t k+1 − − T) in Eq. (15) is to linearly approximate it using two boundary values, thus ; the Moore-Penrose generalized inverse [35] of I − P k+1 can be used to substitute it when it is singular.
Using the discrete map matrix D k defined in Eq. (20), the transition matrix of the system over one periodic time interval can be easily constructed as According to Floquet theory, the stability of the system depends on the eigenvalues of the transition matrix . Only when the moduli of all of eigenvalues are less than 1, the system was stable.
It can be seen from Eq. (21) that the construction of containing a number of high-dimensional matrix multiplication operations. The specific number of multiplications is determined by the discretization of the tooth passing interval T . Since the discrete map matrix D k is a sparse matrix, (17) the computational time can be reduced extensively using the sparse matrices multiplication. According to Eqs. (20) and (21), the transition matrix k+1,1 mapping 1 to k+1 can be expressed as where k,1 denotes the transition matrix mapping 1 to k . As shown in Fig. 3, the non-zero elements in matrix D k can be divided into several partitioned matrices having the following forms S k,1 = Q k+1 (P 0 + P k ), S k,2 = Q k+1 P k+1 Q k+1 P k , S k,3 = I and, matrix k,1 can be divided into T k,1 , T k,2 , and T k,3 accordingly.
For this matrix multiplication, only these partitioned matrices contribute to the result. Thus, the transition matrix k+1,1 can be expressed as In the construction of k+1,1 considering l modes of the system, the multiplication of two 2l(m + 1) × 2l(m + 1) matrices is reduced to two multiplications of lower dimensional matrices by the proposed method. This

Milling stability analysis
In this section, the pose-dependent milling stability of the TriMule robot is investigated. Figure 4 shows definitions of the workspace and some relevant structural and dimensional parameters of the robot. The reference configuration is defined as q 3,1 = q 3,2 = q 3,3 and q 3,4 = 0.5(q min + q max ) . Here, q 3,i ( i = 1 ∼ 4 ) denotes the lengths of actuated and passive limbs; q max and q min are the maximum and minimum lengths of the RP limb. In addition, the tool axis is assumed to be always perpendicular to the x-y plane. As shown in Fig. 4, the middle cross-section of the cylinder is consequently designated as the reference plane. The dimensional and workspace parameters of the robot are given in Table 1.
According to Ref. [19], stability lobes in low spindle speed robot milling are dominated by the lower-order structural modes of the robot. Considering that the lowerorder structural modes of the robot exhibits highly posedependency, the stability lobes at four typical configurations, i.e., configurations P 0 −d H + h 1 , P 1 0 R∕2 − d H + h 1 , P 2 R∕2 −d H + h 1 , and P 3 0 − R∕2 − d H + h 1 are first calculated by the proposed method. The material, cutting force coefficients, and cutting parameters used in the simulation are listed in Table 2. The milling tool is perpendicular to the x-y plane and feeding in x direction at each configuration. Figure 5 shows the predicted stability charts. As is evident, not only the distribution of the stability boundaries, but also the limiting stable depth of cut exhibits strong pose-dependency in Titanium alloy milling process. This phenomenon is mainly resulted from the pose-dependency of lower-order structural modes of the robot, which confirms the necessity for considering the influence of cutting configurations in low spindle speed robot milling.
In order to guide the selection of cutting configuration of the robot under specific cutting condition, Fig. 6 shows the distributions of the minimum stable axial depth of cut over the entire workspace. The milling tool is perpendicular to the x-y plane and feeding in x direction at each configuration. It can be seen that the planar symmetrical structure of the robot makes the distribution of the minimum stable axial depth of cut also planar symmetry. The maximal value of them appears at the center point on the top plane of the workspace. It can also be found that these values decrease monotonically from center axis to the edge of the workspace, indicating that the TriMule robot has better cutting performance nearing the center axis of the workspace. The major features of this distribution are almost in agreement with the distributions of lower-order natural frequencies of the robot according to Ref. [28]. The minimum stable axial depth of cut is inversely proportional to the dynamic compliance of lower-order modes of the robot, indicating that the smaller dynamic compliance is, the larger stable cutting region is.
For further investigations of the influence of feed directions on milling stability, Fig. 7 gives the distribution of feed-direction-dependent minimum stable axial depth of cut at the configurations along the center axis of the workspace (see Fig. 4). The region inside the envelope is stable, while outside is unstable. It can be seen from Fig. 7 that the larger the elongation of robot along z-axis is, the smaller the minimum stable axial depth of cut is. The largest and smallest minimum stable axial depths of cut appear when feeding in x-direction and y-direction at each specific configuration, respectively. This phenomenon is probably resulted from the difference of the 1st-order structural mode in different major forced directions.
For a better illustration of the computational efficiency of the improved FDM, a comparison study is made on a workstation (Intel i5-3320 M CPU and 8 GB RAM) as follows. Compared with 145 s taken by the traditional FDM at a specific configuration, only about 68 s was needed by the proposed improved FDM. Computational efficiency is increased by 53%, leading to drastic time-savings for robot milling stability analysis over entire workspace. These observations fully consolidate the usefulness of the proposed model for predicting pose-dependent stability of the robot.

Experimental validation
This section takes a series of milling experiments with the TriMule robot (see Fig. 1) to verify the efficiency of the proposed method. The contact parameters between tool holder and spindle in Eq. (3) are firstly identified through impact test at TCP. The detailed parameters are shown in Table 3. As a check to verify the effectiveness  Fig. 7 The distribution of the feed-direction-dependent minimum stable axial depth of cut of the proposed dynamic model, a comparison study between predicted and measured TCP FRFs was made as shown in Fig. 8. It can be seen that the major features of predicted FRFs are in good agreement with the measured FRFs, and the first eighth-order modes of the robot are focused on 0-300-Hz spectrum. The 9th-12thorder modes are the local modes of spindle and tool-tool holder.
In low speed milling tests, the material of the workpiece is the same as the simulations in Sect. 4.2, with the relevant K c = diag k u , k v , k w , k ru , k rv , k rw parameters as shown in Table 2. The tool axis is aligned perpendicular to the worktable plane, and the feed direction is chosen in x-direction in all cutting tests. Figure 9a shows the experimental setup, where a LMS dynamic analyzer is used to process vibration signals measured by accelerometer (PCB 356A26). Classical "step-cutting" strategy is applied to check the stability with the radial depth of cut keeping constant. Unstable cutting is assumed when additional dominant frequencies of vibration signals arise in the frequency spectrum.
To investigate the influence of cutting configurations on the stability, four typical configurations (as shown in Fig. 9) the same as simulations were chosen. At each single configuration, 12 cutting tests with spindle speed as 800, 1200, and 1600 rev/min, and axial depth of cut as 1, 2, 3, and 4 mm were conducted. The radial immersion and feed rate were assigned as 10% and 0.05 mm/tooth in all tests, respectively. Figure 10 shows the comparisons of measured and predicted stability results in low spindle speed milling at four different configurations, from which it can be found that the predicted stability lobes match well with those obtained by the experiments. As is evident, though the cutting parameters are identical, the limiting axial depths of cut are different at different configurations. This phenomenon confirms that the posedependent lower-order structural modes of the robot lead to the pose-dependent stability.
For a further investigation of chatter frequencies in different cutting configurations, Fig. 11 gives the finished surface characteristics and vibration signals in frequency domain of four cutting conditions (marked as A-D in Fig. 10). Obviously, conditions A-D with the same cutting parameters but different cutting configurations exhibit different stability. Undesirable finished surface is observed when chatter occurs in conditions C and D. It can be found from vibration signals that the chatter frequency are 31.09 Hz and 30.38 Hz in conditions C and D respectively, which are a little higher than the natural frequency of the 2nd-order structural mode of the robot. It may result from that the cutting stability in this condition is dominated by the lower-order mode having maximum dynamic compliance in feed direction. High-frequency spindle and tool-tool holder's local modes are damped out by cutting process at low speed range cutting (400-2000 rev/min). Only tool passing frequency and its harmonics are observed in conditions A and B without chatter.
Thus, this paper provides an efficient tool for the prediction of pose-dependent chatter stability in robot milling over entire workspace, and lays a solid ground for subsequently optimizations of cutting parameters and cutting configurations.

Conclusion
This paper deals with the pose-dependent milling stability analysis of a novel 5-DOF hybrid robot named TriMule. The conclusions are drawn as follows.
1. By combining RCSA method with experimental modal analysis, a dynamic model of the TriMule robot which can predict its full frequency spectrum dynamic behaviors within the entire workspace is first established. Then, a robot dynamic milling process model considering regenerative effect, structural, and load mode coupling effects is derived. By introducing sparse matrix multiplication operations in the construction of the system transition matrix, an improved full-discretization method for stability analysis is proposed. These treatments lead to dramatic time-savings in the stability anal-ysis over entire workspace. The proposed method can be used to predict the limiting stability boundaries of the robot over entire workspace with sufficient efficiency and accuracy. 2. The milling stability of the TriMule robot has been analyzed using the proposed method. Simulation results show that the stability lobes exhibit highly pose-dependency, which is mainly resulted from the pose-dependency of the lower-order structural modes of the robot. The limiting stable depth of cut decreases monotonically with the elongation of the robot along z-axis, and is largest in the center of top plane of the workspace. At a specific configuration, the largest limiting stable depth of cut appears when feeding in x-direction. The effectiveness of the proposed method has been verified by a comparison study against experimental cutting tests. The chatter frequencies appeared close to the 2nd-order structural