Abstract

In the machining process that material removal involves, the time-delay effect due to the regenerative cutting is the root cause of tool chatter, which is a severe nonlinear vibration that leads to system failures. The Chebyshev collocation method (CCM) can be applied to the stability analysis of the delay-affected machining systems. However, when the degree-of-freedom (DOF) of the system is high, the computational efficiency of the Chebyshev collocation method is far lower than the commonly used semidiscretization and full-discretization methods (FDMs). In this article, a robotic milling model allowing an arbitrarily high degree-of-freedom is proposed as a benchmark to evaluate the computation performance for different stability analysis algorithms. An improved algorithm named the “fast Chebyshev collocation method” (FCCM) is proposed to handle the delay differential equation (DDE) with a high degree-of-freedom. The proposed fast Chebyshev collocation method accelerates the traditional Chebyshev collocation method in two approaches: one is the inversion of the matrix when constructing the transition matrix, and the other is the reduction of the dimension of the transition matrix by applying the Sherman–Morrision–Woodbury formula. Subsequently, both the full-discretization method and the proposed method are applied to the robotic milling system to show their convergence rate, computation efficiency, and accuracy. The results demonstrate that the proposed method is overall advantageous to the full-discretization methods in convergence, efficiency, and accuracy even when the degree-of-freedom is sufficiently high, implying that the proposed fast Chebyshev collocation method can be a potential alternative tool to deal with the stability analysis for time-delay systems.

1 Introduction

Regenerative chatter [1,2] is known as the main factor causing poor surface quality and tool wear in the field of high-speed machining [35]. To avoid chatter, stability analysis of the dynamic system describing the machining system is important. Typically, the dynamic of the machining structure is modeled by modal parameters obtained via experimental modal analysis [6] or finite element analysis (FEA) [7,8]. As for the cutting force model, the linear force model proportional to the chip thickness is commonly used. By considering the wave left on the machined surface by the vibration of the tool in the previous time, the time delay is introduced to the cutting force model. Combining the structure dynamic model and the cutting force model, the governing equation of the machining process is formulated as the delay differential equation (DDE) with constant or time-varying periodic coefficients.

As for the stability analysis of the machining system based on the DDE, since a closed-form analytical solution of DDE is usually impossible to find, many efforts have been taken to find the numerical method to predict the stability of DDE with both high computation efficiency and accuracy. These methods explore the stability of DDE around the equilibrium or periodic orbit under different combinations of parameters, resulting in the stability lobe diagram (SLD). SLD divides the parameter space, usually, the two-dimension parameter space including the spindle speed and the axial depth of the cutting tool, into the stable (chatter-free) and unstable (chatter) zones, then one can easily choose the appropriate parameter to avoid chatter during the practical machining process. Nowadays, many methods have been proposed to predict the SLD of a variety of DDEs. For DDE with constant coefficients, the D-subdivision method [9] can be used by finding the set of parameters under which the characteristic equation has imaginary roots. Sadath and Vyasarayani [10] proposed Galerkin approximations for stability analysis of DDE with periodic coefficients by transforming the original DDE into an equivalent partial differential equation with a time-periodic boundary condition. By implementing the frequency-domain transformation of the DDE, Altintas and Budak [11] first proposed the zero-order approximation method, or single frequency method, to predict the SLD of the milling process, where the critical axial depth has an analytical solution, which greatly enhances the computation efficiency. Then, the multifrequency method [12] was proposed to obtain more accurate results in low radial immersion milling.

Different from the analytical, frequency-domain method mentioned above, Insperger and Stépán [13] proposed the commonly used semidiscretization method (SDM) as a time-domain method to predict the stability of DDE, where only the delay term of the dynamic equation is approximated by zero-order [14] or first-order interpolation polynomials [15], discretizing the original DDE into piecewise ordinary differential equations (ODE). Then, the stability is predicted by calculating the critical eigenvalue of the transition matrix, which is an approximation of the infinite-dimension monodromy operator of the original DDE. Subsequently, Ding et al. [16] proposed the full-discretization method (FDM) as a modified version of SDM, where the coefficients of the state term are decomposed into the time-invariant part and the periodic part, and the state term and the delay term are discretized simultaneously. This modification significantly reduces the calculation times of the matrix exponent, enhancing the computation efficiency. To improve the prediction accuracy of the critical eigenvalue, several high-order SDMs [17] and FDMs [18] were proposed using high-order interpolation polynomial to discretize the state and the delay terms. On the other hand, to improve the computation efficiency, Henninger and Eberhard [19] summarized the numeric techniques to reduce the computation time of SDM, including reducing the dimension of the transition matrix and utilizing the sparsity of the corresponding matrix when implementing matrix multiplying.

In some situations such as interrupted milling with low radial immersion, the tool leaves the workpiece periodically, causing the delay term to vanish during some part of the period. In these time intervals, the original DDE becomes ODE. Based on this observation, Mann et al. [20] proposed the time finite element analysis (TFEA) method to predict the SLD and the surface location error simultaneously. TFEA method separates the DDE into the cutting state and the free state, at the free state, the DDE is replaced by an ODE, improving the prediction accuracy. However, the computation efficiency is low due to the symbolic operation used in the construction of the transition matrix. To overcome this difficulty, Khasawneh and Mann [21] proposed the spectral element method, which can achieve the spectral convergence rate compared to the linear convergence rate in SDM and FDM. Butcher et al. [22] proposed the Chebyshev collocation method (CCM) which is also able to achieve spectral convergence rate. Breda et al. [23,24] proposed a pseudospectral method for the stability analysis of autonomous and periodic linear delayed dynamical systems, which has a fast convergence rate. Given that one of the main drawbacks of the above spectral element method and CCM is that the spectral convergence rate is lost when the coefficients of DDE have discontinuities within the whole period, which limits their use for the machining process when multiple teeth of the tool engaged in cutting simultaneously, Khasawneh et al. [25] proposed the multi-interval Chebyshev collocation method to handle DDE with discontinuities. A comprehensive comparison of the above-mentioned semi-analytical methods can be found in Refs. [26] and [27].

It should be noted that in most literature, only the most flexible mode is chosen to model the structure dynamic of the machining system. But, Wan et al. [28] pointed out that this choice is only valid when the multiple modes are separated far away enough in terms of the modal stiffness, and proposed the lowest envelope method to construct the SLD when two or more modes with the almost modal stiffness co-exist. Moreover, in some cases such as the rotary drilling system with a long, flexible drill pipe, it is convenient to apply finite element model to discretize the original continuous system [7,8]. All these modeling methods result in a DDE with high degrees-of-freedom (DOF). There is little literature discussing and comparing the performance of different semi-analytical methods when applied to DDE with high DOFs, and among these methods, CCM seems to be a potential choice to predict the SLD due to its spectral convergence rate and relatively simple construction of the transition matrix. Totis et al. [29] used both CCM and SDM to predict the SLD in milling with spindle speed variation and compared the performance of those two methods in detail. It was concluded that CCM is superior to SDM both in computation accuracy and efficiency. However, they pointed out in their recent work [30] that the computation performance is not satisfactory especially when the dimension of the transition matrix is large.

