**1. Introduction**

The three main types of vibration that occur during high-speed milling are free vibration, excited vibration, and self-excited vibration, and the regenerative chatter in self-excited vibration is the main cause of instability in the machining process [1]. Because this chatter typically leads to poor surface quality of the machined parts, aggravated tool wear, and even reduced service life of machine tools, it must be avoided to solve or overcome these problems [2]. In the analysis of chatter, dynamic milling processes that consider the delay effect are generally described using delay differential equations with respect to the timeperiodic coefficients [3,4]. The chatter stability based on these differential equations is an important index for achieving high-performance machining [5].

In recent decades, experimental, analytical, and numerical methods have been studied to obtain stability lobe diagrams, which can be used to determine exact values avoiding chatter. Wu et al. [6] utilized the Lyapunov index to measure whether chatter occurred; however, they could only predict the system stability under specific parameters. Davies et al. [7] used the time-domain calculation method to predict unstable regions in the stability lobe diagrams, which was only applicable to small radial depths of cut. Altintas and Budak et al. [8] first presented a zero-order approximation method to quickly obtain the stability lobe diagrams, where the real and imaginary parts of the feature are used to determine the cutting parameters. It is worth mentioning that this method was only applicable to high radial depths of cut. Bayly et al. [9] employed the time-finite element method to calculate

**Citation:** Du, X.; Ren, P.; Zheng, J. Predicting Milling Stability Based on Composite Cotes-Based and Simpson's 3/8-Based Methods. *Micromachines* **2022**, *13*, 810. https://doi.org/10.3390/ mi13050810

Academic Editors: Youqiang Xing, Xiuqing Hao and Duanzhi Duan

Received: 28 April 2022 Accepted: 21 May 2022 Published: 23 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

tool motion and utilized the weighted residual method to determine the Floquet transition matrix. Based on the delay differential equations, Insperger et al. [10] proposed the semidiscretization method (SDM) and the first-order semi-discretization method (1stSDM), which are widely used as the basis for evaluating other stability analysis methods in the time domain. The SDM and 1stSDM only discretized delay states and periodic coefficients. Based on direct integration, Ding et al. [11] presented a full-discretization method (FDM) to predict milling stability. Unlike the SDM and 1stSDM, the FDM carried out the linear interpolation on the state terms to improve computational efficiency. Further, Insperger [12] compared the FDM with the SDM and 1stSDM in terms of the convergence speed and computational efficiency. It was found that the convergence speed of the 1stSDM was faster than those of the SDM and FDM, while the computational efficiency of the FDM was higher than those of the SDM and 1stSDM. After the comparisons, some FDM-based methods were used to predict the milling stability. In addition to the above methods, the numerical integration methods are widely utilized to obtain the stability lobe diagrams. Ding et al. [13] approximated the integral terms with the Newton–Cotes and Gauss integral formulas to predict the milling stability. Lu et al. [14] approximated the solution of the delay differential equation with the direct integration technique. Qin et al. [15] approximated the state terms of the delay differential equation with the Chebyshev wavelets of the second kind. Moreover, Liu et al. [16] combined the numerical solution of the delay differential equation with the Simpson formula to predict milling stability (SEM), which was only suitable for large radial immersion. However, the aforementioned numerical methods cannot simultaneously guarantee fast convergence speed and high computational efficiency.

To this end, we present two numerical analysis methods, namely the composite Cotesbased method (CCM) and Simpson's 3/8-based method (S38M), in this paper to improve the convergence speed and computational efficiency at the same time. The remainder of this paper is organized as follows. Section 2 presents the mathematical model of system motion with two degrees of freedom (DOF) and the state-space. The composite Cortesbased method and its simulation and analysis for the single-DOF model are presented in Section 3. Subsequently, the Simpson's 3/8-based method and its simulation and analysis for the single-DOF and two-DOF models are presented in Section 4. Finally, the conclusions are stated in Section 5.

### **2. Mathematical Model**

Taking down-milling as an example, the dynamic system of end-milling with two degrees of freedom is shown in Figure 1, where the workpiece is assumed to be rigid, and the milling cutter is assumed to be evenly distributed. Note that the helix angle of the milling cutter is not taken into account. Considering the regenerative effect, the mathematical model of system motion can be expressed as a second-order differential equation [17] as follows:

