A Fast Chebyshev Collocation Method for Stability Analysis of a Robotic Machining System With Time Delay

Li, Chenglin; Luo, Qingpeng; Meng, Guang; Liu, Xianbo

doi:10.1115/1.4067062

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 [3–5]. 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

View large Download slide

The robotic milling system

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:

M (q) \ddot{q} + C (q, \dot{q}) + G (q) = τ_{d}

(1)

where

q, \dot{q}, \ddot{q}

are vectors of the size

6 \times 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, \dot{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 (⁠

X_{1} - Y_{1} - Z_{1}

in Fig. 1) as

J {(q)}^{- T} M (q) J {(q)}^{- 1} ((\begin{matrix} {\ddot{e}}_{p} \\ {\ddot{e}}_{o} \end{matrix}) + k_{d} (\begin{matrix} {\dot{e}}_{p} \\ {\dot{e}}_{o} \end{matrix}) + k_{p} (\begin{matrix} e_{p} \\ e_{o} \end{matrix})) = 0

(2)

where $J (q)$ is the Jacobi matrix of size $6 \times 6$ mapping the angular velocity of the joints into the velocity and angular velocity in the end coordinate frame [31]. $e_{p}, e_{o}$ denote the error between the actual and expected displacement and rotation angle in the end coordinate frame, respectively: $e_{p} = {(x - x_{d}, y - y_{d}, z - z_{d})}^{T}$ ⁠, $e_{o} = {(α - α_{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 $k_{d}, k_{p}$ 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

M_{tool} \ddot{X} + C_{tool} \dot{X} + K_{tool} X = F

(3)

where $M_{tool}$ ⁠, $C_{tool}$ ⁠, and $K_{tool}$ 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 $M_{tool}$ and $K_{tool}$ ⁠: $C_{tool} = α_{M} M_{tool} + α_{K} K_{tool}$ ⁠, 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 = {(\begin{matrix} e_{p, 1}^{T} & e_{o, 1}^{T} & \dots & e_{p, n + 1}^{T} & e_{o, n + 1}^{T} \end{matrix})}^{T}$ ⁠, where $e_{p, i}, e_{o, 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 (⁠ $X_{1} Y_{1} Z_{1}$ in Fig. 1, where the spindle and the tool are connected), subscript “ $n + 1$ ” denotes the node at the tool tip (⁠ $X_{n + 1} Y_{n + 1} Z_{n + 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

\begin{matrix} \hat{M} \ddot{X} + \hat{C} \dot{X} + \hat{K} X = F \\ \hat{M} = M_{tool} + (\begin{matrix} J {(q)}^{- T} M (q) J {(q)}_{6 \times 6}^{- 1} & O \\ O & O_{6 n \times 6 n} \end{matrix}) \\ \hat{C} = C_{tool} + (\begin{matrix} k_{d} J {(q)}^{- T} M (q) J {(q)}_{6 \times 6}^{- 1} & O \\ O & O_{6 n \times 6 n} \end{matrix}) \\ \hat{K} = K_{tool} + (\begin{matrix} k_{p} J {(q)}^{- T} M (q) J {(q)}_{6 \times 6}^{- 1} & O \\ O & O_{6 n \times 6 n} \end{matrix}) \end{matrix}

(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 $\hat{M}, \hat{C}$ ⁠, and $\hat{K}$ remain to be constant. Here, assuming that the two-dimensional cutting model is applied, then F is formulated as: $F = {(\begin{matrix} 0_{1 \times 6 n} & F_{x} & F_{y} & 0_{1 \times 4} \end{matrix})}^{T}$ ⁠, where $F_{x}, F_{y}$ 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 $X_{n + 1} Y_{n + 1} Z_{n + 1}$ ⁠, i.e., the tool coordinate frame, is fixed on the tool tip, i.e., the node $n + 1$ ⁠, where $X_{n + 1}, Y_{n + 1}, Z_{n + 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, $a_{e}$ denotes the radial depth of cut, and v denotes the feed velocity of the workpiece relative to the cutting tool. $D_{tool}$ 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 $N_{tool}$ equally distributed cutting teeth, resulting in a constant pitch angle, then $φ_{j} (t) = (2 π Ω / 60) t + (j - 1) (2 π / N_{tool})$ ⁠.

Fig. 2

View large Download slide

The two-dimensional cutting force model: (a) view from the $Z_{n + 1}$ axis and (b) the three-dimensional view

$F_{t, j}, F_{r, j}$ represent the tangential and normal cutting force applied on the jth tooth, and $h_{j}$ 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: $F_{t, j} = a K_{t} h_{j}$ ⁠, $F_{r, j} = a K_{r} h_{j}$ ⁠, and $h_{j} = (x (t) - x (t - T)) sin φ_{j} + (y (t) - y (t - T)) cos φ_{j}$ ⁠, where $K_{t}$ and $K_{r}$ are tangential and normal cutting force coefficients, respectively. T denotes the tooth passing period: $T = 60 / (N_{tool} Ω)$ ⁠, which is also the delay of the dynamic system.

According to the previous assumptions, we have

x_{d} (t) = x_{0} + v t, y_{d} (t) = y_{0}

⁠, where

x_{0}

and

y_{0}

are the initial position of the tool tip in

X_{n + 1}

and

Y_{n + 1}

direction. Then,

x (t) - x (t - T) = e_{p, x} (t) - e_{p, x} (t - T) + f

and

y (t) - y (t - T) = e_{p, y} (t) - e_{p, y} (t - T)

⁠, where

e_{p, x}

and

e_{p, y}

are the terms

x - x_{d}

and

y - y_{d}

in

e_{p}

⁠, respectively. f denotes the feed per tooth, i.e.,

f = v T = 60 v / (N_{tool} Ω)

⁠. 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

X_{n + 1}

and

Y_{n + 1}

direction, the two-dimensional cutting force is represented as

\begin{matrix} \begin{matrix} (\begin{matrix} F_{x} (t) \\ F_{y} (t) \end{matrix}) = H {(t)}_{2 \times 2} (\begin{matrix} e_{p, x} (t) - e_{p, x} (t - T) + f \\ e_{p, y} (t) - e_{p, y} (t - T) \end{matrix}) \end{matrix} \\ H (t) = a \sum_{j = 1}^{N} g (φ_{j}) (\begin{matrix} - cos φ_{j} & - sin φ_{j} \\ sin φ_{j} & - cos φ_{j} \end{matrix}) (\begin{matrix} K_{t} \\ K_{r} \end{matrix}) (\begin{matrix} sin φ_{j} & cos φ_{j} \end{matrix}) \end{matrix}

(5)

where

H (t)

is a time-varying periodic matrix of size

2 \times 2

⁠, i.e.,

H (t) = H (t + T)

⁠. Window function g is a switch function determining whether the tooth is in cutting

\begin{matrix} g (φ) = {\begin{matrix} 1 & φ_{s t} < φ (mod 2 π) \leq φ_{e x} \\ 0 & otherwise \end{matrix} \end{matrix}

(6)

φ_{s t}

and

φ_{e x}

are the entrance and exit angles of the cutting tooth, respectively,

\begin{matrix} (φ_{s t}, φ_{e x}) = {\begin{matrix} (0, \arccos (1 - 2 \frac{a_{e}}{D_{tool}})), & up milling \\ (\arccos (2 \frac{a_{e}}{D_{tool}} - 1), π), & down milling \end{matrix} \end{matrix}

(7)

Then, the final DDE of the robotic milling system is formulated as

\begin{matrix} \begin{matrix} \hat{M} \ddot{X} (t) + \hat{C} \dot{X} (t) + \hat{K} X (t) = {\hat{H}}_{1} (t) X (t) + {\hat{H}}_{2} (t) X (t - τ) + \hat{f} (t), \hat{f} (t) = (\begin{matrix} O_{6 n \times 1} \\ H (t) (\begin{matrix} f \\ 0 \end{matrix}) \\ O_{4 \times 1} \end{matrix}) \end{matrix} \\ \begin{matrix} {\hat{H}}_{1} (t) = (\begin{matrix} O_{6 n \times 6 n} \\ H (t) \\ O_{4 \times 4} \end{matrix}), {\hat{H}}_{2} (t) = (\begin{matrix} O_{6 n \times 6 n} \\ - H (t) \\ O_{4 \times 4} \end{matrix}) \end{matrix} \end{matrix}

(8)

In Eq. (8), $\hat{M}, \hat{C}$ ⁠, and $\hat{K}$ are constant matrices of size $6 (n + 1) \times 6 (n + 1)$ ⁠, ${\hat{H}}_{1} (t)$ and ${\hat{H}}_{2} (t)$ are periodic matrices of size $6 (n + 1) \times 6 (n + 1)$ ⁠: ${\hat{H}}_{1} (t) = {\hat{H}}_{1} (t + T), {\hat{H}}_{2} (t) = {\hat{H}}_{2} (t + T)$ ⁠, and $\hat{f} (t)$ is periodic vector of length $6 (n + 1)$ ⁠: $\hat{f} (t) = \hat{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. is performed. Let

Z (t) = {(X^{T} (t) {(\hat{M} \dot{X} (t) + \frac{1}{2} \hat{C} X (t))}^{T})}^{T}

⁠, Eq. is rewritten as the first-order state-space form

(9)

where

A_{0}

is a constant matrix, and

A (t) = A (t + T), B (t) = B (t + T)

⁠, and

F_{0} (t) = F_{0} (t + T)

⁠. Therefore, Eq. has a periodic solution

\bar{Z} (t)

⁠, i.e.,

\bar{Z} (t) = \bar{Z} (t + T)

⁠. Let

Z (t) = \bar{Z} (t) + \hat{Z} (t)

⁠, where

\hat{Z} (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

\bar{Z} (t)

satisfies

\dot{\bar{Z}} (t) = A_{0} \bar{Z} (t) + F_{0} (t)

(10)

which is an ODE. Combine Eqs. and , we have

\dot{\hat{Z}} (t) = A_{0} \hat{Z} (t) + A (t) \hat{Z} (t) + B (t) \hat{Z} (t - τ)

(11)

which is a DDE only corresponding to the perturbation $\hat{Z} (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

\hat{Z} (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. , 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

[t_{s t}, t_{s t} + T]

⁠, where

t_{0}

is the time when the angular position of one certain cutting tooth is equal to the entrance angle:

φ_{j} (t_{s t}) = φ_{s t}

⁠. This shifting ensures that there is at most one nonsmooth point inside the interval

[t_{s t}, t_{s t} + T]

⁠.

t_{s t}

has infinite possible values, and the one with the smallest absolute value is chosen

t_{s t} = \frac{φ_{s t} - φ_{j_{s t}} (0)}{2 π / N_{tool}} T, j_{s t} = \underset{1 \leq j \leq N_{tool}}{\arg \min} | φ_{s t} - φ_{j} (0) |

(12)

where $j_{s t}$ is the index of the cutting tooth such that the absolute value $| φ_{s t} - φ_{j} (0) |$ reaches the minimum. Then, let $ρ_{cut} = (φ_{e x} - φ_{s t}) / (2 π / N_{tool})$ ⁠, 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 $(t_{s t}, t_{s t} + 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

\hat{Z} (t)

is approximated changes from

(t_{s t}, t_{s t} + T]

to

(t_{s t}, t_{s t} + ρ_{cut} T]

⁠, using the

N + 1

order Chebyshev points of the second kind on [a, b] as the interpolation points:

x_{i} = ((b + a) / 2) + ((b - a) / 2) cos (i π / N)

⁠,

i = 0, 1, \dots, N

⁠, the solution

\hat{Z} (t)

is approximated as

\begin{matrix} \hat{Z} (t) \approx \sum_{i = 0}^{N} \hat{Z} (t_{i}) ϕ_{i} (t) \\ t_{i} = \frac{2 t_{s t} + ρ_{cut} T}{2} + \frac{ρ_{cut} T}{2} cos \frac{i π}{N}, i = 0, 1, \dots, N \end{matrix}

(13)

ϕ_{i} (t)

is the ith Lagrange basis polynomial corresponding to point

t_{i}

ϕ_{i} (t) = \prod_{j = 0, j \neq i}^{N} \frac{(t - t_{0}) \dots (t - t_{N})}{(t_{i} - t_{0}) \dots (t_{i} - t_{N})}

(14)

Then, substitute Eq. into Eq. , we get a series of equations, which is the discrete form of the approximation of Eq.

\sum_{j = 0}^{N} \hat{Z} (t_{j}) {\dot{ϕ}}_{j} (t_{i}) = (A_{0} + A (t_{i})) \hat{Z} (t_{i}) + B (t_{i}) \hat{Z} (t_{i} - τ), i = 0, 1, \dots, N

(15)

then rewrite them as the matrix form

\begin{matrix} D (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) = M_{A} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) + M_{B} (\begin{matrix} \hat{Z} (t_{0} - τ) \\ \hat{Z} (t_{1} - τ) \\ ⋮ \\ \hat{Z} (t_{N} - τ) \end{matrix}) \end{matrix}

(16)

where

D = (\begin{matrix} {\dot{ϕ}}_{0} (t_{0}) & \dots & {\dot{ϕ}}_{N} (t_{0}) \\ ⋮ & ⋱ & ⋮ \\ {\dot{ϕ}}_{0} (t_{N}) & \dots & {\dot{ϕ}}_{N} (t_{N}) \end{matrix}) \otimes I_{l}

(17)

here, the operation

\otimes

is the Kronecker product.

I_{l}

is the identity matrix of size

l \times l

⁠, and l is the length of vector

\hat{Z} (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”:

\begin{matrix} D = \frac{2}{ρ_{cut} T} \bar{D} \otimes I_{l} \\ {\bar{D}}_{i j} = \frac{c_{i}}{c_{j}} \frac{{(- 1)}^{i + j}}{cos \frac{i π}{N} - cos \frac{j π}{N}}, i, j = 0, 1, \dots, N, i \neq j \\ {\bar{D}}_{i i} = \sum_{j = 0, j \neq i}^{N} {\bar{D}}_{i j}, i = 0, 1, \dots, N c_{i} = {\begin{matrix} 2, & i = 0, N \\ 1, & otherwise \end{matrix} \end{matrix}

(18)

In addition, due to the existence of the free state, the relationship of the endpoint, i.e.,

\hat{Z} (t_{N})

and

\hat{Z} (t_{0} - τ)

⁠, becomes

\hat{Z} (t_{N}) = e^{A_{0} (1 - ρ_{cut}) T} \hat{Z} (t_{0} - τ)

[22]. To enforce this condition, the last l rows of the matrices

D, M_{A}

⁠, and

M_{B}

need to be modified, resulting in the following equations:

\begin{matrix} \hat{D} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) = {\hat{M}}_{A} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) + {\hat{M}}_{B} (\begin{matrix} \hat{Z} (t_{0} - τ) \\ \hat{Z} (t_{1} - τ) \\ ⋮ \\ \hat{Z} (t_{N} - τ) \end{matrix}) \end{matrix}

(19)

where

(20)

and

\begin{matrix} \hat{D} = (\begin{matrix} {\bar{D}}_{U L} \otimes I_{l} & {\bar{D}}_{U R} \otimes I_{l} \\ O & I_{l} \end{matrix}) \end{matrix}

(21)

where ${\bar{D}}_{U L}$ denotes the block of $(2 / (ρ_{cut} T)) \bar{D}$ including the first N rows and columns, while ${\bar{D}}_{U R}$ denotes the block of $(2 / (ρ_{cut} T)) \bar{D}$ including the first N rows and the last column. Notice that $\hat{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. , the discrete map mapping state terms within

[- T, 0]

to [0, T] is immediately obtained

\begin{matrix} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) = Φ (\begin{matrix} \hat{Z} (t_{0} - τ) \\ \hat{Z} (t_{1} - τ) \\ ⋮ \\ \hat{Z} (t_{N} - τ) \end{matrix}) \end{matrix}

(22)

where

Φ = {(\hat{D} - {\hat{M}}_{A})}^{- 1} {\hat{M}}_{B}

(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. . Based on the Floquet theory, the stability of the DDE can be determined according to the spectral radius of the transition matrix

ρ (Φ)

\begin{matrix} ρ (Φ) {\begin{matrix} < 1 & stable \\ = 1 & critical stable \\ > 1 & instable \end{matrix} \end{matrix}

(24)

In practical application, the inversion of the matrix ${(\hat{D} - {\hat{M}}_{A})}^{- 1}$ can be avoided by solving the generalized eigenvalue problem, i.e., finding $λ$ such that the dominant $| {\hat{M}}_{B} - λ (\hat{D} - {\hat{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

[a_{min}, a_{max}]

⁠, then the parameter space

(Ω, a)

is discretized into

N_{Ω} \times N_{a}

mesh points

(Ω_{j}, a_{i})

with constant increment

d Ω

and

d a

⁠, respectively,

\begin{matrix} a_{i} = a_{min} + (i - 1) d a, i = 1, 2, \dots, N_{a} \\ Ω_{j} = Ω_{min} + (j - 1) d Ω, i = 1, 2, \dots, N_{Ω} \end{matrix}

(25)

$N_{a}$ 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

for

j = 1

to

N_{Ω}

do

Ω = Ω_{j}

Calculate

\hat{D}

for

i = 1

to

N_{a}

do

a = a_{i}

Calculate

{\hat{M}}_{A}

⁠,

{\hat{M}}_{B}

Solve the generalized eigenvalue problem

| {\hat{M}}_{B} - λ (\hat{D} - {\hat{M}}_{A}) | = 0

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}, a_{i})$ when using CCM and FDM as $t_{eig, CCM}$ and $t_{eig, FDM}$ ⁠, respectively. Assuming that $t_{eig, CCM}$ and $t_{eig, FDM}$ remain the same regardless of the $Ω_{j}$ and $a_{i}$ ⁠, we do not use the time spent on constructing the whole SLD to assess the computation efficiency of CCM and FDM but $t_{eig, CCM}$ and $t_{eig, FDM}$ instead.

The DDE model built in Sec. 2 is adopted to apply the algorithm. The other necessary parameters are $Ω = 500 rpm, a = 8 mm, a_{e} / D_{tool} = 0.2, and N_{tool} = 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 $t_{eig, CCM}$ and $t_{eig, 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} = t_{eig, FDM} / t_{eig, 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 + 2 N$ ⁠, while in CCM it is $l (N + 1)$ ⁠. When the DOF of the system l is high, i.e., $l ≫ 2$ ⁠, then $l (N + 1) ≫ l + 2 N$ ⁠, resulting in a huge computation cost when solving the (generalized) eigenvalue problem compared to the FDM.

Fig. 3

View large Download slide

$t_{eig, CCM}$ and $t_{eig, FDM}$ ⁠: (a) value in parameter space $(l, N)$ and (b) $N = 60$

Fig. 4

View large Download slide

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

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.1–4.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. are almost zeros except for the

2 \times 2

block

H (t)

⁠; therefore, the term

B (t) \hat{Z} (t - τ)

depends only on two elements of the delay term

\hat{Z} (t - τ)

⁠. Then, Eq. can be rewritten as the following form:

\begin{matrix} \dot{\hat{Z}} (t) = A_{0} \hat{Z} (t) + A (t) \hat{Z} (t) + \hat{B} (t) U (t - τ) \\ U (t) = P Z (t) \end{matrix}

(26)

where

(27)

The matrix P can be regarded as a selection matrix choosing the elements of

\hat{Z} (t - τ)

that really contribute to the value of

B (t) \hat{Z} (t - τ)

⁠. Noted that

B (t) = \hat{B} (t) P

⁠, Eq. can be reduced to the following form:

\begin{matrix} \hat{D} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) & = {\hat{M}}_{A} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) + {\hat{M}}_{B}^{r} (\begin{matrix} \hat{Z} (t_{0} - τ) \\ U (t_{1} - τ) \\ ⋮ \\ U (t_{N - 1} - τ) \end{matrix}) \\ (\begin{matrix} \hat{Z} (t_{0}) \\ U (t_{1}) \\ ⋮ \\ U (t_{N - 1}) \end{matrix}) & = M_{P} (\begin{matrix} \hat{Z} (t_{0}) \\ \hat{Z} (t_{1}) \\ ⋮ \\ \hat{Z} (t_{N}) \end{matrix}) \end{matrix}

(28)

where

{\hat{M}}_{B}^{r}

is a

l (N + 1) \times (l + 2 (N - 1))

size matrix, while

M_{P}

is a

(l + 2 (N - 1)) \times l (N + 1)

size matrix

\begin{matrix} {\hat{M}}_{B}^{r} = {(\begin{matrix} \hat{B} (t_{0}) P \\ \hat{B} (t_{1}) \\ ⋱ \\ \hat{B} (t_{N - 1}) \\ Φ_{f} & O & \dots & O \end{matrix})}_{l (N + 1) \times (l + 2 (N - 1))} \\ M_{P} = {(\begin{matrix} I_{l} & O \\ P & O \\ ⋱ & ⋮ \\ P & O \end{matrix})}_{(l + 2 (N - 1)) \times l (N + 1)} \end{matrix}

(29)

Then, the reduced discrete map can be obtained

\begin{matrix} (\begin{matrix} \hat{Z} (t_{0}) \\ U (t_{1}) \\ ⋮ \\ U (t_{N - 1}) \end{matrix}) = Φ_{r} (\begin{matrix} \hat{Z} (t_{0} - τ) \\ U (t_{1} - τ) \\ ⋮ \\ U (t_{N - 1} - τ) \end{matrix}) \end{matrix}

(30)

where the reduced transition matrix

Φ_{r}

can be formulated as

Φ_{r} = M_{P} {(\hat{D} - {\hat{M}}_{A})}^{- 1} {\hat{M}}_{B}^{r}

(31)

Now $Φ_{r}$ is an $(l + 2 (N - 1)) \times (l + 2 (N - 1))$ 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}$ ⁠)

for

j = 1

to

N_{Ω}

do

Ω = Ω_{j}

Calculate

\hat{D}

for

i = 1

to

N_{a}

do

a = a_{i}

Calculate

Φ_{r} = M_{P} {(\hat{D} - {\hat{M}}_{A})}^{- 1} {\hat{M}}_{B}^{r}

Calculate the spectral radius

ρ (Φ_{r})

end for

Define the time spend on constructing $Φ_{r}$ and calculating $ρ (Φ_{r})$ for a single parameter pair $(Ω, a)$ is $t_{eig, CCM, r}$ ⁠, similar to Sec. 3.2, then the results of $t_{eig, CCM}$ and $t_{eig, 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) \times l (N + 1)$ matrix ${(\hat{D} - {\hat{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

View large Download slide

$t_{eig, CCM}$ and $t_{eig, CCM, r}$ ⁠: (a) value in parameter space $(l, N)$ and (b) $N = 60$

4.2 Acceleration of Matrix Inversion.

To simplify the matrix inversion, considering the matrix

\hat{D} - {\hat{M}}_{A}

⁠, it has the following partition:

(32)

then

(33)

Appropriately blocking the matrix

{\hat{M}}_{B}^{r}

and

M_{P}

(34)

Then,

Φ_{r}

can be formulated as

(35)

Now we need to simplify the calculation of

{\hat{D}}_{U L, M_{A}}^{- 1}

⁠. The matrix

{\hat{D}}_{U L, M_{A}}

can be separated into two parts

{\hat{D}}_{U L, M_{A}} = \underset{{\hat{D}}_{ULC}}{\underset{︸}{{\bar{D}}_{U L} \otimes I_{l} - I_{N} \otimes A_{0}}} - A_{t} = {\hat{D}}_{ULC} - A_{t}

(36)

{\hat{D}}_{ULC}

is of size

l N \times l N

and only depends on the spindle speed

Ω

⁠. As for

A_{t}

⁠, similarly, noted that the elements of the coefficient matrix

A (t)

in Eq. are almost zeros except for the

2 \times 2

block

H (t)

⁠, i.e.,

(37)

Then,

A_{t}

can be reformulated as

\begin{matrix} A_{t} = \underset{U_{A}}{\underset{︸}{(I_{N} \otimes U_{H})}} \underset{\hat{H}}{\underset{︸}{(\begin{matrix} H (t_{0}) \\ ⋱ \\ H (t_{N - 1}) \end{matrix})}} \underset{V_{A}}{\underset{︸}{(I_{N} \otimes V_{H})}} = U_{A} \hat{H} V_{A} \end{matrix}

(38)

where

U_{A}

is of size

l N \times 2 N

⁠,

V_{A}

is of size

2 N \times l N

⁠, and

\hat{H}

is of size

2 N \times 2 N

⁠. By applying the Sherman–Morrison–Woodbury formula, the inverse matrix

{\hat{D}}_{U L, M_{A}}^{- 1}

can be derived

\begin{matrix} {\hat{D}}_{U L, M_{A}}^{- 1} = {({\hat{D}}_{ULC} - U_{A} \hat{H} V_{A})}^{- 1} \\ \begin{matrix} = {\hat{D}}_{ULC}^{- 1} + {\hat{D}}_{ULC}^{- 1} U_{A} {\hat{H}}^{- 1} {(I_{2 N} - V_{A} {\hat{D}}_{ULC}^{- 1} U_{A} \hat{H})}^{- 1} V_{A} {\hat{D}}_{ULC}^{- 1} \end{matrix} \end{matrix}

(39)

Substitute Eq. into Eq. , we have

\begin{matrix} Φ_{r} = (M_{P L} {\hat{D}}_{ULC}^{- 1} + M_{P L} {\hat{D}}_{ULC}^{- 1} U_{A} \hat{H} {(I_{2 N} - V_{A} {\hat{D}}_{ULC}^{- 1} U_{A} \hat{H})}^{- 1} V_{A} {\hat{D}}_{ULC}^{- 1}) \\ ({\hat{M}}_{B U}^{r} - {\hat{D}}_{U R} (\begin{matrix} Φ_{f} & O \end{matrix})) \end{matrix}

(40)

Notice that ${\hat{D}}_{ULC}, {\hat{D}}_{U R}$ ⁠, and $Φ_{f}$ only depend on the spindle speed $Ω$ ⁠, meaning that the matrices $M_{P L} {\hat{D}}_{ULC}^{- 1}$ ⁠, $M_{P L} {\hat{D}}_{ULC}^{- 1} U_{A}$ ⁠, $V_{A} {\hat{D}}_{ULC}^{- 1}$ ⁠, $V_{A} {\hat{D}}_{ULC}^{- 1} U_{A}$ ⁠, and ${\hat{D}}_{U R} (\begin{matrix} Φ_{f} & O \end{matrix})$ 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

for

j = 1

to

N_{Ω}

do

Ω = Ω_{j}

Calculate

M_{P L} {\hat{D}}_{ULC}^{- 1}

⁠,

M_{P L} {\hat{D}}_{ULC}^{- 1} U_{A}

⁠,

V_{A} {\hat{D}}_{ULC}^{- 1} U_{A}

⁠,

V_{A} {\hat{D}}_{ULC}^{- 1}

⁠, and

{\hat{D}}_{U R} (\begin{matrix} Φ_{f} & O \end{matrix})

for

i = 1

to

N_{a}

do

a = a_{i}

Calculate

\begin{matrix} Φ_{r} = (M_{P L} {\hat{D}}_{ULC}^{- 1} + M_{P L} {\hat{D}}_{ULC}^{- 1} U_{A} \hat{H} \\ {(I_{2 N} - V_{A} {\hat{D}}_{ULC}^{- 1} U_{A} \hat{H})}^{- 1} V_{A} {\hat{D}}_{ULC}^{- 1}) ({\hat{M}}_{B U}^{r} - {\hat{D}}_{U R} (\begin{matrix} Φ_{f} & O \end{matrix})) \end{matrix}

Calculate the spectral radius

ρ (Φ_{r})

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 $l N \times l N$ matrix ${\hat{D}}_{ULC}$ is performed once; while inside the loop, a decomposition of the $2 N \times 2 N$ matrix is performed for the construction of the transition matrix $Φ_{r}$ of dimension $l + 2 (N - 1)$ then following the calculation of the spectral radius $ρ (Φ_{r})$ for every a. We denote the time for the decomposition of ${\hat{D}}_{ULC}$ as $t_{inv, FCCM}$ and the time spent on the spectral radius as $t_{eig, FCCM}$ ⁠. The comparison of $t_{eig, FCCM}$ and $t_{eig, FDM}$ is shown in Fig. 6.

Fig. 6

View large Download slide

$t_{eig, FCCM}$ and $t_{eig, FDM}$ ⁠: (a) value in parameter space $(l, N)$ and (b) $N = 120$

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

{\hat{D}}_{ULC}

outside the loop, it is not valid anymore to use the

t_{eig, FCCM}

and

t_{eig, 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

A_{0}

in

{\hat{D}}_{ULC}

⁠. Notice that

\begin{matrix} \begin{matrix} A_{0} = (\begin{matrix} - \frac{1}{2} {\hat{M}}^{- 1} \hat{C} & {\hat{M}}^{- 1} \\ \frac{1}{4} \hat{C} {\hat{M}}^{- 1} \hat{C} - \hat{K} & - \frac{1}{2} \hat{C} {\hat{M}}^{- 1} \end{matrix}) \end{matrix} \\ = \underset{U_{M}}{\underset{︸}{(\begin{matrix} I_{l / 2} \\ - \frac{1}{2} \hat{C} \end{matrix})}} {\hat{M}}^{- 1} \underset{V_{M}}{\underset{︸}{(\begin{matrix} - \frac{1}{2} \hat{C} & I_{l / 2} \end{matrix})}} - \underset{U_{K}}{\underset{︸}{(\begin{matrix} O \\ I_{l / 2} \end{matrix})}} \hat{K} \underset{V_{K}}{\underset{︸}{(\begin{matrix} I_{l / 2} & O \end{matrix})}} \\ = U_{M} {\hat{M}}^{- 1} V_{M} + U_{K} \hat{K} V_{K} \end{matrix}

(41)

Again, by allying the Sherman–Morrison–Woodbury formula, the inverse matrix

{\hat{D}}_{ULC}^{- 1}

has the following expression:

\begin{matrix} \begin{matrix} {\hat{D}}_{ULC}^{- 1} = {\hat{D}}_{ULC, K}^{- 1} + {\hat{D}}_{ULC, K}^{- 1} U_{M A} {(I_{N} \otimes \hat{M} + {\bar{D}}_{U L}^{- 1} \otimes \hat{C} + {\bar{D}}_{U L}^{- 2} \otimes \hat{K})}_{\frac{l N}{2} \times \frac{l N}{2}}^{- 1} V_{M A} {\hat{D}}_{ULC, K}^{- 1} \end{matrix} \\ \begin{matrix} {\hat{D}}_{ULC, K}^{- 1} = {\bar{D}}_{U L}^{- 1} \otimes I_{l} - {\bar{D}}_{U L}^{- 2} \otimes {(\begin{matrix} O & O \\ \hat{K} & O \end{matrix})}_{l \times l}, U_{M A} = I_{N} \otimes U_{M} \end{matrix} \end{matrix}

(42)

the derivation in detail is provided in Appendix C. Now, the inversion of ${\hat{D}}_{ULC}$ follows two steps: first, calculate ${\hat{D}}_{ULC, K}^{- 1}$ ⁠; next, calculate ${\hat{D}}_{ULC}^{- 1}$ ⁠. The main source of the computation cost changes from the decomposition of ${\hat{D}}_{ULC}$ to that of $(I_{N} \otimes \hat{M} + {\bar{D}}_{U L}^{- 1} \otimes \hat{C} + {\bar{D}}_{U L}^{- 2} \otimes \hat{K})$ ⁠, which is only half the size of ${\hat{D}}_{ULC}$ ⁠.

Define the time for decomposing the matrix ${\hat{D}}_{ULC}$ using Eq. (42) as $t_{inv, FCCM, r}$ ⁠, the comparison of $t_{inv, FCCM}$ and $t_{inv, FCCM, r}$ is shown in Fig. 7. Because the dimension of the matrix needed to be decomposed is only half of ${\hat{D}}_{ULC}$ ⁠, the decomposition time should be approximately a quarter of that of ${\hat{D}}_{ULC}$ ⁠, i.e., $t_{inv, FCCM, r} \approx \frac{1}{4} t_{inv, 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 $(I_{N} \otimes \hat{M} + {\bar{D}}_{U L}^{- 1} \otimes \hat{C} + {\bar{D}}_{U L}^{- 2} \otimes \hat{K})$ (denoted as $\hat{D_{ULC, r}}$ in the remaining of this article) compared to ${\hat{D}}_{ULC}$ ⁠. The distribution of the elements of ${\hat{D}}_{ULC}$ and ${\hat{D}}_{ULC, r}$ is shown in Fig. 8. In ${\hat{D}}_{ULC}$ ⁠, most blocks are identity matrices, however, along the main diagonal are the dense matrix $A_{0}$ ⁠, where almost all the elements are no-zero. On contrast, all the blocks of ${\hat{D}}_{ULC, r}$ are the linear combination of the matrices $\hat{M}$ ⁠, $\hat{C}$ ⁠, and $\hat{K}$ ⁠, which retain the sparse band structure due to discretization using the finite element analysis. The density of the original $A_{0}$ is in a sense “equally distributed” to the other blocks in ${\hat{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

View large Download slide

$t_{inv, FCCM}$ and $t_{inv, FCCM, r}$ ⁠: (a) value in parameter space $(l, N)$ and (b) $N = 120$

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

View large Download slide

Distribution of the nonzero elements in the sparse matrix ${\hat{D}}_{ULC}$ and ${\hat{D}}_{ULC, r}$ ⁠: (a) ${\hat{D}}_{ULC}$ and (b) ${\hat{D}}_{ULC, r}$

After this dimension reduction, in Algorithm 3, we use Eq. (42) to calculate ${\hat{D}}_{ULC}^{- 1}$ ⁠, 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) = (5000 rpm, 2 mm)$ 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

View large Download slide

Comparison of the convergence rate between FDM and FCCM: (a) $l = 24$ ⁠, (b) $l = 84$ ⁠, (c) $l = 132$ ⁠, and (d) $l = 192$

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 $10^{- 14}$ ⁠, 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 $10^{- 7}$ 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

t_{eig, FCCM}

and

t_{eig, FDM}

⁠, we define the new value, called

t_{loop, FCCM}

and

t_{loop, FDM}

to implement the comparison of computation efficiency, which means the time spent on one fixed spindle speed

Ω

when constructing the SLD

t_{loop, FCCM, r} = t_{inv, FCCM, r} + N_{a} * t_{eig, FCCM}, t_{loop, FDM} = N_{a} * t_{eig, FDM}

(43)

N_{a}

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

N_{a} = 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,

t_{loop, FCCM}

exceeds

t_{loop, FDM}

⁠, meaning that FCCM is slower. The reason can possibly be explained as follows. From Eq. ,

N_{a}

is also an essential parameter that decides the relative size of

t_{loop, FCCM}

and

t_{loop, FDM}

⁠. For a fixed

(l, N)

⁠, the relation between them and

N_{a}

is shown in Fig. 11. Denote the value of

N_{a}

when

t_{loop, FCCM}

and

t_{loop, FDM}

are equal as

N_{a, critical}

N_{a, critical} = [\frac{t_{inv, FCCM, r}}{t_{eig, FDM} - t_{eig, FCCM}}]

(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 $N_{a}$ exceeds $N_{a, critical}$ ⁠, meaning that FCCM is more computationally efficient only when the mesh of the SLD is sufficiently dense. The value of $N_{a, critical}$ under different $(l, N)$ is shown in Fig. 11. Fortunately, for most of the $(l, N)$ ⁠, $N_{a, critical}$ is below 100 except for an abrupt peak at the up-right corner of the $(l, N)$ space, where the maximum $N_{a, critical}$ reaches 382. By increasing $N_{a}$ ⁠, 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

View large Download slide

$t_{loop, FDM}$ and $t_{loop, FCCM}$ ⁠: (a) value in parameter space $(l, N)$ ⁠, (b) $N = 120$ ⁠, (c) the ratio $ρ_{t} = t_{loop, FDM} / t_{loop, FCCM}$ in parameter space $(l, N)$ ⁠, and (d) the distribution in the $(l, N)$ space whether FCCM is faster or not

Fig. 11

View large Download slide

$t_{loop, FCCM}$ and $t_{loop, FDM}$ as the function of $N_{a}$ ⁠: (a) the function and the critical point $N_{a, critical}$ and (b) $N_{a, critical}$ in space of $(l, N)$

Fig. 12

View large Download slide

Distribution in the $(l, N)$ space whether FCCM is faster or not with larger $N_{a}$ ⁠: (a) $N_{a} = 150$ and (b) $N_{a} = 200$

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:

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.
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.
In most cases considered in this article, the FCCM performs better than the original CCM and the FDM in computation efficiency. The values $t_{loop, FCCM}$ and $t_{loop, 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} = t_{loop, FDM} / t_{loop, FCCM}$ is below one except in several cases when l and N are both high (⁠ $l \geq 120$ ⁠, $N \geq 100$ ⁠). Increasing $N_{a}$ ⁠, 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

\begin{matrix} \begin{matrix} τ_{d} = M (q) q^{*} + H (q, \dot{q}) + G (q) \end{matrix} \\ \begin{matrix} q^{*} = J {(q)}^{- 1} ((\begin{matrix} {\dot{v}}_{d} \\ {\dot{ω}}_{d} \end{matrix}) + k_{d} (\begin{matrix} \begin{matrix} v_{d} \\ ω_{d} \end{matrix} \end{matrix}) - k_{p} (\begin{matrix} e_{p} \\ e_{o} \end{matrix}) - \dot{J} (q, \dot{q}) \dot{q}) - k_{d} \dot{q} \end{matrix} \end{matrix}

(A1)

where

\dot{J} (q, \dot{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 = {(n^{T}, s^{T}, a^{T})}^{T}

⁠, n, s, a all refer to

3 \times 1

vector, then

e_{p}

and

e_{o}

represent the error of position and orientation, respectively,

\begin{matrix} e_{p} = (\begin{matrix} x - x_{d} \\ y - y_{d} \\ z - z_{d} \end{matrix}), e_{o} = - \frac{1}{2} (n \times n_{d} + s \times s_{d} + a \times a_{d}) \end{matrix}

(A2)

according to Ref. [32], when the absolute value of the elements of $e_{o}$ is small, we approximately have ${\dot{e}}_{o} \approx ω - ω_{d}, {\ddot{e}}_{o} \approx \dot{ω} - {\dot{ω}}_{d}$ ⁠.

Since

(\begin{matrix} v \\ ω \end{matrix}) = J (q) \dot{q}

[31], deriving it on both sides

\begin{matrix} (\begin{matrix} \dot{v} \\ \dot{ω} \end{matrix}) = J (q) \ddot{q} + \dot{J} (q, \dot{q}) \dot{q}, \ddot{q} = J {(q)}^{- 1} ((\begin{matrix} \dot{v} \\ \dot{ω} \end{matrix}) - \dot{J} (q, \dot{q}) \dot{q}) \end{matrix}

(A3)

Substitute Eqs. and into Eq. , we have

\begin{matrix} M (q) (\ddot{q} - q^{*}) = 0 \\ M (q) J {(q)}^{- 1} ((\begin{matrix} \dot{v} - {\dot{v}}_{d} \\ \dot{ω} - {\dot{ω}}_{d} \end{matrix})) + k_{d} (J (q) \dot{q} - (\begin{matrix} v_{d} \\ ω_{d} \end{matrix})) + k_{p} (\begin{matrix} e_{p} \\ e_{o} \end{matrix})) = 0 \\ M (q) J {(q)}^{- 1} ((\begin{matrix} \dot{v} - {\dot{v}}_{d} \\ \dot{ω} - {\dot{ω}}_{d} \end{matrix}) + k_{d} (\begin{matrix} v - v_{d} \\ ω - ω_{d} \end{matrix}) + k_{p} (\begin{matrix} e_{p} \\ e_{o} \end{matrix})) = 0 \\ M (q) J {(q)}^{- 1} ((\begin{matrix} {\ddot{e}}_{p} \\ {\ddot{e}}_{o} \end{matrix}) + k_{d} (\begin{matrix} \begin{matrix} {\dot{e}}_{p} \\ {\dot{e}}_{o} \end{matrix} \end{matrix}) + k_{p} (\begin{matrix} \begin{matrix} e_{p} \\ e_{o} \end{matrix} \end{matrix})) = 0 \end{matrix}

(A4)

Then, Eq. is transformed into the Cartesian space of the end by left-multiplying

J {(q)}^{- T}

J {(q)}^{- T} M (q) J {(q)}^{- 1} ((\begin{matrix} {\ddot{e}}_{p} \\ {\ddot{e}}_{o} \end{matrix}) + k_{d} (\begin{matrix} {\dot{e}}_{p} \\ {\dot{e}}_{o} \end{matrix}) + k_{p} (\begin{matrix} e_{p} \\ e_{o} \end{matrix})) = 0

(A5)

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

Table 1

The modified D–H parameters of the robot

Number of joint	d (mm)	a (mm)	α (deg)
1	815	0	0
2	0	350	−90
3	0	850	0
4	820	145	−90
5	0	0	90
6	170	0	−90

Table 2

Dynamic parameters of the robot

Number of joint	Mass m (kg)	Center of mass (x_c, y_c, z_c) (mm)	Center inertia tensor (I_xx, I_yy, I_zz) (kg m²)
1	121.09	(140.890, 15.739, −123.659)	(2.577, 6.170, 6.725)
2	182.99	(401.195, −0.501, −248.760)	(17.503, 1.640, 17.070)
3	126.74	(60.468, −4.014, 13.270)	(5.235, 2.031, 4.968)
4	87.06	(−3.556, 0.12, −169.963)	(0.430, 0.430, 0.487)
5	39.72	(0.059, 32.403, 0.013)	(0.183, 0.120, 0.174)
6	44.25	(−3.177, 0.801, 113.311)	(0.247, 0.175, 0.136)

Number of joint	Mass m (kg)	Center of mass (x_c, y_c, z_c) (mm)	Center inertia tensor (I_xx, I_yy, I_zz) (kg m²)
1	121.09	(140.890, 15.739, −123.659)	(2.577, 6.170, 6.725)
2	182.99	(401.195, −0.501, −248.760)	(17.503, 1.640, 17.070)
3	126.74	(60.468, −4.014, 13.270)	(5.235, 2.031, 4.968)
4	87.06	(−3.556, 0.12, −169.963)	(0.430, 0.430, 0.487)
5	39.72	(0.059, 32.403, 0.013)	(0.183, 0.120, 0.174)
6	44.25	(−3.177, 0.801, 113.311)	(0.247, 0.175, 0.136)

Table 3

Other parameters of the robot

Quantities	Symbol	Value	Units
Coordinate of the end of the spindle in the frame of the last link of the robot	p_end	(−200, 0, 200)	mm
The direction cosine matrix of the end frame relative to the frame of the last link	A_end	$(\begin{matrix} 0 & 0 & 1 \\ 0 & - 1 & 0 \\ 1 & 0 & 0 \end{matrix})$	—
The RMAC control parameter	k_p	10⁸	s⁻²
	k_d	$2 \sqrt{k_{p}}$	s⁻¹

Quantities	Symbol	Value	Units
Coordinate of the end of the spindle in the frame of the last link of the robot	p_end	(−200, 0, 200)	mm
The direction cosine matrix of the end frame relative to the frame of the last link	A_end	$(\begin{matrix} 0 & 0 & 1 \\ 0 & - 1 & 0 \\ 1 & 0 & 0 \end{matrix})$	—
The RMAC control parameter	k_p	10⁸	s⁻²
	k_d	$2 \sqrt{k_{p}}$	s⁻¹

Table 4

Parameters of the cutting tool and the milling force (the parameters in Tables 1–3 are coming from Ref. [35], parameters $K_{t}, K_{r}$ in Table 4 are coming from Ref. [6], and the other parameters in Table 4 are coming from Ref. [36])

Quantities	Symbol	Value	Units
Number of the cutting teeth	N_tool	4	—
Diameter of the cutting tool	D_tool	25	mm
Length of the cutting tool	L	121	mm
Young's modulus of elasticity	E	231.932	GPa
Poisson ration	μ	0.324	—
Shear correction factor	k	$\frac{6 (1 + μ)}{7 + 6 μ}$	—
Density of the tool	ρ	8563.55	kg/m³
The mass coefficient of the Rayleigh damping	α_M	35.372	—
The stiffness coefficient of the Rayleigh damping	α_K	2.061 × 10⁻¹⁰	—
The direction cosine matrix of the local frame of the element relative to the end frame	T	$(\begin{matrix} 0 & 0 & - 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{matrix})$	—
Tangential coefficient of the cutting force	K_t	1835	MPa
Normal coefficient of the cutting force	K_r	1137	MPa

Quantities	Symbol	Value	Units
Number of the cutting teeth	N_tool	4	—
Diameter of the cutting tool	D_tool	25	mm
Length of the cutting tool	L	121	mm
Young's modulus of elasticity	E	231.932	GPa
Poisson ration	μ	0.324	—
Shear correction factor	k	$\frac{6 (1 + μ)}{7 + 6 μ}$	—
Density of the tool	ρ	8563.55	kg/m³
The mass coefficient of the Rayleigh damping	α_M	35.372	—
The stiffness coefficient of the Rayleigh damping	α_K	2.061 × 10⁻¹⁰	—
The direction cosine matrix of the local frame of the element relative to the end frame	T	$(\begin{matrix} 0 & 0 & - 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{matrix})$	—
Tangential coefficient of the cutting force	K_t	1835	MPa
Normal coefficient of the cutting force	K_r	1137	MPa

Appendix C: Derivation of Eq. (42)

According to Eqs. and , we have

\begin{matrix} \begin{matrix} {\hat{D}}_{ULC} = {\bar{D}}_{U L} \otimes I_{l} - I_{N} \otimes A_{0} \end{matrix} \\ = \underset{{\bar{D}}_{U L, l}}{\underset{︸}{{\bar{D}}_{U L} \otimes I_{l}}} + \underset{U_{K A}}{\underset{︸}{(I_{N} \otimes U_{K})}} (I_{N} \otimes \hat{K}) \underset{V_{K A}}{\underset{︸}{(I_{N} \otimes V_{K})}} \\ - \underset{U_{M A}}{\underset{︸}{(I_{N} \otimes U_{M})}} (I_{N} \otimes {\hat{M}}^{- 1}) \underset{V_{M A}}{\underset{︸}{(I_{N} \otimes V_{M})}} \\ = \underset{{\hat{D}}_{ULC, K}}{\underset{︸}{{\bar{D}}_{U L, l} + U_{K A} (I_{N} \otimes \hat{K}) V_{K A}}} - U_{M A} (I_{N} \otimes {\hat{M}}^{- 1}) V_{M A} \\ = {\hat{D}}_{ULC, K} - U_{M A} (I_{N} \otimes {\hat{M}}^{- 1}) V_{M A} \end{matrix}

(C1)

By applying the Sherman–Morrison–Woodbury formula, the associated property of the Kronecker product,

(A \otimes B) (C \otimes D) = (A C) \otimes (B D)

⁠, and the inversion of the Kronecker product,

{(A \otimes B)}^{- 1} = A^{- 1} \otimes B^{- 1}

\begin{matrix} {\hat{D}}_{ULC, K}^{- 1} = {\bar{D}}_{U L, l}^{- 1} - {\bar{D}}_{U L, l}^{- 1} U_{K A} {(I_{N} \otimes {\hat{K}}^{- 1} + V_{K A} {\bar{D}}_{U L, l}^{- 1} U_{K A})}^{- 1} V_{K A} {\bar{D}}_{U L, l}^{- 1} \\ = {\bar{D}}_{U L}^{- 1} \otimes I_{l} - {\bar{D}}_{U L, l}^{- 1} U_{K A} (I_{N} \otimes K^{- 1} \\ + D_{N}^{- 1} \otimes ((\begin{matrix} I & O \end{matrix}) (\begin{matrix} O \\ I \end{matrix})))^{- 1} V_{K A} {\bar{D}}_{U L, l}^{- 1}] \\ = {\bar{D}}_{U L}^{- 1} \otimes I_{l} - {\bar{D}}_{U L, l}^{- 1} U_{K A} (I_{N} \otimes \hat{K}) V_{K A} {\bar{D}}_{U L, l}^{- 1} \\ = {\bar{D}}_{U L}^{- 1} \otimes I_{l} - {\bar{D}}_{U L}^{- 2} \otimes ((\begin{matrix} O \\ I \end{matrix}) \hat{K} (\begin{matrix} I & O \end{matrix})) \\ = {\bar{D}}_{U L}^{- 1} \otimes I_{l} - {\bar{D}}_{U L}^{- 2} \otimes {(\begin{matrix} O & O \\ \hat{K} & O \end{matrix})}_{l \times l} \end{matrix}

(C2)

Then

\begin{matrix} \begin{matrix} {\hat{D}}_{ULC}^{- 1} = {\hat{D}}_{ULC, K}^{- 1} + {\hat{D}}_{ULC, K}^{- 1} U_{M A} (I_{N} \otimes M \end{matrix} \\ - V_{M A} {\hat{D}}_{ULC, K}^{- 1} U_{M A})^{- 1} V_{M A} {\hat{D}}_{ULC, K}^{- 1} \\ = {\hat{D}}_{ULC, K}^{- 1} + {\hat{D}}_{ULC, K}^{- 1} U_{M A} (I_{N} \otimes M - {\bar{D}}_{U L}^{- 1} \otimes (\begin{matrix} - \frac{1}{2} \hat{C} & I_{l / 2} \end{matrix}) (\begin{matrix} I_{l / 2} \\ - \frac{1}{2} \hat{C} \end{matrix}) \\ + {\bar{D}}_{U L}^{- 2} \otimes (\begin{matrix} - \frac{1}{2} \hat{C} & I_{l / 2} \end{matrix}) (\begin{matrix} O & O \\ \hat{K} & O \end{matrix}) (\begin{matrix} I_{l / 2} \\ - \frac{1}{2} \hat{C} \end{matrix}))^{- 1} V_{M A} {\hat{D}}_{ULC, K}^{- 1} \\ = {\hat{D}}_{ULC, K}^{- 1} + {\hat{D}}_{ULC, K}^{- 1} U_{M A} (I_{N} \otimes \hat{M} \\ + {\bar{D}}_{U L}^{- 1} \otimes \hat{C} + {\bar{D}}_{U L}^{- 2} \otimes \hat{K})_{\frac{l N}{2} \times \frac{l N}{2}}^{- 1} V_{M A} {\hat{D}}_{ULC, K}^{- 1} \end{matrix}

(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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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

Google Scholar

Crossref

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.

A Fast Chebyshev Collocation Method for Stability Analysis of a Robotic Machining System With Time Delay

Abstract

1 Introduction

2 Modeling of the Robotic Milling System

2.1 Structural Dynamics of the Robotic Milling System.

2.2 Cutting Force Model.

3 Chebyshev Collocation Method and the Evaluation of Computational Efficiency

3.1 Chebyshev Collocation Method.

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

4 The Accelerated Algorithm: Fast Chebyshev Collocation Method

4.1 Dimension Reduction of the Transition Matrix.

4.2 Acceleration of Matrix Inversion.

4.3 Convergence Analysis of Fast Chebyshev Collocation Method.

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

5 Concluding Remarks

Acknowledgment

Funding Data

Data Availability Statement

Appendix A: Derivation of Eq. (4)

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

Appendix C: Derivation of Eq. (42)

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

A Fast Chebyshev Collocation Method for Stability Analysis of a Robotic Machining System With Time Delay

Abstract

1 Introduction

2 Modeling of the Robotic Milling System

2.1 Structural Dynamics of the Robotic Milling System.

2.2 Cutting Force Model.

3 Chebyshev Collocation Method and the Evaluation of Computational Efficiency

3.1 Chebyshev Collocation Method.

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

4 The Accelerated Algorithm: Fast Chebyshev Collocation Method

4.1 Dimension Reduction of the Transition Matrix.

4.2 Acceleration of Matrix Inversion.

4.3 Convergence Analysis of Fast Chebyshev Collocation Method.

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

5 Concluding Remarks

Acknowledgment

Funding Data

Data Availability Statement

Appendix A: Derivation of Eq. (4)

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

Appendix C: Derivation of Eq. (42)

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only