Inspired by the work in Ref. [19], this article proposes a fast Chebyshev collocation method (FCCM) which greatly reduces the dimension of the transition matrix. Since it has been shown that FDM computes faster than SDM when the dimensions of the transition matrix are the same [16], this article uses FDM as the reference and compares the computation performance of the proposed method and the FDM. The remainder of this article is organized as follows. In Sec. 2, the dynamic model of a robotic milling system is built; this model allows the use of a finite element model of the cutting tool, providing arbitrarily high DOFs by increasing the number of elements. Combining this dynamic model with the classic orthogonal cutting force model, the DDE of the full system is derived which will be used in Secs. 3 and 4 to test the computation performance of stability prediction algorithms. In Sec. 3, the original algorithm of CCM is briefly introduced, and the computation efficiency of CCM and FDM is compared. In Sec. 4, an acceleration algorithm named fast Chebyshev collocation method is proposed in detail, including the reduction of the dimension of the transition matrix and the acceleration of the matrix inversion. The convergence rate and the computation efficiency of the proposed method are then analyzed comparing to the FDM. Finally, the main conclusions are drawn in Sec. 5.

2 Modeling of the Robotic Milling System

To evaluate the performance of the SLD prediction method under dynamic systems with high DOFs, a robotic milling model allowing arbitrarily high DOFs is proposed in this section. The modeling process of the structural dynamic of the system and the cutting force is introduced in Secs. 2.1 and 2.2, respectively.

2.1 Structural Dynamics of the Robotic Milling System.

A robot milling system consists of three parts, the robotic arm, the spindle, and the cutting tool, as shown in Fig. 1. Since the spindle is fixed on the last joint of the robotic arm, it can be considered as a part of the robotic arm. The modeling contains two parts: the robot-spindle assembly's dynamic model and the cutting tool's finite element model.

Fig. 1
The robotic milling system
Fig. 1
The robotic milling system
Close modal
As for the robot-spindle assembly, by applying the Lagrange equation of the second kind [31], the governing equation of the dynamic system is represented as the following nonlinear form:
(1)
where q,q˙,q¨ are vectors of the size 6×1, denoting the rotation angle, angular velocity, and acceleration of the joints of the robotic arm, respectively. M(q) is the inertia matrix, C(q,q˙) is the torque vector of the centrifuge and Coriolis force matrix, and G(q) is the torque vector of the gravitational force. τd represents the driving torque. For simplicity, we assume that during the milling process, the tool is always perpendicular to the surface of the workpiece and the feed trajectory is a straight line. Furthermore, we assume that the feed speed is constant. By applying the resolved motion acceleration control (RMAC) [32], the governing equation can be expressed in the Cartesian coordinate frame fixed to the end of the spindle, i.e., the end coordinate frame (X1Y1Z1 in Fig. 1) as
(2)

where J(q) is the Jacobi matrix of size 6×6 mapping the angular velocity of the joints into the velocity and angular velocity in the end coordinate frame [31]. ep,eo denote the error between the actual and expected displacement and rotation angle in the end coordinate frame, respectively: ep=(xxd,yyd,zzd)T, eo=(ααd,ββd,γγd)T, where x, y, z denote the displacement and α,β,γ denote the rotation angle. The subscript “d” denotes “expected,” i.e., displacement and rotation angle along the ideal feed trajectory. Scalars kd,kp denote the control parameters used in the RMAC method. The deriving process of Eq. (2) in detail is given in Appendix  A.

As for the flexible cutting tool, it is regarded as a uniform Timosheko beam, and a finite element model is used to discretize it. After assembling the mass and stiffness matrices of all the elements (already have been transformed from the local coordinate frame of the beam element to the end coordinate), the structural dynamic equation of the cutting tool can be written as
(3)

where Mtool, Ctool, and Ktool denote the mass, damping, and stiffness matrix of the whole tool. The vector F is the external force related to the cutting process, refer to Sec. 2.2. Hence, time-delay effects due to the regenerative cutting are involved in this system. Assuming that the structural damping of the tool is Rayleigh damping, the damping matrix is directly obtained by the linear combination of Mtool and Ktool: Ctool=αMMtool+αKKtool, where αM,αK are the corresponding mass and stiffness coefficients. X denotes the vector consisting of all the nodal displacements and rotation angles. Assuming that the control (the RMAC method) applied to the joints has already been stable, and the rotation angle of the joints and the deformation of the tool are all small. Then, X can approximately have the following expression: X=(ep,1Teo,1Tep,n+1Teo,n+1T)T, where ep,i,eo,i denote the error of the displacement and rotation angle of the ith node, respectively. n is the number of elements. Specifically, the subscript “1” denotes the node at the end of the tool (X1Y1Z1 in Fig. 1, where the spindle and the tool are connected), subscript “n+1” denotes the node at the tool tip (Xn+1Yn+1Zn+1 in Fig. 1).

Finally, the dynamic equation of the robot-spindle assembly and the structural dynamic equation of the cutting tool are assembled together, resulting in the governing equation of the whole robotic milling system
(4)

Here, we assume that during the milling process, changes in the elements of the matrix M(q) and J(q) due to the changes in q are small and can be ignored; therefore, the matrices M̂,Ĉ, and K̂ remain to be constant. Here, assuming that the two-dimensional cutting model is applied, then F is formulated as: F=(01×6nFxFy01×4)T, where Fx,Fy denote the milling force applied on the tool tip, i.e., the node n+1, in x and y direction, whose expression in detail is derived in Sec. 2.2.

2.2 Cutting Force Model.

A classic two-dimensional cutting force model applied on the tool tip [33] is shown in Fig. 2. The coordinate frame Xn+1Yn+1Zn+1, i.e., the tool coordinate frame, is fixed on the tool tip, i.e., the node n+1, where Xn+1,Yn+1,Zn+1 denote the axis along the feed direction, perpendicular to the feed direction, and along the tool axis, respectively. a denotes the axial depth of cut, ae denotes the radial depth of cut, and v denotes the feed velocity of the workpiece relative to the cutting tool. Dtool is the diameter of the tool. Ω denotes the spindle speed (unit: rotations per minute (rpm)). ϕj represents the angular position of the tooth j. For simplicity, assuming that the cutting tool has Ntool equally distributed cutting teeth, resulting in a constant pitch angle, then φj(t)=(2πΩ/60)t+(j1)(2π/Ntool).

Fig. 2
The two-dimensional cutting force model: (a) view from the Zn+1 axis and (b) the three-dimensional view
Fig. 2
The two-dimensional cutting force model: (a) view from the Zn+1 axis and (b) the three-dimensional view
Close modal

Ft,j,Fr,j represent the tangential and normal cutting force applied on the jth tooth, and hj denotes the chip thickness of the jth tooth; by using the linear cutting force model and the classic circular tool path model [34], the cutting force components are formulated as: Ft,j=aKthj, Fr,j=aKrhj, and hj=(x(t)x(tT))sinφj+(y(t)y(tT))cosφj, where Kt and Kr are tangential and normal cutting force coefficients, respectively. T denotes the tooth passing period: T=60/(NtoolΩ), which is also the delay of the dynamic system.