$$\mathbf{M}\ddot{\mathbf{q}}(t) + \mathbf{C}\dot{\mathbf{q}}(t) + \mathbf{K}\mathbf{q}(t) = -a\_p \mathbf{K}\_c(t)[\mathbf{q}(t) - \mathbf{q}(t - T)] \tag{1}$$

where **q**(*t*), . **<sup>q</sup>**(*t*), and .. **q**(*t*) represent the displacement vector, the first derivative of **q**(*t*) w.r.t *t*, and the second derivative of **q**(*t*) w.r.t *t*, respectively. **M**, **C**, **K**, and *a<sup>p</sup>* represent the mass matrix, damping matrix, stiffness matrix, and depth of cut, respectively. *T* = 60/(*Nv*) represents the delay period of the delay differential equations, where *N* represents the cutter tooth number, and *v* represents the spindle speed. The radial cutting force coefficient matrix **K***c*(*t*) can be expressed as follows:

$$\mathbf{K}\_c(t) = \begin{bmatrix} h\_{xx}(t) & h\_{xy}(t) \\ h\_{yx}(t) & h\_{yy}(t) \end{bmatrix} \tag{2}$$

where

$$h\_{\rm xx}(t) = \sum\_{j=1}^{N} g\left(\phi\_{\vec{j}}(t)\right) \sin\left(\phi\_{\vec{j}}(t)\right) \left[K\_{l}\cos\left(\phi\_{\vec{j}}(t)\right) + K\_{\rm n}\sin\left(\phi\_{\vec{j}}(t)\right)\right] \tag{3}$$

$$h\_{xy}(t) = \sum\_{j=1}^{N} g\left(\phi\_j(t)\right) \cos\left(\phi\_j(t)\right) \left[K\_l \cos\left(\phi\_j(t)\right) + K\_{ll} \sin\left(\phi\_j(t)\right)\right] \tag{4}$$

() () ( ) ( ) ( ) ( ) ( ) ( ) ( )

() () ( ) ( ) ( ) ( ) () ( )

() () ( ) ( ) ( ) ( ) () ( )

() () ( ) ( ) ( ) ( ) () ( )

*t* represents the angular position of tooth *j,* as shown in Figure 1, and is de-

*xx j j t j n j*

*xy j j t j n j*

*yx j j t j n j*

*yy j j t j n j*

*ht g t t K t K t*

where *Kt* and *Kn* represent the tangential and normal cutting force coefficients, respec-

*ht g t t K t K t*

*ht g t t K t K t*

*ht g t t K t K t*

sin cos sin

 φ

= + (3)

cos cos sin ( )

sin sin cos ( )

cos sin cos ( )

= −+ (6)

 φ

= −+ (5)

 φ

= + (4)

 φ  φ

> φ

> > φ

> > > φ

$$h\_{\rm yx}(t) = \sum\_{j=1}^{N} \mathcal{g}\left(\phi\_{j}(t)\right) \sin\left(\phi\_{j}(t)\right) \left[-\mathcal{K}\_{l} \sin\left(\phi\_{j}(t)\right) + \mathcal{K}\_{n} \cos\left(\phi\_{j}(t)\right)\right] \tag{5}$$

and and is defined as follows:

where

and

tively. ( ) *<sup>j</sup>* φ

φ

φ

cess, arccos(2 1) *st*

$$h\_{\mathcal{Y}}(t) = \sum\_{j=1}^{N} \mathcal{g}\left(\phi\_{\hat{j}}(t)\right) \cos\left(\phi\_{\hat{j}}(t)\right) \left[-K\_{\text{f}} \sin\left(\phi\_{\hat{j}}(t)\right) + K\_{\text{ll}} \cos\left(\phi\_{\hat{j}}(t)\right)\right] \tag{6}$$

where *K<sup>t</sup>* and *K<sup>n</sup>* represent the tangential and normal cutting force coefficients, respectively. *φj*(*t*) represents the angular position of tooth *j,* as shown in Figure 1, and is defined as follows: where φ*st* and φ*ex* represent the start and exit cutting angles of the tooth *j*, respectively. In the up-milling process, 0 φ*st* = and arccos 1 2 ( ) *ex* φ= − *a D* ; in the down-milling pro-

$$\phi\_{\bar{l}}(t) = \frac{2\pi v}{60}t + \frac{2\pi(\bar{j} - 1)}{N} \tag{7}$$

