*4.2. Numerical Algorithm*

The Lagrange polynomial denoted by *PL*(*x*) [16,24] can be used to approximate *f*(*x*) as follows:

$$P\_L(\mathbf{x}) = \sum\_{k=0}^{L} f(\mathbf{x}\_k) \eta\_{L,k}(\mathbf{x}) \tag{39}$$

where

$$\eta\_{L,k}(\mathbf{x}) = \prod\_{\substack{0 \le m \le k \\ m \ne L}} \frac{\mathbf{x} - \mathbf{x}\_m}{\mathbf{x}\_L - \mathbf{x}\_m} \tag{40}$$

It is worth noting that *x<sup>i</sup>* , *xi*+1, and *xi*+2 only exist when *L* = 2, and then, the Simpson's 3/8 formula (41) can be derived using *P*2(*x*).

$$\int\_{\mathbf{x}\_{i}}^{\mathbf{x}\_{i+2}} f(\mathbf{x})d\mathbf{x} = \frac{\hbar}{3}(f(\mathbf{x}\_{i}) + 4f(\mathbf{x}\_{i+1}) + f(\mathbf{x}\_{i+2})) \tag{41}$$

Based on (41), (38) can be rewritten as follows:

$$\begin{aligned} \mathbf{x}(t\_3) &= \begin{aligned} \varepsilon^{2\mathbf{A}h} \mathbf{x}(t\_1) + \frac{a\_2 h}{3} \Big\{ \varepsilon^{2\mathbf{A}h} \mathbf{B}(t\_1) [\mathbf{x}(t\_1) - \mathbf{x}(t\_1 - T)] + 4 \mathbf{e}^{\mathbf{A}h} \mathbf{B}(t\_2) [\mathbf{x}(t\_2) - \mathbf{x}(t\_2 - T)] \\ &+ \mathbf{B}(t\_3) [\mathbf{x}(t\_3) - \mathbf{x}(t\_3 - T)] \end{aligned} \tag{42}$$

which can be expressed as follows:

$$\begin{array}{lcl} \mathbf{x}(t\_{i+2}) &= \mathbf{e}^{2\mathbf{A}h} \mathbf{x}(t\_i) + \frac{a\_l h}{3} \Big\{ \mathbf{e}^{2\mathbf{A}h} \mathbf{B}(t\_i) [\mathbf{x}(t\_i) - \mathbf{x}(t\_i - T)] + 4 \mathbf{e}^{\mathbf{A}h} \mathbf{B}(t\_{i+1}) [\mathbf{x}(t\_{i+1}) - \mathbf{x}(t\_{i+1} - T)] \\ &+ \mathbf{B}(t\_{i+2}) [\mathbf{x}(t\_{i+2}) - \mathbf{x}(t\_{i+2} - T)] \Big\} \end{array} \tag{43}$$

As **B**(*ti+*2) cannot be solved directly when *i* = *n*, the classical numerical integration method is introduced to determine **x**(*tn*+1) as follows:

$$\mathbf{x}(t\_{n+1}) = e^{\mathbf{A}t}\mathbf{x}(t\_n) + \frac{a\_p \hbar}{2} \left\{ e^{\mathbf{A}t} \mathbf{B}(t\_n) [\mathbf{x}(t\_n) - \mathbf{x}(t\_n - T)] + \mathbf{B}(t\_{n+1}) [\mathbf{x}(t\_{n+1}) - \mathbf{x}(t\_{n+1} - T)] \right\} \tag{44}$$

Combining (17), (43), and (44), the transition matrix can be constructed as follows:

$$\mathbf{C}\_{2} = \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ -e^{2\mathbf{A}h} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -e^{2\mathbf{A}h} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & -e^{2\mathbf{A}h} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & -e^{\mathbf{A}h} & \mathbf{0} \end{bmatrix} \tag{45}$$

$$\mathbf{D}\_{2} = \begin{bmatrix} \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ -\varepsilon^{2\text{Ah}} \mathbf{B}\_{1} & -4\varepsilon^{\text{Ah}} \mathbf{B}\_{2} & -\mathbf{B}\_{3} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -\varepsilon^{2\text{Ah}} \mathbf{B}\_{2} & -4\varepsilon^{\text{Ah}} \mathbf{B}\_{3} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & -4\varepsilon^{\text{Ah}} \mathbf{B}\_{i-1} & -\mathbf{B}\_{i} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & -\varepsilon^{2\text{Ah}} \mathbf{B}\_{i-1} & -4\varepsilon^{\text{Ah}} \mathbf{B}\_{i} & -\mathbf{B}\_{i+1} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \frac{3}{2}\varepsilon^{\text{Ah}} \mathbf{B}\_{i} & \frac{3}{2}\mathbf{B}\_{i+1} \end{bmatrix} \tag{46}$$

$$\mathbf{E}\_2 = \begin{bmatrix} \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & e^{\mathbf{A}t\_f} \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \end{bmatrix} \tag{47}$$

$$\mathbf{I}\_2 = \begin{bmatrix} \mathbf{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{1} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \\ \end{bmatrix} \tag{48}$$

According to (45)–(48), the dynamic mapping of the discrete system can be determined as follows:

$$\mathbf{x}\left(\mathbf{I}\_2 + \mathbf{C}\_2 + \frac{a\_p h}{3}\mathbf{D}\_2\right) \begin{bmatrix} \mathbf{x}(t\_1) \\ \mathbf{x}(t\_2) \\ \vdots \\ \mathbf{x}(t\_{n+1}) \end{bmatrix} = \begin{pmatrix} a\_p h \\ \frac{a\_p h}{3}\mathbf{D}\_2 + \mathbf{E}\_2 \end{pmatrix} \begin{bmatrix} \mathbf{x}(t\_1 - T) \\ \mathbf{x}(t\_2 - T) \\ \vdots \\ \mathbf{x}(t\_{n+1} - T) \end{bmatrix} \tag{49}$$

Therefore, the transition matrix of the dynamic system in a tool pass cycle is obtained as follows:

$$\boldsymbol{\Phi}\_{2} = \left(\mathbf{I}\_{2} + \mathbf{C}\_{2} + \frac{a\_{p}h}{3}\mathbf{D}\_{2}\right)^{-1} \left(\frac{a\_{p}h}{3}\mathbf{D}\_{2} + \mathbf{E}\_{2}\right) \tag{50}$$

### *4.3. Simulation and Analysis*

4.3.1. Single-DOF Milling Model

A benchmark example of the single-DOF milling model is utilized to validate and analyze S38M. Its state-space expression is introduced as (31), where

$$\mathbf{A} = \begin{bmatrix} 0 & 1 \\ -\omega\_n^2 & -2\zeta\omega\_n \end{bmatrix} \tag{51}$$

$$\mathbf{B}(t) = \begin{bmatrix} 0 & 0 \\ -\frac{\hbar(t)}{m\_l} & 0 \end{bmatrix} \tag{52}$$

To analyze the convergence rate of S38M, the radial immersion ratio *a*/*D* is set as 1 to avoid intermittent milling. The depths of cut are set as 0.0003 m, 0.0008 m, and 0.0015 m. Figure 7 shows the convergence rate comparisons of the S38M, CCM, SDM, and SEM when the spindle speed is 8000 rpm. It is worth noting that the S38M exhibits the highest convergence rate compared with other chatter analysis methods. According to the enlarged figures on the right, the S38M has a good convergence rate even when the discrete number is large. Additionally, the results show that the convergence rate of the SEM fluctuates evidently, and the SEM is not as stable as the other methods. To further analyze the convergence rate of the S38M, the spindle speed is raised to 10,000 rpm, and the other parameters remain unchanged. Figure 8 shows the resultant convergence rates. It can be seen that these results are consistent with the previous results.

The low-immersion condition can be utilized to verify the stability of the convergence rate [25]. The radial immersion ratio is set as 0.05; that is, *a*/*D* = 0.05. The depths of cut are set as 0.0001 m, 0.0002 m, and 0.0003 m. The spindle speed was set to 5000 rpm. Figure 9 shows the convergence rate comparisons of S38M, CCM, SDM, and SEM under low immersion. The convergence rates arranged from high to low are in the following order: S38M, CCM, SDM, and SEM. Moreover, it is seen that the convergence rate of SEM fluctuated remarkably.

**Figure 7.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when *v* = 8000 rpm and *a*/*D* = 1. (**a**) *ap* = 0.0003 m and *n*∈[30,100]; (**b**) *ap* = 0.0003 m and *n*∈[75,100]; (**c**) *ap* = 0.0008 m and *n*∈ [30,100]; (**d**) *ap* = 0.0008 m and *n*∈[75,100]; (**e**) *ap* = 0.0015 m and *n*∈[30,100]; (**f**) *ap* = 0.0015 m and *n* ∈[75,100]. **Figure 7.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when *v* = 8000 rpm and *a*/*D* = 1. (**a**) *a<sup>p</sup>* = 0.0003 m and *n* ∈ [30, 100]; (**b**) *a<sup>p</sup>* = 0.0003 m and *n* ∈ [75, 100]; (**c**) *a<sup>p</sup>* = 0.0008 m and *n* ∈ [30, 100]; (**d**) *a<sup>p</sup>* = 0.0008 m and *n* ∈ [75, 100]; (**e**) *a<sup>p</sup>* = 0.0015 m and *n* ∈ [30, 100]; (**f**) *a<sup>p</sup>* = 0.0015 m and *n* ∈ [75, 100].

**Figure 8.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when *v* = 10000 rpm and *a*/*D* = 1. (**a**) *ap* = 0.0003 m and *n*∈[30,100]; (**b**) *ap* = 0.0003 m and *n*∈[75,100]; (**c**) *ap* = 0.0008 m and *n*∈[30,100]; (**d**) *ap* = 0.0008 m and *n*∈[75,100]; (**e**) *ap* = 0.0015 m and *n*∈[30,100]; (**f**) *ap* = 0.0015 m and *n*∈[75,100]. **Figure 8.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when *v* = 10,000 rpm and *a*/*D* = 1. (**a**) *a<sup>p</sup>* = 0.0003 m and *n* ∈ [30, 100]; (**b**) *a<sup>p</sup>* = 0.0003 m and *n* ∈ [75, 100]; (**c**) *a<sup>p</sup>* = 0.0008 m and *n* ∈ [30, 100]; (**d**) *a<sup>p</sup>* = 0.0008 m and *n* ∈ [75, 100]; (**e**) *a<sup>p</sup>* = 0.0015 m and *n* ∈ [30, 100]; (**f**) *a<sup>p</sup>* = 0.0015 m and *n* ∈ [75, 100].

The low-immersion condition can be utilized to verify the stability of the convergence rate [25]. The radial immersion ratio is set as 0.05; that is, *a*/*D* = 0.05. The depths of cut are tuated remarkably.

**Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m.

set as 0.0001 m, 0.0002 m, and 0.0003 m. The spindle speed was set to 5000 rpm. Figure 9 shows the convergence rate comparisons of S38M, CCM, SDM, and SEM under low immersion. The convergence rates arranged from high to low are in the following order: S38M, CCM, SDM, and SEM. Moreover, it is seen that the convergence rate of SEM fluc-

The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as

accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The **Figure 9.** Convergence rate comparisons of the S38M, CCM, SDM, and SEM when and *a*/*D* = 0.05. (**a**) *ap* = 0.0001 m; (**b**) *ap* = 0.0002 m; (**c**) *ap* = 0.0003 m. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when *a*/*D* = 1. However, as *a*/*D* decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The *AMRE* and *MSE* are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when *a*/*D* = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when *a*/*D* = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The

As described above, the fluctuations appear in Figure 9. It is found that for the S38M, the convergence rate is often a local maximum when the discrete number is an odd one; the convergence rate is often a local minimum when the discrete number is an even one. To further study this problem, the discrete numbers are set as 30, 40, 50, 60, 70, and 80. The corresponding *AMRE* and *MSE* are plotted in Figure 12. The results show that the S38M attains a better computational accuracy when the discrete number is an even one. 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large. **Table 4.** Stability lobe diagrams of the S38M.

**Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. **Table 4.** Stability lobe diagrams of the S38M. *n* **25 40 55**  *n* **25 40 55**  *n* **25 40 55** 

**Table 4.** Stability lobe diagrams of the S38M.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 21 of 25

**Figure 10.** *AMRE* and *MSE* of the S38M, CCM, SDM, and SEM. (**a**) *AMRE* when *a*/*D* = 0.05; (**b**) *AMRE* when *a*/*D* = 1; (**c**) *MSE* when *a*/*D* = 0.05; (**d**) *MSE* when *a*/*D* = 1. **Figure 10.** *AMRE* and *MSE* of the S38M, CCM, SDM, and SEM. (**a**) *AMRE* when *a*/*D* = 0.05; (**b**) *AMRE* when *a*/*D* = 1; (**c**) *MSE* when *a*/*D* = 0.05; (**d**) *MSE* when *a*/*D* = 1. **Figure 10.** *AMRE* and *MSE* of the S38M, CCM, SDM, and SEM. (**a**) *AMRE* when *a*/*D* = 0.05; (**b**) *AMRE* when *a*/*D* = 1; (**c**) *MSE* when *a*/*D* = 0.05; (**d**) *MSE* when *a*/*D* = 1. **Figure 10.** *AMRE* and *MSE* of the S38M, CCM, SDM, and SEM. (**a**) *AMRE* when *a*/*D* = 0.05; (**b**) *AMRE* when *a*/*D* = 1; (**c**) *MSE* when *a*/*D* = 0.05; (**d**) *MSE* when *a*/*D* = 1.

**Figure 11.** Computational time of the S38M, CCM, SDM, and SEM. (**a**) *a*/*D* = 0.05; (**b**) *a*/*D* = 1. **Figure 11.** Computational time of the S38M, CCM, SDM, and SEM. (**a**) *a*/*D* = 0.05; (**b**) *a*/*D* = 1. **Figure 11. Figure 11.**  Computational time of the S38M, CCM, SDM, and SEM. ( Computational time of the S38M, CCM, SDM, and SEM. (**a**) **a***a*/*D* = 0.05; (**b**) *a*/*D* = 1. ) *a*/*D* = 0.05; (**b**) *a*/*D* = 1.

( )

**Figure 12.** *AMRE* and *MSE* of the S38M for the fluctuation analysis when *a*/*D* = 0.05. (**a**) *AMRE*; (**b**) *MSE.* **Figure 12.** *AMRE* and *MSE* of the S38M for the fluctuation analysis when *a*/*D* = 0.05. (**a**) *AMRE*; (**b**) *MSE*.

0 20 0

*m m x t xt xt m*

( )

2

ω

( )

As described above, the fluctuations appear in Figure 9. It is found that for the S38M, the convergence rate is often a local maximum when the discrete number is an odd one; the convergence rate is often a local minimum when the discrete number is an even one. To further study this problem, the discrete numbers are set as 30, 40, 50, 60, 70, and 80. The corresponding *AMRE* and *MSE* are plotted in Figure 12. The results show that the S38M attains a better computational accuracy when the discrete number is an even one.

### 4.3.2. Two-DOF milling model 4.3.2. Two-DOF Milling Model

The two-DOF milling mathematical model is expressed as follows [17]: The two-DOF milling mathematical model is expressed as follows [17]:

*t tn t n*

ζω

$$
\begin{aligned}
\begin{bmatrix}m\_l & 0\\0 & m\_l\end{bmatrix} \begin{bmatrix}\ddot{\mathbf{x}}(t)\\ \ddot{y}(t)\end{bmatrix} + \begin{bmatrix}2m\_l\zeta\omega\_n & 0\\0 & 2m\_l\zeta\omega\_n\end{bmatrix} \begin{bmatrix}\dot{\mathbf{x}}(t)\\ \dot{y}(t)\end{bmatrix} + \begin{bmatrix}m\_l\omega\_n^2 & 0\\0 & m\_l\omega\_n^2\end{bmatrix} \begin{bmatrix}\mathbf{x}(t)\\ y(t)\end{bmatrix} \\
= -a\_p \begin{bmatrix}h\_{\text{xx}}(t) & h\_{\text{xy}}(t)\\ h\_{\text{yx}}(t) & h\_{\text{yy}}(t)\end{bmatrix} \times \left\{\begin{bmatrix}\mathbf{x}(t)\\ y(t)\end{bmatrix} - \begin{bmatrix}\mathbf{x}(t-T)\\ y(t-T)\end{bmatrix}\right\}\end{aligned} \tag{53}
$$

where the tool is symmetrical by default. *mt*, *ζ*, and *ωn* in the *x*-direction are identical to those in the *y*-direction. *hxx*(*t*), *hxy*(*t*), *hyx*(*t*), and *hyy*(*t*) are the same as those defined in (3)– (6). where the tool is symmetrical by default. *m<sup>t</sup>* , *ζ*, and *ω<sup>n</sup>* in the *x*-direction are identical to those in the *y*-direction. *hxx*(*t*), *hxy*(*t*), *hyx*(*t*), and *hyy*(*t*) are the same as those defined in (3)–(6).

A new state vector is defined as () () () () () ,,, *<sup>T</sup> t xt* <sup>=</sup> *<sup>y</sup> t xt <sup>y</sup> <sup>t</sup>* **<sup>r</sup>** . Then, (53) can be expressed as follows: A new state vector is defined as . **r**(*t*) = - *x*(*t*), *y*(*t*), . *x*(*t*), . *y*(*t*) *T* . Then, (53) can be expressed as follows: .

$$\dot{\mathbf{r}}(t) = \mathbf{A}\mathbf{r}(t) + a\_p \mathbf{B}(t)[\mathbf{r}(t) - \mathbf{r}(t-T)] \tag{54}$$

Where where

$$\mathbf{A} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\omega\_n^2 & 0 & -2\zeta\omega\_n & 0 \\ 0 & -\omega\_n^2 & 0 & -2\zeta\omega\_n \end{bmatrix} \tag{55}$$

and

and

$$\mathbf{B} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -h\_{\rm xx}(t)/m\_t & -h\_{\rm xy}(t)/m\_t & 0 & 0 \\ -h\_{\rm yx}(t)/m\_t & -h\_{\rm yy}(t)/m\_t & 0 & 0 \end{bmatrix} \tag{56}$$

The following deduction is the same as that in Section 4.1.

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

2

*n n*

**Table 5.** Comparisons of S38M and SDM for the two-DOF milling model. **Table 5.** Comparisons of S38M and SDM for the two-DOF milling model. **Table 5.** Comparisons of S38M and SDM for the two-DOF milling model. **Table 5.** Comparisons of S38M and SDM for the two-DOF milling model. **Table 5.** Comparisons of S38M and SDM for the two-DOF milling model.

SDM. Additionally, the computation time is remarkably saved using the S38M.

SDM. Additionally, the computation time is remarkably saved using the S38M.

SDM. Additionally, the computation time is remarkably saved using the S38M.

SDM. Additionally, the computation time is remarkably saved using the S38M.

SDM. Additionally, the computation time is remarkably saved using the S38M.

SDM. Additionally, the computation time is remarkably saved using the S38M.

() () () ()

() () () ()

<sup>=</sup> − −

<sup>=</sup> − −

<sup>=</sup> − −

<sup>=</sup> − −

<sup>=</sup> − −

<sup>=</sup> − −

The following deduction is the same as that in Section 4.1.

The following deduction is the same as that in Section 4.1.

The following deduction is the same as that in Section 4.1.

The following deduction is the same as that in Section 4.1.

The following deduction is the same as that in Section 4.1.

The following deduction is the same as that in Section 4.1.

*htm htm htm htm*

*htm htm htm htm*

() () () ()

*htm htm htm htm*

*xx t xy t yx t yy t*

*xx t xy t yx t yy t*

*xx t xy t yx t yy t*

0 0 00 0 0 00 / / 00 / / 00

0 0 00 0 0 00 / / 00 / / 00

0 0 00 0 0 00 / / 00 / / 00

 

 

 

− −

− −

− −

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of

0 0 00 0 0 00 / / 00 / / 00

0 0 00 0 0 00 / / 00 / / 00

0 0 00 0 0 00 / / 00 / / 00

**B** (56)

**B** (56)

**B** (56)

**B** (56)

**B** (56)

**B** (56)

 

 

 

− −

− −

− −

() () () ()

() () () ()

() () () ()

*htm htm htm htm*

*htm htm htm htm*

*htm htm htm htm*

*xx t xy t yx t yy t*

*xx t xy t yx t yy t*

*xx t xy t yx t yy t*

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of

The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 23 of 25

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 23 of 25

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 23 of 25

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 23 of 25

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 23 of 25

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 23 of 25

#### **5. Conclusions 5. Conclusions 5. Conclusions 5. Conclusions**

**5. Conclusions** 

**5. Conclusions** 

**5. Conclusions** 

In this study, we focused on the prediction of milling stability using composite Cotesbased and Simpson's 3/8-based methods. First, the composite Cotes-based method is proposed for preventing chatter in milling processes. The three steps for obtaining the milling stability lobe diagrams are as follows: (1) establish the time-delay differential equation while considering the regenerative effects; (2) determine the transition matrix using the integral equations; (3) analyze the system stability according to the Floquet theory. To measure the calculation accuracy, the *AMRE* and *MSE* are adopted in our work. The results demonstrate that for the single-DOF model, when the discrete number is 40, the *AMRE* and *MSE* can be respectively reduced from 0.076 to 0.041 and from 4.24 <sup>×</sup> <sup>10</sup>−<sup>8</sup> to 2.62 <sup>×</sup> <sup>10</sup>−<sup>8</sup> , and the calculation time can be reduced from 55.63s to 20.21s by using the proposed composite Cotes-based method. In addition, the Simpson's 3/8-based method is proposed to further improve the calculation efficiency. The results demonstrate that for the single-DOF model, the proposed Simpson's 3/8-based method significantly improves the convergence speed while sharing the same computation time with the composite Cotes-based method when the radial immersion ratio is large; for the two-DOF model,

the accuracy of the stability lobe diagrams of the Simpson's 3/8-based method are better than those of the SDM, and the computation time is remarkably saved using the Simpson's 3/8-based method.

**Author Contributions:** Conceptualization, X.D. and J.Z.; methodology, X.D.; software, P.R.; validation, X.D., P.R. and J.Z.; formal analysis, P.R.; investigation, X.D.; resources, X.D.; data curation, P.R.; writing—original draft preparation, X.D.; writing—review and editing, J.Z.; visualization, X.D.; supervision, J.Z.; project administration, J.Z.; funding acquisition, J.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [National Natural Science Foundation of China] grant number [52005142] and [Open Fund of State Key Laboratory of Robotics and System] grant number [SKLRS-2021-KF-08]. The authors would like to thank the Mechatronic Institute, Zhejiang Sci-Tech University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.