According to the previous assumptions, we have xd(t)=x0+vt,yd(t)=y0, where x0 and y0 are the initial position of the tool tip in Xn+1 and Yn+1 direction. Then, x(t)x(tT)=ep,x(t)ep,x(tT)+f and y(t)y(tT)=ep,y(t)ep,y(tT), where ep,x and ep,y are the terms xxd and yyd in ep, respectively. f denotes the feed per tooth, i.e., f=vT=60v/(NtoolΩ). Adding up all the cutting force applied on the tooth engaged in cutting, transforming all the cutting force to the center of the tool tip section (ignoring the resulting equivalent torque) and projecting them to the Xn+1 and Yn+1 direction, the two-dimensional cutting force is represented as
(5)
where H(t) is a time-varying periodic matrix of size 2×2, i.e., H(t)=H(t+T). Window function g is a switch function determining whether the tooth is in cutting
(6)
φst and φex are the entrance and exit angles of the cutting tooth, respectively,
(7)
Then, the final DDE of the robotic milling system is formulated as
(8)

In Eq. (8), M̂,Ĉ, and K̂ are constant matrices of size 6(n+1)×6(n+1), Ĥ1(t) and Ĥ2(t) are periodic matrices of size 6(n+1)×6(n+1): Ĥ1(t)=Ĥ1(t+T),Ĥ2(t)=Ĥ2(t+T), and f̂(t) is periodic vector of length 6(n+1): f̂(t)=f̂(t+T). τ is the delay of the system which is equal to the period T.

As the number of elements n increases, the DOF of Eq. (8) can be arbitrarily high, allowing the evaluation of the performance of the SLD prediction algorithm applied to DDE with high DOFs. All the parameters used to derive Eq. (8), including the parameters of the robotic arm, the flexible cutting tool, and the cutting force are all listed in Appendix  B.

3 Chebyshev Collocation Method and the Evaluation of Computational Efficiency

In this section, the stability analysis of Eq. (8) is performed. Let Z(t)=(XT(t)(M̂X˙(t)+12ĈX(t))T)T, Eq. (8) is rewritten as the first-order state-space form
(9)
where A0 is a constant matrix, and A(t)=A(t+T),B(t)=B(t+T), and F0(t)=F0(t+T). Therefore, Eq. (9) has a periodic solution Z¯(t), i.e., Z¯(t)=Z¯(t+T). Let Z(t)=Z¯(t)+Ẑ(t), where Ẑ(t) represents the perturbation around the periodic solution. Because the delay τ is exactly equal to the period of the coefficient matrix T and A(t) and B(t) have opposite signs, the periodic solution Z¯(t) satisfies
(10)
which is an ODE. Combine Eqs. (9) and (10), we have
(11)

which is a DDE only corresponding to the perturbation Ẑ(t).

Equation (11) is then used to predict the stability of the original DDE. In Sec. 3.1, the stability analysis process using the Chebyshev collocation method is introduced in detail.

3.1 Chebyshev Collocation Method.

The basic idea of traditional CCM is to use an interpolation polynomial to approximate the solution Ẑ(t) in the time interval [0, T]. Due to the cutting tooth entering and leaving the cutting zone, the coefficients corresponding to the cutting force of Eq. (11), i.e., the matrices A(t) and B(t) often have nonsmooth, even discontinuous points. So first the overall time interval [0, T] should be shifted to the new one [tst,tst+T], where t0 is the time when the angular position of one certain cutting tooth is equal to the entrance angle: φj(tst)=φst. This shifting ensures that there is at most one nonsmooth point inside the interval [tst,tst+T]. tst has infinite possible values, and the one with the smallest absolute value is chosen
(12)

where jst is the index of the cutting tooth such that the absolute value |φstφj(0)| reaches the minimum. Then, let ρcut=(φexφst)/(2π/Ntool), which represents the ratio of the arc length of the cutting zone to the pitch angle. When ρcut=1, the coefficients are sufficiently smooth over (tst,tst+T], which means that the original CCM is valid, so only cases when ρcut<1 need to be considered. (For convenience, the situation when ρcut>1 is not considered in this article, but according to Ref. [25], generalization to that case is straightforward.)

When ρcut<1, the milling process is divided into two states, namely, the cutting state when there is only one tooth in contact with the workpiece, and the free state when no tooth is engaged in cutting and the tool is under free vibration. Then, CCM is applied only within the cutting state, since in the free state, the coefficients of A(t) and B(t) are all zeros and the DDE degenerates into an ODE. Therefore, the time interval where Ẑ(t) is approximated changes from (tst,tst+T] to (tst,tst+ρcutT], using the N+1 order Chebyshev points of the second kind on [a, b] as the interpolation points: xi=((b+a)/2)+((ba)/2)cos(iπ/N), i=0,1,,N, the solution Ẑ(t) is approximated as
(13)
ϕi(t) is the ith Lagrange basis polynomial corresponding to point ti
(14)
Then, substitute Eq. (13) into Eq. (11), we get a series of equations, which is the discrete form of the approximation of Eq. (11)
(15)
then rewrite them as the matrix form
(16)
where
(17)
here, the operation is the Kronecker product. Il is the identity matrix of size l×l, and l is the length of vector Ẑ(t). In the context of this article, l=12(n+1). And according to Ref. [6], the matrix D has the following explicit expression, also called the “Chebyshev differential matrix”:
(18)
In addition, due to the existence of the free state, the relationship of the endpoint, i.e., Ẑ(tN) and Ẑ(t0τ), becomes Ẑ(tN)=eA0(1ρcut)TẐ(t0τ) [22]. To enforce this condition, the last l rows of the matrices D,MA, and MB need to be modified, resulting in the following equations:
(19)
where
(20)
and
(21)

where D¯UL denotes the block of (2/(ρcutT))D¯ including the first N rows and columns, while D¯UR denotes the block of (2/(ρcutT))D¯ including the first N rows and the last column. Notice that D̂ is a matrix that only depends on the geometry of the cutting tool, and the period T, which only depends on the spindle speed Ω.

From Eq. (19), the discrete map mapping state terms within [T,0] to [0, T] is immediately obtained
(22)
where
(23)
Φ is the final transition matrix used to predict the stability, which is a finite dimension approximation of the infinite dimension monodromy operator of the DDE Eq. (11). Based on the Floquet theory, the stability of the DDE can be determined according to the spectral radius of the transition matrix ρ(Φ)
(24)

In practical application, the inversion of the matrix (D̂M̂A)1 can be avoided by solving the generalized eigenvalue problem, i.e., finding λ such that the dominant |M̂Bλ(D̂M̂A)|=0.

3.2 Comparison Between Chebyshev Collocation Method and Full-Discretization Method.

Now that the transition matrix Φ can be obtained, the SLD can be constructed by calculating the spectral of Φ under different combinations of the parameters, which are often, in milling, the spindle speed Ω and the depth of cut a. In this article, we apply the method commonly used in SDM and FDM to approximately construct the SLD. Given the range of Ω and a that one is interested in, i.e., [Ωmin,Ωmax] and [amin,amax], then the parameter space (Ω,a) is discretized into NΩ×Na mesh points (Ωj,ai) with constant increment dΩ and da, respectively,
(25)

Na and NΩ represent the number of discrete points along the axis of a and Ω, respectively.

The algorithm of construction of SLD using CCM is shown in Algorithm 1. After performing the algorithm, the spectral radius of the transition matrix ρΦ(Ω,a) is obtained on all the mesh points. Then, the stability border where ρΦ(Ω,a)=1 can be found, dividing the parameter space into stable and unstable zones.