**Figure 1.** Dynamic system of end-milling with two degrees of freedom. **Figure 1.** Dynamic system of end-milling with two degrees of freedom.

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

φφ

φφ

φφ

φφ

1

1

1

=

1

=

*j*

*N*

*j*

=

*j*

*N*

=

*N*

*j*

*N*

**3. Composite Cotes-based method**  *g φj*(*t*) is a piecewise function used to judge whether the tooth *j* is cutting or not and is defined as follows:

$$\mathcal{g}\left(\phi\_{\dot{\!\!/}}(t)\right) = \begin{cases} 1 & \phi\_{\text{st}} < \phi\_{\dot{\!\!/}}(t) < \phi\_{\text{ex}} \\ 0 & \text{otherwise} \end{cases} \tag{8}$$

where *φst* and *φex* represent the start and exit cutting angles of the tooth *j*, respectively. In the up-milling process, *φst* = 0 and *φex* = arccos(1 − 2*a*/*D*); in the down-milling process, *φst* = arccos(2*a*/*D* − 1) and *φex* = π, where *a*/*D* represents the radial immersion ratio.

### **3. Composite Cotes-Based Method**

*3.1. State-Space Expression*

**x**(*t*) can be defined as **x**(*t*) = **q**(*t*) **M** . **<sup>q</sup>**(*t*) <sup>+</sup> **Cq**(*t*)/2 . Then, the state expression of (1) can be expressed as follows:

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

where

$$\mathbf{A} = \begin{bmatrix} -\mathbf{M}^{-1}\mathbf{C}/2 & \mathbf{M}^{-1} \\ \mathbf{C}\mathbf{M}^{-1}\mathbf{C}/4 - \mathbf{K} & -\mathbf{C}\mathbf{M}^{-1}/2 \end{bmatrix} \tag{10}$$

$$\mathbf{B}(t) = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{K}\_{\mathbf{c}}(t) & \mathbf{0} \end{bmatrix} \tag{11}$$

If *<sup>a</sup>p***B**(*t*)[**x**(*t*) <sup>−</sup> **<sup>x</sup>**(*<sup>t</sup>* <sup>−</sup> *<sup>T</sup>*)] is regarded as the inhomogeneous term of . **x**(*t*) = **Ax**(*t*), then the following equation can be derived:

$$\mathbf{x}(t) = e^{\mathbf{A}(t-t\_0)}\mathbf{x}(t\_0) + a\_p \int\_{t\_0}^{t} \left\{ e^{\mathbf{A}(t-\xi)} \mathbf{B}(\xi) [\mathbf{x}(\xi) - \mathbf{x}(\xi - T)] \right\} d\xi \tag{12}$$

where *t*<sup>0</sup> represents the start time. It is worth noting that **x**(*ξ*) is unknown such that (12) is a Volterra integral equation of the second kind.

### *3.2. Numerical Algorithm*

To solve (9) using the composite Cotes formula, the integral equation *f*(*x*) is assumed to be continuous in the interval [*c*, *d*]. Then, [*c*, *d*] is divided into *m* isometric intervals; that is, the span of each interval *h* is equal to (*d* − *c*)/*m*. Based on the isometric nodes *u<sup>k</sup>* = *c* + (*k* − 1)/*h*, where *k* = 1, 2, . . . , *m* + 1, we obtain the following:

$$I\_m = (d - c) \sum\_{k=0}^{m} \mathbb{C}\_k^{(m)} f(u\_k) \tag{13}$$

where *C* (*m*) *k* represents the Cotes coefficients [18], expressed as follows:

$$\mathbb{C}\_{k}^{(m)} = \frac{(-1)^{m-k}}{nk!(n-k)} \int\_{0}^{m} \prod\_{j=0, j\neq k}^{m} (t-j) \mathbf{d}t , k = 0, 1, \ldots, m \tag{14}$$

It is worth noting that *C* (*m*) *k* only depends on *m*. When *m* = 4, (14) can be expressed as follows:

$$\mathcal{C} = \frac{1}{90} [7f(u\_0) + 32f(u\_1) + 12f(u\_2) + 32f(u\_3) + 7f(u\_4)] \tag{15}$$

where *T* consists of the durations of the free and excited vibrations [15]. The durations of the free vibration and excited vibration can be defined as *t<sup>f</sup>* and *tc*, respectively. For the free vibration, **B**(*ξ*) = 0, and (12) can be expressed as follows:

$$\mathbf{x}(t) = e^{\mathbf{A}(t-t\_0)}\mathbf{x}(t\_0) \tag{16}$$

Therefore, the state equation can be solved when *t = t*<sup>0</sup> + *t<sup>f</sup>* , and it can be expressed as follows:

$$\mathbf{x}\left(t\_0 + t\_f\right) = e^{\mathbf{A}t\_f}\mathbf{x}(t\_0)\tag{17}$$

For the excited vibration, its time lies in [*t*<sup>0</sup> + *t<sup>f</sup>* , *t*<sup>0</sup> + *T*]. The interval is divided into *n* isometric intervals, and the span of each interval *l* is equal to (*T* − *t<sup>f</sup>* )/*n*. The discrete time can be expressed as *t<sup>i</sup>* = *t*<sup>0</sup> + *t<sup>f</sup>* + (*i* − 1)*h*, where *i* = 1, 2, . . . , *n* + 1. According to the numerical integral solution of the Volterra integral equation of the second kind [19], we obtain the following:

$$\mathbf{x}(t\_{i+1}) = e^{\mathbf{A}(t\_{i+1} - t\_i)}\mathbf{x}(t\_i) + a\_p \int\_{t\_i}^{t\_{i+1}} \left\{ e^{\mathbf{A}(t\_{i+1} - \xi)} \mathbf{B}(\xi) [\mathbf{x}(\xi) - \mathbf{x}(\xi - t)] \right\} d\xi \tag{18}$$

If [*t<sup>i</sup>* , *t<sup>i</sup>* + 1] is divided into four isometric intervals, then *t<sup>i</sup>* , *ti+*1/4, *ti+*1/2, *ti+*3/4, and *ti+*<sup>1</sup> can be derived. Combining (13) and (15), the integral equation can be expressed as follows:

$$\int\_{t\_i}^{t\_{i+1}} f(\xi)d\xi = \frac{t\_{i+1} - t\_i}{90} [7f(t\_i) + 32f(t\_{i+1/4}) + 12f(t\_{i+1/2}) + 32f(t\_{i+3/4}) + 7f(t\_{i+1})] \tag{19}$$

Upon substituting (19) into (18), we can see that **x**(*ti*+1/4), **x**(*ti*+1/2), and **x**(*ti*+3/4) cannot be solved directly. The barycentric Lagrange interpolation method [20] is introduced here to approximate **x**(*ti*+1/4), **x**(*ti*+1/2), and **x**(*ti*+3/4) and obtain the following:

$$\mathbf{x}(t\_{i+1/4}) \simeq \frac{2\mathbf{1}\mathbf{x}(t\_i) + 1\mathbf{4}\mathbf{x}(t\_{i+1}) - 3\mathbf{x}(t\_{i+2})}{32} \tag{20}$$

$$\mathbf{x}(t\_{i+1/2}) \simeq \frac{\mathbf{3}\mathbf{x}(t\_i) + \mathbf{6}\mathbf{x}(t\_{i+1}) - \mathbf{x}(t\_{i+2})}{8} \tag{21}$$

$$\mathbf{x}(t\_{i+3/4}) \cong \frac{5\mathbf{x}(t\_i) + 30\mathbf{x}(t\_{i+1}) - 3\mathbf{x}(t\_{i+2})}{32} \tag{22}$$

Substituting (19)–(22) into (18), **x**(*t<sup>i</sup>* ) can be expressed as follows:

$$\begin{cases} \mathbf{x}(t\_{i+1}) &= \varepsilon^{\mathbf{A}h} \mathbf{x}(t\_i) + \frac{a\_p h}{6} \left\{ \frac{5}{2} \varepsilon^{\mathbf{A}h} \mathbf{B}(t\_i) (\mathbf{x}(t\_i) - \mathbf{x}(t\_i - T)) + 4 \mathbf{B}(t\_{i+1}) (\mathbf{x}(t\_{i+1}) - \mathbf{x}(t\_{i+1} - T)) \\ &- \frac{1}{2} \varepsilon^{-\mathbf{A}h} \mathbf{B}(t\_{i+2}) (\mathbf{x}(t\_{i+2}) - \mathbf{x}(t\_{i+2} - T)) \right\} \end{cases} \tag{23}$$

It is worth noting that **B**(*tn*+2) cannot be calculated directly using (9) when *i* = *n*. A Newton Cotes formula [13] is introduced here to calculate **x**(*tn*+1) and is expressed as follows:

$$\mathbf{x}(t\_{n+1}) = e^{\mathbf{A}h}\mathbf{x}(t\_n) + \frac{a\_p h}{2} \left[ e^{\mathbf{A}h} \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{24}$$

Combining (17), (23), and (24), the transition matrix can be constructed as follows:

$$\mathbf{C}\_{1} = \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ -e^{\mathbf{A}h} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -e^{\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^{\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{25}$$

$$\mathbf{D}\_{1} = \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ -\frac{75}{2}e^{\mathbf{A}\mathbf{i}}\mathbf{B}\_{1} & -60\mathbf{B}\_{2} & \frac{15}{2}e^{-\mathbf{A}\mathbf{i}}\mathbf{B}\_{3} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -\frac{75}{2}e^{\mathbf{A}\mathbf{i}}\mathbf{B}\_{2} & -60\mathbf{B}\_{3} & \cdots & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & -60\mathbf{B}\_{l-1} & \frac{15}{2}e^{-\mathbf{A}\mathbf{i}}\mathbf{B}\_{l} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & -\frac{75}{2}e^{\mathbf{A}\mathbf{i}}\mathbf{B}\_{l-1} & -60\mathbf{B}\_{l} & \frac{15}{2}e^{-\mathbf{A}\mathbf{i}}\mathbf{B}\_{l+1} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & 45e^{\mathbf{A}\mathbf{B}\_{l}} & 45\mathbf{B}\_{l+1} \end{bmatrix} \tag{26}$$

$$\mathbf{E}\_1 = \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{27}$$

where **B***<sup>i</sup>* represents **B**(*t<sup>i</sup>* ), *i* = 1, 2, . . . , *n* + 1.

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

$$\mathbf{x}\left(\mathbf{I}\_{1} + \mathbf{C}\_{1} + \frac{a\_{p}h}{90}\mathbf{D}\_{1}\right)\begin{bmatrix}\mathbf{x}(t\_{1})\\\mathbf{x}(t\_{2})\\\vdots\\\mathbf{x}(t\_{n+1})\end{bmatrix} = \begin{pmatrix}a\_{p}h\\90\\90\end{pmatrix}\mathbf{D}\_{1} + \mathbf{E}\_{1}\begin{bmatrix}\mathbf{x}(t\_{1}-T)\\\mathbf{x}(t\_{2}-T)\\\vdots\\\mathbf{x}(t\_{n+1}-T)\end{bmatrix} \tag{28}$$

where **I**<sup>1</sup> is a (2*n* + 1) × (2*n* + 1) identity matrix.

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

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

It is worth noting that the chatter stability must be determined using the Floquet theory. If the modulus of any eigenvalue in **Φ**<sup>1</sup> exceeds 1, the system is unstable; if the moduli of all eigenvalues in **Φ**<sup>1</sup> are less than 1, the system is stable [21,22]. Therefore, the boundary curve between the unstable and stable regions in the lobe diagram can be used as the criterion for judging whether chatter occurs.

### *3.3. Simulation and Analysis*

In this section, SDM [10] and SEM [16] are compared. For objective comparison, the CCM, SDM, and SEM share the same parameters and machining conditions. A benchmark example of the single-DOF milling model is utilized to validate and analyze the CCM. The single-DOF milling mathematical model is expressed as follows [23]:

$$
\ddot{\mathbf{x}}(t) + 2\zeta\omega\_n\dot{\mathbf{x}}(t) + \omega\_n^2 \mathbf{x}(t) = -\frac{a\_ph(t)}{m\_t}(\mathbf{x}(t) - \mathbf{x}(t-T))\tag{30}
$$

where *ζ* denotes the relative damping, *ω<sup>n</sup>* denotes the angular natural frequency, and *h*(*t*) denotes the specific cutting force coefficient.

The state-space for single-DOF milling mathematical model is expressed as follows:

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

where

$$\mathbf{A} = \begin{bmatrix} -\mathcal{J}\omega\_n & \frac{1}{m\_t} \\ m\_l \omega\_n^2 \left(\zeta^2 - 1\right) & -\mathcal{J}\omega\_n \end{bmatrix} \tag{32}$$

$$\mathbf{B}(t) = \begin{bmatrix} 0 & 0 \\ -h(t) & 0 \end{bmatrix} \tag{33}$$

where *h*(*t*) is equal to *hxx*(*t*), defined in (3).

The parameters in (32) and (33) are defined as follows: *f<sup>n</sup>* = 922 Hz, *ξ* = 0.011, *<sup>m</sup><sup>t</sup>* = 0.3993 kg, *<sup>K</sup><sup>t</sup>* = 6 <sup>×</sup> <sup>10</sup><sup>8</sup> N/m<sup>2</sup> and *<sup>K</sup><sup>t</sup>* = 2 <sup>×</sup> <sup>10</sup><sup>8</sup> N/m<sup>2</sup> . The adopted machining condition is down-milling. The radial immersion ratio *a*/*D* is set as 1 to avoid intermittent milling. The spindle speed is set as 5000 rpm (*v* = 5000 rpm). The depths of cut are set as 0.001 m, 0.0005 mm, and 0.0002 mm, respectively. All programs in this study are executed in MATLAB R2019a and run on a personal computer (AMD Ryzen 5 5600H; CPU 4.0 GHz, 16 GB). The maximum modulus of the eigenvalues of **Φ** is labelled as |*λ*|, and the maximum modulus of the eigenvalues of **Φ** with SDM is labelled as |*λ*0|. In this study, |*λ*0| is treated as the exact value.