Table

Algorithm 1. CCM

 forj=1 to NΩdo
  Ω=Ωj
  Calculate D̂
  fori=1 to Nado
   a=ai
   Calculate M̂A, M̂B
   Solve the generalized eigenvalue problem |M̂Bλ(D̂M̂A)|=0
  end for
 end for
 forj=1 to NΩdo
  Ω=Ωj
  Calculate D̂
  fori=1 to Nado
   a=ai
   Calculate M̂A, M̂B
   Solve the generalized eigenvalue problem |M̂Bλ(D̂M̂A)|=0
  end for
 end for

To evaluate the computation efficiency of CCM, the commonly used FDM is chosen as the reference method in this article. Define the time spent on constructing the transition matrix Φ and calculate the spectral radius of Φ in single point (Ωj,ai) when using CCM and FDM as teig,CCM and teig,FDM, respectively. Assuming that teig,CCM and teig,FDM remain the same regardless of the Ωj and ai, we do not use the time spent on constructing the whole SLD to assess the computation efficiency of CCM and FDM but teig,CCM and teig,FDM instead.

The DDE model built in Sec. 2 is adopted to apply the algorithm. The other necessary parameters are Ω=500rpm,a=8mm,ae/Dtool=0.2,andNtool=4. All codes of the algorithm are written with matlab 2023a and implemented on the same computer (12th Gen Intel(R) Core(TM) i5-12400 CPU @2.50 GHz 16.0 GB RAM, Windows 11 OS). When applying FDM, the method to improve the computation efficiency proposed in Ref. [12] is adopted including reducing the dimension of the transition matrix and utilizing the sparse matrix multiplying. For reliability of the results, we run the algorithm proposed here and in the remains of this article ten times and use the average value of the running time as the final result.

The result of teig,CCM and teig,FDM in different system DOF l and the number of discrete segments N are shown in Fig. 3. The ratio of the computation time of the single point ρt=teig,FDM/teig,CCM and distribution of whether CCM is faster in the parameter space (l,N) is also shown in Fig. 4. The computation time for both FDM and CCM increases with the DOF l or number of segments N. In most cases except for (l,N)=(24,20), the CCM computes slower than FDM. The ratio of the computation time also increases with the increase of l or N, when (l,N)=(132,60), the ratio ρt is 54.639, meaning that CCM is almost 50 times slower than FDM. The main reason for this inefficiency lies in the enormous difference in the dimension of the transition matrix when applying two kinds of methods. In FDM, according to Ref. [19], the final dimension of the transition matrix Φ is l+2N, while in CCM it is l(N+1). When the DOF of the system l is high, i.e., l2, then l(N+1)l+2N, resulting in a huge computation cost when solving the (generalized) eigenvalue problem compared to the FDM.

Fig. 3
teig,CCM and teig,FDM: (a) value in parameter space (l,N) and (b) N=60
Fig. 3
teig,CCM and teig,FDM: (a) value in parameter space (l,N) and (b) N=60
Close modal
Fig. 4
The comparison of the computation time of the single point (Ω,a) in the parameter space (l,N): (a) the ratio of the computation time and (b) the distribution in the (l,N) space whether CCM is faster or not
Fig. 4
The comparison of the computation time of the single point (Ω,a) in the parameter space (l,N): (a) the ratio of the computation time and (b) the distribution in the (l,N) space whether CCM is faster or not
Close modal

Hence, to improve the computation efficiency of CCM, the first attempt is to reduce the dimension of the transition matrix Φ, which is one of the main topics of Sec. 4.

4 The Accelerated Algorithm: Fast Chebyshev Collocation Method

Based on the traditional CCM, an accelerated algorithm named FCCM is proposed in this section. The acceleration approach includes the reduction of the dimension of the transition matrix and the acceleration of the matrix inversion. The convergence rate and the computation efficiency of the proposed FCCM are then analyzed and compared to the FDM, referring to Secs. 4.14.4.

4.1 Dimension Reduction of the Transition Matrix.

Inspired by the work in Ref. [19] commonly applied in SDM and FDM, the elements of the coefficient matrix B(t) in Eq. (11) are almost zeros except for the 2×2 block H(t); therefore, the term B(t)Ẑ(tτ) depends only on two elements of the delay term Ẑ(tτ). Then, Eq. (11) can be rewritten as the following form:
(26)
where
(27)
The matrix P can be regarded as a selection matrix choosing the elements of Ẑ(tτ) that really contribute to the value of B(t)Ẑ(tτ). Noted that B(t)=B̂(t)P, Eq. (11) can be reduced to the following form:
(28)
where M̂Br is a l(N+1)×(l+2(N1)) size matrix, while MP is a (l+2(N1))×l(N+1) size matrix
(29)
Then, the reduced discrete map can be obtained
(30)
where the reduced transition matrix Φr can be formulated as
(31)

Now Φr is an (l+2(N1))×(l+2(N1)) size matrix that has almost the same dimension as the transition matrix Φ constructed in SDM or FDM when the number of the discrete segments N are equal. The new algorithm constructing the SLD is shown as below.

Table

Algorithm 2 CCM (using the reduced transition matrix Φr)

 forj=1 to NΩdo
  Ω=Ωj
  Calculate D̂
  fori=1 to Nado
   a=ai
   Calculate Φr=MP(D̂M̂A)1M̂Br
   Calculate the spectral radius ρ(Φr)
  end for
 end for
 forj=1 to NΩdo
  Ω=Ωj
  Calculate D̂
  fori=1 to Nado
   a=ai
   Calculate Φr=MP(D̂M̂A)1M̂Br
   Calculate the spectral radius ρ(Φr)
  end for
 end for

Define the time spend on constructing Φr and calculating ρ(Φr) for a single parameter pair (Ω,a) is teig,CCM,r, similar to Sec. 3.2, then the results of teig,CCM and teig,CCM,r are shown in Fig. 5. The spindle speed Ω, the depth of cut a, and other parameters used here and in the remainder of this section are all the same as in Sec. 3.2 unless otherwise noted. Although the dimension of the transition matrix has been reduced significantly, the time spent on calculating the spectral radius of Φr is longer than that when using the original Φ to calculate. When using Eq. (23) to calculate the spectral radius, the inversion of the l(N+1)×l(N+1) matrix (D̂M̂A)1 is avoided by solving the generalized eigenvalue problem instead. However, using Eq. (31) reintroduces the inversion. The time spent on inversion and solving the generalized eigenvalue problem has the same order corresponding to the dimension l(N+1), which prevents the improvement of the computation efficiency. Therefore, effort must be taken to further accelerate the matrix inversion, which is discussed in Sec. 4.2.

Fig. 5
teig,CCM and teig,CCM,r: (a) value in parameter space (l,N) and (b) N=60
Fig. 5
teig,CCM and teig,CCM,r: (a) value in parameter space (l,N) and (b) N=60
Close modal

4.2 Acceleration of Matrix Inversion.

To simplify the matrix inversion, considering the matrix D̂M̂A, it has the following partition:
(32)
then
(33)
Appropriately blocking the matrix M̂Br and MP
(34)
Then, Φr can be formulated as
(35)
Now we need to simplify the calculation of D̂UL,MA1. The matrix D̂UL,MA can be separated into two parts
(36)
D̂ULC is of size lN×lN and only depends on the spindle speed Ω. As for At, similarly, noted that the elements of the coefficient matrix A(t) in Eq. (11) are almost zeros except for the 2×2 block H(t), i.e.,
(37)
Then, At can be reformulated as
(38)
where UA is of size lN×2N, VA is of size 2N×lN, and Ĥ is of size 2N×2N. By applying the Sherman–Morrison–Woodbury formula, the inverse matrix D̂UL,MA1 can be derived
(39)
Substitute Eq. (39) into Eq. (35), we have
(40)

Notice that D̂ULC,D̂UR, and Φf only depend on the spindle speed Ω, meaning that the matrices MPLD̂ULC1, MPLD̂ULC1UA, VAD̂ULC1, VAD̂ULC1UA, and D̂UR(ΦfO) can all be calculated only once outside the loop of the axial depth of cut a. Within the loop, the dimension of the matrix that needs to be inversed changes from lN to 2N, which is a significant reduction. The above improvement then leads to the following FCCM algorithm.

Table

Algorithm 3 FCCM

  forj=1 to NΩdo
  Ω=Ωj
  Calculate MPLD̂ULC1, MPLD̂ULC1UA, VAD̂ULC1UA, VAD̂ULC1, and D̂UR(ΦfO)
  fori=1 to Nado
   a=ai
   Calculate
  
   Calculate the spectral radius ρ(Φr)
  end for
 end for
  forj=1 to NΩdo
  Ω=Ωj
  Calculate MPLD̂ULC1, MPLD̂ULC1UA, VAD̂ULC1UA, VAD̂ULC1, and D̂UR(ΦfO)
  fori=1 to Nado
   a=ai
   Calculate
  
   Calculate the spectral radius ρ(Φr)
  end for
 end for

Now for every fixed spindle speed Ω, the algorithm is divided into two parts: outside the loop of the depth of cut a, a decomposition for the lN×lN matrix D̂ULC is performed once; while inside the loop, a decomposition of the 2N×2N matrix is performed for the construction of the transition matrix Φr of dimension l+2(N1) then following the calculation of the spectral radius ρ(Φr) for every a. We denote the time for the decomposition of D̂ULC as tinv,FCCM and the time spent on the spectral radius as teig,FCCM. The comparison of teig,FCCM and teig,FDM is shown in Fig. 6.

Fig. 6
teig,FCCM and teig,FDM: (a) value in parameter space (l,N) and (b) N=120
Fig. 6
teig,FCCM and teig,FDM: (a) value in parameter space (l,N) and (b) N=120
Close modal
By reducing the dimension of the matrix needs to be inversed inside the loop of a, the FCCM spends less time on the construction of the transition matrix Φ and calculating the spectral radius than FDM. However, due to the existence of the decomposition of D̂ULC outside the loop, it is not valid anymore to use the teig,FCCM and teig,FDM to replace the time of the construction of the whole SLD for comparison. Before we perform the comparison, the computation efficiency can be further improved by utilizing the special structure of A0 in D̂ULC. Notice that
(41)
Again, by allying the Sherman–Morrison–Woodbury formula, the inverse matrix D̂ULC1 has the following expression:
(42)

the derivation in detail is provided in Appendix  C. Now, the inversion of D̂ULC follows two steps: first, calculate D̂ULC,K1; next, calculate D̂ULC1. The main source of the computation cost changes from the decomposition of D̂ULC to that of (INM̂+D¯UL1Ĉ+D¯UL2K̂), which is only half the size of D̂ULC.

Define the time for decomposing the matrix D̂ULC using Eq. (42) as tinv,FCCM,r, the comparison of tinv,FCCM and tinv,FCCM,r is shown in Fig. 7. Because the dimension of the matrix needed to be decomposed is only half of D̂ULC, the decomposition time should be approximately a quarter of that of D̂ULC, i.e., tinv,FCCM,r14tinv,FCCM. From Fig. 7, however, the actual decomposition time is shorter than predicted, meaning a better efficiency improvement. The possible reason for this enhancement is the rearrangement of the elements in the sparse matrix (INM̂+D¯UL1Ĉ+D¯UL2K̂) (denoted as DULC,r̂ in the remaining of this article) compared to D̂ULC. The distribution of the elements of D̂ULC and D̂ULC,r is shown in Fig. 8. In D̂ULC, most blocks are identity matrices, however, along the main diagonal are the dense matrix A0, where almost all the elements are no-zero. On contrast, all the blocks of D̂ULC,r are the linear combination of the matrices M̂, Ĉ, and K̂, which retain the sparse band structure due to discretization using the finite element analysis. The density of the original A0 is in a sense “equally distributed” to the other blocks in D̂ULC,r, possibly easing the burden of the sparse matrix solver applied in matlab, resulting in a further improvement of the computation efficiency not only due to the reduction of the dimension.

Fig. 7
tinv,FCCM and tinv,FCCM,r: (a) value in parameter space (l,N) and (b) N=120
Fig. 7
tinv,FCCM and tinv,FCCM,r: (a) value in parameter space (l,N) and (b) N=120
Close modal
Fig. 8
Distribution of the nonzero elements in the sparse matrix D̂ULC and D̂ULC,r: (a) D̂ULC and (b) D̂ULC,r
Fig. 8
Distribution of the nonzero elements in the sparse matrix D̂ULC and D̂ULC,r: (a) D̂ULC and (b) D̂ULC,r
Close modal

After this dimension reduction, in Algorithm 3, we use Eq. (42) to calculate D̂ULC1, then Eq. (40) for Φr.

4.3 Convergence Analysis of Fast Chebyshev Collocation Method.

One of the important advantages of CCM is its spectral convergence rate, to demonstrate this superiority, the spectral radius of the transition matrix at the point (Ω,a)=(5000rpm,2mm) is calculated (other parameter remains the same as Sec. 3.2) under different system DOF l and number of the discretized segments N. The absolute error of the spectral radius and the reference value |ρ(Φ)ρ| is shown in Fig. 9. Where for every fixed l, ρ is the value of ρ(Φ) when N=160 (for the bigger value, the computation meets the memory usage problem, so in this article, we just regard it as the reference value).

Fig. 9
Comparison of the convergence rate between FDM and FCCM: (a) l=24, (b) l=84, (c) l=132, and (d) l=192
Fig. 9
Comparison of the convergence rate between FDM and FCCM: (a) l=24, (b) l=84, (c) l=132, and (d) l=192
Close modal

It is clear that when N is the same, the magnitude of the error of FCCM is often several orders higher than FDM. And when l=24, the FCCM shows the apparent spectral convergence rate: when N reaches 90, the error shows a rapid decrease with the exponential rate, when N=120, the error order reaches 1014, near the machine precise, and when N gets larger, this order remains unchanged. However, when l gets larger, this property seems to gradually disappear, and the order that the error can finally reach get lower, when l=192, the error order can only reach 107 at last. Two possible reasons may be able to explain this decay of the spectral convergence rate: one is that when l gets larger, the numeric error introduced by the decomposition of the sparse matrix of quite a huge dimension significantly disturbs and influences the numeric property of the collocation method; therefore, the spectral convergence property itself disappears under extremely high DOF; another possibility is that the spectral convergence property still exists even with very high DOF, but N used in this article is not large enough to observe it. But anyhow, FCCM is still proven to be superior to FDM in terms of computation accuracy.