The convergence rate comparisons of the CCM, SDM, and SEM are shown in Figure 2. It is seen from the figure that when *n* is small, the convergence rate of the CCM is significantly higher than those of the SDM and SEM regardless of the depth of cut. As *n* increases, the convergence rates of the three methods approach gradually, and the convergence rates of the CCM and SEM are higher than that of SDM. The figure (a) is enlarged to observe the

differences better. It is worth noting that the convergence rate of the CCM is significantly higher than those of the SDM and SEM when the discrete number changes from 75 to 100. The results demonstrate that the convergence rate of the CCM outweighs those of the SDM and SEM. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 8 of 25

**Figure 2.** Convergence rate comparisons of the CCM, SDM and SEM when *a*/*D* = 1. (**a**) *ap* = 0.0002 m and *n*∈[30,100]; (**b**) *ap* = 0.0002 m and *n*∈[75,100]; (**c**) *ap* = 0.0005 m and *n*∈[30,100]; (**d**) *ap* = 0.001 m and *n*∈[30,100]. **Figure 2.** Convergence rate comparisons of the CCM, SDM and SEM when *a*/*D* = 1. (**a**) *ap* = 0.0002 m and *n* ∈ [30, 100]; (**b**) *a<sup>p</sup>* = 0.0002 m and *n* ∈ [75, 100]; (**c**) *a<sup>p</sup>* = 0.0005 m and *n* ∈ [30, 100]; (**d**) *a<sup>p</sup>* = 0.001 m and *n* ∈ [30, 100].