4.4 Comparison Between Fast Chebyshev Collocation Method and Full-Discretization Method.

As is pointed out in Sec. 4.2, instead of teig,FCCM and teig,FDM, we define the new value, called tloop,FCCM and tloop,FDM to implement the comparison of computation efficiency, which means the time spent on one fixed spindle speed Ω when constructing the SLD
(43)
Na is the number of mesh points along the axis of the depth of cut a. The comparison of the computation efficiency of the newly proposed FCCM and FDM can be explored using these two variables. First, the results when Na=100 is shown in Fig. 10. In most cases, the overall computation time of FCCM successfully becomes shorter than FDM when l and N are equal. However, in the area where l and N are both large, tloop,FCCM exceeds tloop,FDM, meaning that FCCM is slower. The reason can possibly be explained as follows. From Eq. (43), Na is also an essential parameter that decides the relative size of tloop,FCCM and tloop,FDM. For a fixed (l,N), the relation between them and Na is shown in Fig. 11. Denote the value of Na when tloop,FCCM and tloop,FDM are equal as Na,critical
(44)

where [] denotes the smallest integer that is greater than or equal to the value “.” The FCCM is faster than FDM for the fixed (l,N) when Na exceeds Na,critical, meaning that FCCM is more computationally efficient only when the mesh of the SLD is sufficiently dense. The value of Na,critical under different (l,N) is shown in Fig. 11. Fortunately, for most of the (l,N), Na,critical is below 100 except for an abrupt peak at the up-right corner of the (l,N) space, where the maximum Na,critical reaches 382. By increasing Na, as shown in Fig. 12, the number of the points where FCCM is slower gets lower. In most situations, the FCCM shows superiority to FDM both in computation accuracy and efficiency even with high DOF.

Fig. 10
tloop,FDM and tloop,FCCM: (a) value in parameter space (l,N), (b) N=120, (c) the ratio ρt=tloop,FDM/tloop,FCCM in parameter space (l,N), and (d) the distribution in the (l,N) space whether FCCM is faster or not
Fig. 10
tloop,FDM and tloop,FCCM: (a) value in parameter space (l,N), (b) N=120, (c) the ratio ρt=tloop,FDM/tloop,FCCM in parameter space (l,N), and (d) the distribution in the (l,N) space whether FCCM is faster or not
Close modal
Fig. 11
tloop,FCCM and tloop,FDM as the function of Na: (a) the function and the critical point Na,critical and (b) Na,critical in space of (l,N)
Fig. 11
tloop,FCCM and tloop,FDM as the function of Na: (a) the function and the critical point Na,critical and (b) Na,critical in space of (l,N)
Close modal
Fig. 12
Distribution in the (l,N) space whether FCCM is faster or not with larger Na: (a) Na=150 and (b) Na=200
Fig. 12
Distribution in the (l,N) space whether FCCM is faster or not with larger Na: (a) Na=150 and (b) Na=200
Close modal

5 Concluding Remarks

A novel fast Chebyshev collocation method is proposed to perform the stability analysis of delay differential equations with linear, periodic coefficients of the system with high DOF. As an application, a robotic milling model allowing arbitrarily high DOF is used to validate the performance of the proposed method. By reducing the dimension of the transition matrix and acceleration of the inversion of the corresponding matrix component when constructing the transition matrix, the proposed method reaches high computation efficiency under high DOF. Numeric experiments are implemented to evaluate the computation efficiency and accuracy, and a comparison with the widely used FDM is performed. The following conclusions are drawn:

  1. In most cases considered in this article, the CCM is slower than FDM in terms of constructing the SLD. In the worst case, CCM is nearly 50 times slower than FDM, proving that CCM is not suitable for the stability analysis of systems with high DOF.

  2. In all the cases considered in this article, the error of the spectral radius of the transition matrix is smaller than FDM, usually several orders higher. FCCM retains the spectral convergent of the original CCM under low DOF (l=24). When DOF is higher, this spectral convergent seems to disappear. However, given that FCCM still has higher computation accuracy than FDM, it is a potential choice for the stability analysis under high system DOF.

  3. In most cases considered in this article, the FCCM performs better than the original CCM and the FDM in computation efficiency. The values tloop,FCCM and tloop,FDM which denote the time spent on constructing the SLD at a fixed spindle speed Ω are used as the comparing criteria. With different DOF l and the number of the discrete segments N, the ratio ρt=tloop,FDM/tloop,FCCM is below one except in several cases when l and N are both high (l120, N100). Increasing Na, the number of mesh points along the axis of the depth of cut a in the SLD does reduce the occurrence of cases where FCCM is slower than FDM.

Acknowledgment

The authors gratefully acknowledge the support received from the NSFC (12272218) and the Ministry of Science and Technology of China (2023YFF0713400).

Funding Data

  • NSFC (12272218; Funder ID: 10.13039/501100001809).

  • Ministry of Science and Technology of China (2023YFF0713400; Funder ID: 10.13039/501100002855).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix A: Derivation of Eq. (4)

Assuming that no external force is applied on the robot-spindle assembly, by applying the RMAC, the drive torque writes
(A1)
where J˙(q,q˙) denotes the derivative of J(q) to time, and v and ω represent the velocity and angular velocity of the end of the assembly, respectively. The subscript d means expected, i.e., the corresponding value when the system goes exactly on the expected trajectory. Suppose that the expected displacement of the end is (x,y,z)T, the direction cosine matrix of the end is A=(nT,sT,aT)T, n, s, a all refer to 3×1 vector, then ep and eo represent the error of position and orientation, respectively,
(A2)

according to Ref. [32], when the absolute value of the elements of eo is small, we approximately have e˙oωωd,e¨oω˙ω˙d.

Since (vω)=J(q)q˙ [31], deriving it on both sides
(A3)
Substitute Eqs. (A1) and (A3) into Eq. (1), we have
(A4)
Then, Eq. (A4) is transformed into the Cartesian space of the end by left-multiplying J(q)T
(A5)

Appendix B: Parameters Used to Derive Eq. (8)

Table 1

The modified D–H parameters of the robot

Number of jointd (mm)a (mm)α (deg)
181500
20350−90
308500
4820145−90
50090
61700−90
Number of jointd (mm)a (mm)α (deg)
181500
20350−90
308500
4820145−90
50090
61700−90
Table 2

Dynamic parameters of the robot

Number of jointMass m (kg)Center of mass (xc, yc, zc) (mm)Center inertia tensor (Ixx, Iyy, Izz) (kg m2)
1121.09(140.890, 15.739, −123.659)(2.577, 6.170, 6.725)
2182.99(401.195, −0.501, −248.760)(17.503, 1.640, 17.070)
3126.74(60.468, −4.014, 13.270)(5.235, 2.031, 4.968)
487.06(−3.556, 0.12, −169.963)(0.430, 0.430, 0.487)
539.72(0.059, 32.403, 0.013)(0.183, 0.120, 0.174)
644.25(−3.177, 0.801, 113.311)(0.247, 0.175, 0.136)
Number of jointMass m (kg)Center of mass (xc, yc, zc) (mm)Center inertia tensor (Ixx, Iyy, Izz) (kg m2)
1121.09(140.890, 15.739, −123.659)(2.577, 6.170, 6.725)
2182.99(401.195, −0.501, −248.760)(17.503, 1.640, 17.070)
3126.74(60.468, −4.014, 13.270)(5.235, 2.031, 4.968)
487.06(−3.556, 0.12, −169.963)(0.430, 0.430, 0.487)
539.72(0.059, 32.403, 0.013)(0.183, 0.120, 0.174)
644.25(−3.177, 0.801, 113.311)(0.247, 0.175, 0.136)
Table 3

Other parameters of the robot

QuantitiesSymbolValueUnits
Coordinate of the end of the spindle in the frame of the last link of the robotpend(−200, 0, 200)mm
The direction cosine matrix of the end frame relative to the frame of the last linkAend(001010100)
The RMAC control parameterkp108s−2
kd2kps−1
QuantitiesSymbolValueUnits
Coordinate of the end of the spindle in the frame of the last link of the robotpend(−200, 0, 200)mm
The direction cosine matrix of the end frame relative to the frame of the last linkAend(001010100)
The RMAC control parameterkp108s−2
kd2kps−1
Table 4

Parameters of the cutting tool and the milling force (the parameters in Tables 13 are coming from Ref. [35], parameters Kt,Kr in Table 4 are coming from Ref. [6], and the other parameters in Table 4 are coming from Ref. [36])

QuantitiesSymbolValueUnits
Number of the cutting teethNtool4
Diameter of the cutting toolDtool25mm
Length of the cutting toolL121mm
Young's modulus of elasticityE231.932GPa
Poisson rationμ0.324
Shear correction factork6(1+μ)7+6μ
Density of the toolρ8563.55kg/m3
The mass coefficient of the Rayleigh dampingαM35.372
The stiffness coefficient of the Rayleigh dampingαK2.061 × 10−10
The direction cosine matrix of the local frame of the element relative to the end frameT(001010100)
Tangential coefficient of the cutting forceKt1835MPa
Normal coefficient of the cutting forceKr1137MPa
QuantitiesSymbolValueUnits
Number of the cutting teethNtool4
Diameter of the cutting toolDtool25mm
Length of the cutting toolL121mm
Young's modulus of elasticityE231.932GPa
Poisson rationμ0.324
Shear correction factork6(1+μ)7+6μ
Density of the toolρ8563.55kg/m3
The mass coefficient of the Rayleigh dampingαM35.372
The stiffness coefficient of the Rayleigh dampingαK2.061 × 10−10
The direction cosine matrix of the local frame of the element relative to the end frameT(001010100)
Tangential coefficient of the cutting forceKt1835MPa
Normal coefficient of the cutting forceKr1137MPa

Appendix C: Derivation of Eq. (42)

According to Eqs. (36) and (41), we have
(C1)
By applying the Sherman–Morrison–Woodbury formula, the associated property of the Kronecker product, (AB)(CD)=(AC)(BD), and the inversion of the Kronecker product, (AB)1=A1B1
(C2)
Then
(C3)

References