The stability lobe diagrams can be used to compare the calculation accuracy among the CCM, SDM, and SEM. As a reference, the discrete number is set as 500. The other parameters are set as follows: the discrete numbers are set as 25, 40, and 55; the radial immersion ratios are set as 0.05, 0.5, and 1; the spindle speed varies from 0.5×104 to 2.5×104 rpm, and 200 equally distributed sampling points are selected within this range; the depth of cut varies from 0 to 0.01 m, and 100 equally distributed sampling points are selected within this range. The stability lobe diagrams can be used to compare the calculation accuracy among the CCM, SDM, and SEM. As a reference, the discrete number is set as 500. The other parameters are set as follows: the discrete numbers are set as 25, 40, and 55; the radial immersion ratios are set as 0.05, 0.5, and 1; the spindle speed varies from 0.5 <sup>×</sup> <sup>10</sup><sup>4</sup> to 2.5 <sup>×</sup> <sup>10</sup><sup>4</sup> rpm, and 200 equally distributed sampling points are selected within this range; the depth of cut varies from 0 to 0.01 m, and 100 equally distributed sampling points are selected within this range.

In the single-DOF system, the stability lobe diagrams obtained with the CCM, SDM, and SEM when *a*/*D* = 1 are given in Table 1. It can be seen that the stability lobe diagrams of the CCM are closer to the reference ones than the SDM and SEM. Note that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25. Overall, the CCM is more accurate than the SDM and SEM. To further compare the calculation accuracy, two indicators are introduced, i.e., the arithmetic mean of In the single-DOF system, the stability lobe diagrams obtained with the CCM, SDM, and SEM when *a*/*D* = 1 are given in Table 1. It can be seen that the stability lobe diagrams of the CCM are closer to the reference ones than the SDM and SEM. Note that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25. Overall, the CCM is more accurate than the SDM and SEM. To further compare the calculation accuracy, two indicators are introduced, i.e., the arithmetic mean of