1.
Merritt
,
H.
,
1965
, “
Theory of Self-Excited Machine-Tool Chatter: Contribution to Machine-Tool Chatter Research—1
,”
ASME J. Eng. Ind.
, 87(4), pp.
447
454
.10.1115/1.3670861
2.
Balachandran
,
B.
,
Kalmár-Nagy
,
T.
, and
Gilsinn
,
D. E.
,
2009
,
Delay Differential Equations
,
Springer
, Berlin, German.
3.
Yan
,
Y.
,
Xu
,
J.
, and
Wiercigroch
,
M.
,
2014
, “
Chatter in a Transverse Grinding Process
,”
J. Sound Vib.
,
333
(
3
), pp.
937
953
.10.1016/j.jsv.2013.09.039
4.
Zhang
,
Z.
,
Li
,
H.
,
Liu
,
X.
,
Zhang
,
W.
, and
Meng
,
G.
,
2018
, “
Chatter Mitigation for the Milling of Thin-Walled Workpiece
,”
Int. J. Mech. Sci.
,
138–139
, pp.
262
271
.10.1016/j.ijmecsci.2018.02.014
5.
Zhang
,
Z.
,
Luo
,
M.
,
Wu
,
B.
, and
Zhang
,
D.
,
2021
, “
Dynamic Modeling and Stability Prediction in Milling Process of Thin-Walled Workpiece With Multiple Structural Modes
,”
Proc. Inst. Mech. Eng., Part B: J. Eng. Manuf.
,
235
(
14
), pp.
2205
2218
.10.1177/0954405420933710
6.
Cordes
,
M.
,
Hintze
,
W.
, and
Altintas
,
Y.
,
2019
, “
Chatter Stability in Robotic Milling
,”
Rob. Comput.-Integr. Manuf.
,
55
, pp.
11
18
.10.1016/j.rcim.2018.07.004
7.
Liu
,
X.
,
Vlajic
,
N.
,
Long
,
X.
,
Meng
,
G.
, and
Balachandran
,
B.
,
2014
, “
State Dependent Delay Influenced Drill-String Oscillations and Stability Analysis
,”
ASME J. Vib. Acoust.
,
136
(
5
), p.
051008
.10.1115/1.4027958
8.
Liu
,
X.
,
Long
,
X.
,
Zheng
,
X.
,
Meng
,
G.
, and
Balachandran
,
B.
,
2020
, “
Spatial Temporal Dynamics of a Drill String With Complex Time-Delay Effects: Bit Bounce and Stick-Slip Oscillations
,”
Int. J. Mech. Sci.
,
170
, p.
105338
.10.1016/j.ijmecsci.2019.105338
9.
Gu
,
K.
, and
Zheng
,
X.
,
2014
, “
An Overview of Stability Crossing Set for Systems With Scalar Delay Channels
,”
Proceedings of the 33rd Chinese Control Conference
, Nanjing, China, July 28–30, pp.
2676
2681
.10.1109/ChiCC.2014.6897059
10.
Sadath
,
A.
, and
Vyasarayani
,
C.
,
2015
, “
Galerkin Approximations for Stability of Delay Differential Equations With Time Periodic Coefficients
,”
ASME J. Comput. Nonlinear Dyn.
,
10
(
2
), p.
021011
.10.1115/1.4026989
11.
Altintaş
,
Y.
, and
Budak
,
E.
,
1995
, “
Analytical Prediction of Stability Lobes in Milling
,”
CIRP Ann.
,
44
(
1
), pp.
357
362
.10.1016/S0007-8506(07)62342-7
12.
Merdol
,
S.
, and
Altintas
,
Y.
,
2004
, “
Multi Frequency Solution of Chatter Stability for Low Immersion Milling
,”
ASME J. Manuf. Sci. Eng.
,
126
(
3
), pp.
459
466
.10.1115/1.1765139
13.
Insperger
,
T.
, and
Stépán
,
G.
,
2002
, “
Semi-Discretization Method for Delayed Systems
,”
Int. J. Numer. Methods Eng.
,
55
(
5
), pp.
503
518
.10.1002/nme.505
14.
Insperger
,
T.
, and
Stépán
,
G.
,
2004
, “
Updated Semi-Discretization Method for Periodic Delay-Differential Equations With Discrete Delay
,”
Int. J. Numer. Methods Eng.
,
61
(
1
), pp.
117
141
.10.1002/nme.1061
15.
Insperger
,
T.
,
Stépán
,
G.
, and
Turi
,
J.
,
2008
, “
On the Higher-Order Semi Discretizations for Periodic Delayed Systems
,”
J. Sound Vib.
,
313
(
1–2
), pp.
334
341
.10.1016/j.jsv.2007.11.040
16.
Ding
,
Y.
,
Zhu
,
L.
,
Zhang
,
X.
, and
Ding
,
H.
,
2010
, “
A Full-Discretization Method for Prediction of Milling Stability
,”
Int. J. Mach. Tools Manuf.
,
50
(
5
), pp.
502
509
.10.1016/j.ijmachtools.2010.01.003
17.
Jiang
,
S.
,
Sun
,
Y.
,
Yuan
,
X.
, and
Liu
,
W.
,
2017
, “
A Second-Order Semi Discretization Method for the Efficient and Accurate Stability Prediction of Milling Process
,”
Int. J. Adv. Manuf. Technol.
,
92
(
1–4
), pp.
583
595
.10.1007/s00170-017-0171-y
18.
Ding
,
Y.
,
Zhu
,
L.
,
Zhang
,
X.
, and
Ding
,
H.
,
2010
, “
Second-Order Full Discretization Method for Milling Stability Prediction
,”
Int. J. Mach. Tools Manuf.
,
50
(
10
), pp.
926
932
.10.1016/j.ijmachtools.2010.05.005
19.
Henninger
,
C.
, and
Eberhard
,
P.
,
2008
, “
Improving the Computational Efficiency and Accuracy of the Semi-Discretization Method for Periodic Delay-Differential Equations
,”
Eur. J. Mech.-A/Solids
,
27
(
6
), pp.
975
985
.10.1016/j.euromechsol.2008.01.006
20.
Mann
,
B. P.
,
Young
,
K. A.
,
Schmitz
,
T. L.
, and
Dilley
,
D. N.
,
2005
, “
Simultaneous Stability and Surface Location Error Predictions in Milling
,”
ASME J. Manuf. Sci. Eng.
, 127(3), pp.
446
453
.10.1115/1.1948394
21.
Khasawneh
,
F. A.
, and
Mann
,
B. P.
,
2011
, “
A Spectral Element Approach for the Stability of Delay Systems
,”
Int. J. Numer. Methods Eng.
,
87
(
6
), pp.
566
592
.10.1002/nme.3122
22.
Butcher
,
E. A.
,
Bobrenkov
,
O. A.
,
Bueler
,
E.
, and
Nindujarla
,
P.
,
2009
, “
Analysis of Milling Stability by the Chebyshev Collocation Method: Algorithm and Optimal Stable Immersion Levels
,”
ASME J. Comput. Nonlinear Dyn.
,
4
(
3
), p.
031003
.10.1115/1.3124088
23.
Breda
,
D.
,
Maset
,
S.
, and
Vermiglio
,
R.
,
2012
, “
Approximation of Eigenvalues of Evolution Operators for Linear Retarded Functional Differential Equations
,”
SIAM J. Numer. Anal.
,
50
(
3
), pp.
1456
1483
.10.1137/100815505
24.
Breda
,
D.
,
Maset
,
S.
, and
Vermiglio
,
R.
,
2014
, “
Pseudospectral Methods for Stability Analysis of Delayed Dynamical Systems
,”
Int. J. Dyn. Control
,
2
(
2
), pp.
143
153
.10.1007/s40435-013-0041-x
25.
Khasawneh
,
F. A.
,
Mann
,
B. P.
, and
Butcher
,
E. A.
,
2011
, “
A Multi-Interval Chebyshev Collocation Approach for the Stability of Periodic Delay Systems With Discontinuities
,”
Commun. Nonlinear Sci. Numer. Simul.
,
16
(
11
), pp.
4408
4421
.10.1016/j.cnsns.2011.03.025
26.
Tweten
,
D. J.
,
Lipp
,
G. M.
,
Khasawneh
,
F. A.
, and
Mann
,
B. P.
,
2012
, “
On the Comparison of Semi-Analytical Methods for the Stability Analysis of Delay Differential Equations
,”
J. Sound Vib.
,
331
(
17
), pp.
4057
4071
.10.1016/j.jsv.2012.04.009
27.
Khasawneh
,
F. A.
,
Mann
,
B. P.
, and
Butcher
,
E. A.
,
2010
, “
Comparison Between Collocation Methods and Spectral Element Approach for the Stability of Periodic Delay Systems
,”
IFAC Proc. Vol.
,
43
(
2
), pp.
69
74
.10.3182/20100607-3-CZ-4010.00014
28.
Wan
,
M.
,
Ma
,
Y.-C.
,
Zhang
,
W.-H.
, and
Yang
,
Y.
,
2015
, “
Study on the Construction Mechanism of Stability Lobes in Milling Process With Multiple Modes
,”
Int. J. Adv. Manuf. Technol.
,
79
(
1–4
), pp.
589
603
.10.1007/s00170-015-6829-4
29.
Totis
,
G.
,
Albertelli
,
P.
,
Sortino
,
M.
, and
Monno
,
M.
,
2014
, “
Efficient Evaluation of Process Stability in Milling With Spindle Speed Variation by Using the Chebyshev Collocation Method
,”
J. Sound Vib.
,
333
(
3
), pp.
646
668
.10.1016/j.jsv.2013.09.043
30.
Totis
,
G.
,
Insperger
,
T.
,
Sortino
,
M.
, and
Stépán
,
G.
,
2019
, “
Symmetry Breaking in Milling Dynamics
,”
Int. J. Mach. Tools Manuf.
,
139
, pp.
37
59
.10.1016/j.ijmachtools.2019.01.002
31.
Spong
,
M. W.
, and
Vidyasagar
,
M.
,
2008
,
Robot Dynamics and Control
,
Wiley
, Hoboken, NJ.
32.
Luh
,
J.
,
Walker
,
M.
, and
Paul
,
R.
,
1980
, “
Resolved-Acceleration Control of Mechanical Manipulators
,”
IEEE Trans. Autom. Control
,
25
(
3
), pp.
468
474
.10.1109/TAC.1980.1102367
33.
Altintas
,
Y.
, and
Ber
,
A.
,
2001
, “
Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design
,”
ASME Appl. Mech. Rev.
,
54
(
5
), p.
B84
.10.1115/1.1399383
34.
Faassen
,
R.
,
Van de Wouw
,
N.
,
Nijmeijer
,
H.
, and
Oosterling
,
J.
,
2007
, “
An Improved Tool Path Model Including Periodic Delay for Chatter Prediction in Milling
,”
ASME J. Comput. Nonlinear Dynam.
, 2(2), pp.
167
179
.10.1115/1.2447465
35.
Wu
,
J.
,
2020
, “
Research on the Stiffness and Milling Chatter of Industrial Robot in Different Configurations
,” Master's thesis,
Jilin University
, Jilin, China.
36.
Li
,
Y.
,
2018
, “
Rapid Dynamics Prediction of Tool Tip and Analysis of the Chatter Stability in Robotic Milling System
,” Master's thesis,
Huazhong University of Science and Technology
, Wuhan, China.