relative error (*AMRE*) and mean squared error (*MSE*) [2]. The *AMRE* and *MSE* are defined

*AMRE*

0

<sup>−</sup> <sup>=</sup> (34)

1 0 1 *<sup>n</sup> i i i i*

*n a* <sup>=</sup>

*a a*

as follows [2]:

relative error (*AMRE*) and mean squared error (*MSE*) [2]. The *AMRE* and *MSE* are defined as follows [2]: *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25

*AMRE* = 1 *n n* ∑ *i*=1 |*a<sup>i</sup>* − *ai*0| *ai*0 *MSE* = 1 *n n* ∑ *i*=1 (*a<sup>i</sup>* − *ai*0) 2 *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 25 ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i iMSE a a n* <sup>=</sup> = − (35) ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. ( )<sup>2</sup> 0 1 1 *<sup>n</sup> i i i MSE a a n* <sup>=</sup> = − (35) where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. ( )<sup>2</sup> 0 11 *n i i i MSE a a n*<sup>=</sup> = − (35) where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively.

where *a<sup>i</sup>* and *ai*<sup>0</sup> denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 <sup>×</sup> <sup>10</sup>−<sup>8</sup> , 4.24 <sup>×</sup> <sup>10</sup>−<sup>8</sup> , and 1.07 <sup>×</sup> <sup>10</sup>−<sup>7</sup> , respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. where *ai* and *ai*0 denote the predicted axial depth and reference axial depth, respectively. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. *<sup>n</sup>*denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* <sup>=</sup> 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM. *n* denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the *AMRE* and *MSE* are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the *AMRE* and *MSE* of the SEM cannot be attained. Take *n* = 40 as an example. The *AMRE* obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The *MSE* obtained with the CCM, SDM, and SEM is 2.62 × 10−8, 4.24 × 10−8, and 1.07× 10−7, respectively. The results show that the CCM is closer to the reference values than the SDM and SEM.

**Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. **Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. **Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. **Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. **Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. **Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. **Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. *n* **25 40 55 Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. *n* **25 40 55 Table 1.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 1. *n* **25 40 55** 

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe

the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* <sup>=</sup> 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

Furthermore, the computational efficiency is an important criterion for measuring

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06

vantages in the computational efficiency.

vantages in the computational efficiency.

**Figure 3.** *AMRE* and *MSE* of the CCM, SDM, and SEM when *a*/*D* = 1. (**a**) *AMRE*; (**b**) *MSE.*  **Figure 3.** *AMRE* and *MSE* of the CCM, SDM, and SEM when *a*/*D* = 1. (**a**) *AMRE*; (**b**) *MSE*.

Furthermore, the computational efficiency is an important criterion for measuring the effectiveness of the chatter analysis methods. Figure 4 shows the computational time of the CCM, SDM, and SEM under different discrete numbers when *a*/*D* = 1. We observe that the computational time can be reduced from 55.63 s to 20.21 s, 62.45 s to 25.79 s, 72.06 s to 30.13 s, and 81.21 s to 35.18 s when the discrete numbers are 40, 45, 50, and 55, respectively. Furthermore, the stability lobe diagrams, *AMRE* and *MSE,* as well as the computational time when *a*/*D* = 0.05 and *a*/*D* = 0.5 are given in Tables 2 and 3 and Figures 5 and 6. It could be noted that the stability lobe diagrams obtained with the CCM are closer to the reference ones than those obtained with the SDM and SEM, and the CCM has the highest calculation accuracy. The results also demonstrate that the CCM has significant advantages in the computational efficiency. **Figure 3.** *AMRE* and *MSE* of the CCM, SDM, and SEM when *a*/*D* = 1. (**a**) *AMRE*; (**b**) *MSE.* 

s to 30.13 s, and 81.21 s to 35.18 s when the discrete numbers are 40, 45, 50, and 55, respectively. Furthermore, the stability lobe diagrams, *AMRE* and *MSE,* as well as the computational time when *a*/*D =* 0.05 and *a*/*D =* 0.5 are given in Tables 2 and 3 and Figures 5 and 6. It could be noted that the stability lobe diagrams obtained with the CCM are closer to the reference ones than those obtained with the SDM and SEM, and the CCM has the highest calculation accuracy. The results also demonstrate that the CCM has significant ad-

s to 30.13 s, and 81.21 s to 35.18 s when the discrete numbers are 40, 45, 50, and 55, respectively. Furthermore, the stability lobe diagrams, *AMRE* and *MSE,* as well as the computational time when *a*/*D =* 0.05 and *a*/*D =* 0.5 are given in Tables 2 and 3 and Figures 5 and 6. It could be noted that the stability lobe diagrams obtained with the CCM are closer to the reference ones than those obtained with the SDM and SEM, and the CCM has the highest calculation accuracy. The results also demonstrate that the CCM has significant ad-

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

**Figure 4.** Computational time of the CCM, SDM, and SEM when *a*/*D* = 1. **Figure 4.** Computational time of the CCM, SDM, and SEM when *a*/*D* = 1.

SDM

SDM

SDM

SEM

SEM

SEM

**Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. **Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. **Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. **Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. **Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. *n* **25 40 55 Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. *n* **25 40 55 Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05. *n* **25 40 55**  *n* **25 40 55**  *n* **25 40 55**  *n* **25 40 55** 

**Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05.

**Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* =0.05.

**Table 2.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.05.

**2022**, *13*, FOR PEER REVIEW 11 of 25

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

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

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

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

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

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

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

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

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5. **Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5. **Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5. **Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

CCM

CCM

CCM

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

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

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

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

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

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

*n* **25 40 55** 

*n* **25 40 55** 

*n* **25 40 55** 

*n* **25 40 55** 

*n* **25 40 55** 

*n* **25 40 55** 

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

**Table 3.** Stability lobe diagrams of the CCM, SDM, and SEM when *a*/*D* = 0.5.

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

**a b**

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

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

**c**

*a*/*D* = 0.5; (**c**) *MSE* when *a*/*D* = 0.05; (**d**) *MSE* when *a*/*D* = 0.5.

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

### **4. Simpson's 3/8-Based Method**

*4.1. State-Space Expression*

To match with S38M, a new state-space expression is defined as follows [17]:

**Figure 5.** *AMRE* and *MSE* of the CCM, SDM, and SEM. (**a**) *AMRE* when *a*/*D* = 0.05; (**b**) *AMRE* when

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

where

**a b**

**d** <sup>7</sup> 10<sup>−</sup> × <sup>8</sup> 10<sup>−</sup> ×

$$\mathbf{A} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix} \tag{35}$$

and

$$\mathbf{B}(t) = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ -\mathbf{M}^{-1}\mathbf{K}\_c(t) & \mathbf{0} \end{bmatrix} \tag{36}$$

Based on the numerical integration method of the Volterra integral equation of the second kind [19], the state-space equation can be deduced as follows:

$$\mathbf{x}(t) = e^{\mathbf{A}(t - t\_{st})} \mathbf{x}(t\_{st}) + a\_p \int\_{t\_{st}}^{t} \left\{ e^{\mathbf{A}(t - \tilde{\boldsymbol{\xi}})} \mathbf{B}(\tilde{\boldsymbol{\xi}}) [\mathbf{x}(\tilde{\boldsymbol{\xi}}) - \mathbf{x}(\tilde{\boldsymbol{\xi}} - \boldsymbol{T})] \right\} d\boldsymbol{\xi} \tag{37}$$

where *tst* represents the start time and defaults to *t*1.

When *t = t*3, we obtain the following:

$$\mathbf{x}(t\_3) = e^{\mathbf{A}(t\_3 - t\_1)}\mathbf{x}(t\_1) + a\_p \int\_{t\_M}^{t\_1} \left\{ e^{\mathbf{A}(t\_1 - \xi)} \mathbf{B}(\xi) [\mathbf{x}(\xi) - \mathbf{x}(\xi - T)] \right\} d\xi \tag{38}$$
