Experimental Research and Significance Analysis of Advanced Interpolation Methods for Optimizing System State Items and Processing Parameters

Abstract

The Finite-Difference Method (FDM) plays a pivotal role in the field of stability prediction, particularly in the modeling and stability analysis of cutting process dynamics. However, traditional approaches to optimizing the FDM often treat system state terms and time-delay terms as a monolithic entity, failing to explicitly distinguish between them, which leads to a lack of specificity in selecting optimization targets. In this study, an innovative approach is introduced by incorporating the third-order Newton interpolation method and the fourth-order Hermite interpolation method. By comparing the computational accuracy and convergence speed, it is found that the 3N-FDM (third-order Newton Finite-Difference Method) exhibits superior overall performance, and it is clearly pointed out that increasing the order does not always result in better outcomes. Additionally, this study selects different discretization numbers, denoted as m, for comparative analysis to thoroughly evaluate their impact on computational accuracy. Experimental validation demonstrates the high accuracy of the 3N-FDM. Through a one-way ANOVA (analysis of variance) of tool wear and workpiece surface roughness, it is revealed that changes in system state terms have the most significant impact on the feed rate $f$, followed by the cutting depth $a_p$, and finally the spindle speed $n$. Based on the experimental results and analysis mentioned above, this study concludes that optimizing system state terms can more effectively explore the combined influences of processing parameters on processing quality, production efficiency, and tool wear.

1. Introduction

Machining vibrations can primarily be categorized into two types: forced vibrations and self-excited vibrations [1]. Forced vibrations are typically caused by external harmonic excitation. Once identified and the source of excitation is determined, appropriate control measures can often effectively prevent, mitigate, or even completely eliminate such vibrations. For self-excited vibrations, according to their different generation mechanisms, chatter can be further refined into the following four categories: force-thermal chatter, modal-coupling chatter, friction chatter, and regenerative chatter [2]. In the field of milling, regenerative chatter is widely regarded as the dominant factor leading to instability in the machining process [3].

The root cause of regenerative chatter lies in the interaction between the milling force and the wavy cutting marks left on the workpiece surface by the preceding tool tooth. Based on this principle, a theoretical cutting diagram is depicted, as shown specifically in Figure 1. Upon completion of the cutting process, tool tooth 4 leaves wavy cutting marks on the workpiece surface. When the subsequent tool tooth 1 engages in cutting, these wavy marks cause fluctuations in the chip thickness and cutting force, which in turn induce vibrations in the tool. This vibration forms a closed-loop feedback mechanism with the cutting process, known as the regenerative effect, and it is precisely this effect that ultimately results in the occurrence of milling chatter [4].

From a mathematical perspective, the milling dynamic process considering the regenerative effect can be modeled as a time-delay differential equation system with infinite-dimensional characteristics that includes a time-periodic matrix. In this system, the time-periodic matrix plays a crucial role as it determines the stable state of the milling process [5]. Therefore, the problem of determining the stability boundary essentially boils down to solving the corresponding time-delay differential equation of the system [6].

The stability prediction methods for this problem currently include various approaches, such as frequency domain methods, discrete methods, and numerical integration methods. Among these, the full-discrete method occupies a dominant position in the field of stability prediction. The full-discrete method cleverly uses a finite-dimensional transfer matrix to approximate the originally infinite-dimensional single-valued operator, thereby significantly reducing the complexity and computation time of the solution. This method not only has the advantage of high computational accuracy but also high efficiency and wide applicability, thus demonstrating significant advantages in practical applications.

The full-discrete method is widely adopted in the field of stability prediction due to its strategy of using finite-dimensional transfer matrices to approximate infinite-dimensional single-valued operators, which significantly reduces the solution complexity and computation time. It offers high computational accuracy, high efficiency, and wide applicability, demonstrating prominent advantages in practical applications. The full-discrete method evolved from the semi-discrete method to address the issue of insufficient computational efficiency. Ding et al. pioneered the use of second-order Lagrange polynomials for interpolating the system state items, known as the second-order FDM [7]. Later, Tang et al. further applied second-order Lagrange polynomial interpolation to the time-lag items, enhancing the computational accuracy of stability prediction [8].

To further optimize the computational performance of discrete methods, researchers continue to explore and adopt various interpolation methods or approximation theories, such as Lagrange polynomials, Newton polynomials, Hermite polynomials, and the least-squares method, to estimate the system’s time-lag items and state items, thereby constructing higher-order full-discrete methods. Specifically, second-order Lagrange polynomials, second-order Newton polynomials, and third-order Hermite polynomials have been used to estimate the system state items, while the corresponding time-lag items are interpolated with second-order Lagrange polynomials, second-order Hermite polynomials, and third-order Newton polynomials, respectively [9,10,11,12,13,14]. Additionally, Ozoegwu et al. used the least-squares method to predict the stability lobe diagram, achieving a convergence accuracy of up to the fourth order [15,16]; Zhou et al. employed fourth-order Lagrange polynomials to interpolate the time-lag items, proposing the fourth-order FDM [17].

However, it is not difficult to find from the above studies that researchers often choose to use higher-order interpolation methods or approximation theories in pursuit of higher computational accuracy. While this does enhance accuracy to some extent, it also causes the computational workload to increase exponentially, leading to extended computation times. Therefore, there is an urgent need to find new optimization ideas to improve the accuracy of the full-discrete method.

During the optimization of the full-discrete method, a commonly overlooked issue is the accuracy of the discrete number m. When m is too small, the computational accuracy significantly decreases because too few discrete points may not adequately capture the detailed characteristics of the system’s dynamic behavior, introducing significant errors. These errors can accumulate and amplify, eventually leading to prediction results deviating from reality. Conversely, when m is too large, although it can improve computational accuracy to some extent, it also leads to a significant reduction in computational efficiency. Too many discrete points increase the computational workload, making the solution process complex and time-consuming, not only increasing computational costs but also potentially limiting the feasibility and practicality of the algorithm in actual applications. Therefore, the value of the discrete number m has a vital impact on the overall computational performance.

More importantly, the current optimization of the full-discrete method mainly targets the system state items and system time-lag items, but it does not optimize according to actual needs. These optimization measures mainly focus on improving the computational accuracy of the full-discrete method, neglecting the impact of the system state items and system time-lag items themselves on stability.

In response to the above issues, this paper optimized the system state items using high-order interpolation methods, analyzed the impact of different interpolation methods on computational accuracy, and determined the optimal full-discrete method. Additionally, this paper investigated the influence of different discrete numbers m on accuracy and computational efficiency, summarized general patterns, and aimed to maximize predictive accuracy while maintaining computational efficiency. Furthermore, this paper employed experimental methods to study the relationship between processing parameters and system state items and summarized the advantages brought by optimizing the system state items.

2. Milling Stability Prediction Based on the Improved Full-Discrete Method

The full-discrete method ingeniously employs a linear periodic system model with a single discrete time lag to accurately describe the dynamic process of milling. In this method, the original system is discretized within a single period. Within each discrete interval, the system state items, system time-lag items, and periodic coefficient items are precisely approximated using linear interpolation techniques, and the system response is obtained directly through integration. According to the system transfer matrix and theory, the stability of the system can be predicted.

The milling dynamics equation considering regenerative chatter can be expressed in an n-dimensional state space as follows:

$$\dot{U}\left(t\right)=AU\left(t\right)+B\left(t\right)\left[U\left(t\right)-U\left(t-\tau \right)\right]$$

(1)

$A$ is the constant matrix representing the time-invariant nature of the system, and $B$ is the periodic coefficient matrix, satisfying the condition $B\left(t+T\right)=B\left(t\right)$.

To discretize the time delay $T$ into $m$ equal small time intervals, denoted as $T=m\tau , m\in N^{*}$, we obtain a series of small time segments $k\tau -n\tau \le t\le k\tau +n\tau$ (where $k$ is a natural number) on each of which, under the initial condition $z_k=z\left(k\tau \right)$, the response of Equation (1) is given by the following:

$$z\left(t\right)=e^{A_0\left(t-k\tau \right)}z\left(k\tau \right)+{\displaystyle {\int }_{k\tau }^t\left(e^{A_0\left(t-\xi \right)}\left(A\left(\xi \right)z\left(\xi \right)+B\left(\xi \right)z\left(\xi -T\right)\right)\right)d\xi }$$

(2)

When the time is $t=\left(k+1\right)\tau$, the corresponding response $x_{k+1}$ is described as follows:

$$x_{k+1}=e^{A_{0}r}x\left(k\tau \right)+{\displaystyle {\int }_0^{k\tau }\left\{e^{A_0\xi }\left[\begin{array}{c}A\left(k\tau +\tau -\xi \right)x\left(k\tau +\tau -\xi \right)+\\ A\left(k\tau +\tau -\xi \right)x\left(k\tau +\tau -\xi -T\right)\end{array}\right]\right\}d\xi }$$

(3)

Then, different interpolating polynomials are used to approximate the system state item $x\left(k\tau +\tau -\xi \right)$ and the system periodic coefficient item $A\left(k\tau +\tau -\xi \right)$ and $B\left(k\tau +\tau -\xi \right)$, as well as the time-lag item $x\left(k\tau +\tau -\xi -T\right)$.

2.1. Applying Higher-Order Interpolating Polynomials to Optimize the System State Item

The system state item is a comprehensive indicator that reflects the current status of the system in real time, encompassing a range of fundamental physical parameters such as displacement, velocity, and acceleration. In the field of mechanical processing, this concept is concretized as the feed rate $f$, the spindle speed, and the cutting trajectory of the tool on the workpiece surface. These parameters collectively provide a detailed description of the dynamic behavior of the system during the mechanical processing process, playing a pivotal role in ensuring processing accuracy, improving production efficiency, and optimizing process parameters.

Currently, for the optimization of the system state item, second-order Lagrange interpolating polynomials and second-order Newton interpolating polynomials are commonly used as approximation tools. However, to further explore the potential for optimization, this study employs third-order Newton interpolating polynomials and fourth-order Hermite interpolating polynomials to approximate the system state item.

This study aims to compare the differences between these two interpolating polynomial approximation methods in terms of computational accuracy and convergence speed, analyze the impact of the discrete number m on processing accuracy and computation time, and summarize the general patterns. Furthermore, exploring whether using higher-order interpolating polynomials can lead to more excellent computational methods. Through this exploration, it is expected to reveal the potential of higher-order interpolating polynomials in enhancing computational efficiency, providing new insights and strategies for optimizing numerical calculation methods, and improving computational efficiency. Finally, the optimal method is experimentally validated to verify the accuracy of its theoretical calculation results. During the trial process, different parameter combinations are selected for testing, and whether chatter occurs during processing is judged through cutting force curves and data.

2.1.1. Approximating the State Item with the Third-Order Newton Interpolation Method

The third-order Newton polynomial, selecting four endpoint values, namely, $x_{k+1}$, $x_k$, $x_{k-1}$, and $x_{k-2}$ is used to interpolate the system state item $x\left(k\tau +\tau -\xi \right)$ in Equation (3), yielding the following expression:

$$x\left(k\tau +\tau -\xi \right)=a_0{x}_{k+1}+b_0{x}_k+c_0{x}_{k-1}+d_0{x}_{k-2}$$

(4)

$$\begin{array}{l}{a}_0=-{\displaystyle \frac{{\xi }^3}{6{\tau }^3}}+{\displaystyle \frac{{\xi }^2}{{\tau }^2}}-{\displaystyle \frac{11\xi }{6\tau }}+1b_0={\displaystyle \frac{{\xi }^3}{2{\tau }^3}}-{\displaystyle \frac{5{\xi }^2}{2{\tau }^2}}+{\displaystyle \frac{3\xi }{\tau }}\\ c_0=-{\displaystyle \frac{{\xi }^3}{2{\tau }^3}}+{\displaystyle \frac{2{\xi }^2}{{\tau }^2}}-{\displaystyle \frac{3\xi }{2\tau }}{d}_0={\displaystyle \frac{{\xi }^3}{6{\tau }^3}}-{\displaystyle \frac{{\xi }^2}{2{\tau }^2}}+{\displaystyle \frac{\xi }{3\tau }}\end{array}$$

(5)

The first-order Lagrange polynomial is used to approximate the time-lag item $x\left(k\tau +\tau -\xi -T\right)$ in Equation (3), which is specifically described as follows:

$$x\left(k\tau +\tau -\xi -T\right)={\displaystyle \frac{\xi }{2}}{x}_{k-m}+(1-{\displaystyle \frac{\xi }{2}})x_{k-m+1}$$

(6)

The first-order Lagrange polynomial uses the endpoints of the interval to approximate the periodic coefficient items $A\left(k\tau +\tau -\xi \right)$ and $B\left(k\tau +\tau -\xi \right)$ in Equation (3), which can be expressed as follows:

$$\begin{array}{l}A\left(k\tau +\tau -\xi \right)\doteq A_{k+1}+\left(A_k-A_{k+1}\right)\xi /\tau \\ B\left(k\tau +\tau -\xi \right)\doteq B_{k+1}+\left(B_k-B_{k+1}\right)\xi /\tau \end{array}$$

(7)

In the formula, $A_k$ and $A_{k+1}$ represent the endpoint values at times $t=k\tau$ and $t=\left(k+1\right)\tau$, respectively.

Subsequently, substituting Equations (4), (6), and (7) into Equation (3), the following is obtained:

$$\left(1-F_{k1}\right)x_{k+1}=\left(F_0+F_k\right)x_k+F_{kp1}{x}_{k-1}+F_{kp2}{x}_{k-2}+F_{k1m}{x}_{k+1-m}+F_{km}{x}_{k-m}$$

(8)

When matrix ${\left(I-F_{kI}\right)}^{-1}$ exists, the discrete mapping can be represented as follows:

$$D_k=\left[\begin{array}{lllllll}{H}_{11}\hfill & H_{p1}\hfill & H_{p2}\hfill & \cdots \hfill & 0\hfill & H_{1m}\hfill & H_m\hfill \\ I\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & I\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & I\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & I\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & I\hfill & 0\hfill \end{array}\right]$$

(9)

$$y_{k+1}=D_k{y}_k$$

(10)

Then, the state transition matrix $\mathsf{\Phi }$ is represented by a series of matrices $D_k$, which is specifically expressed as follows:

$$y_m=D_{m-1}{D}_{m-2}\cdots D_1{D}_0{y}_0$$

(11)

Finally, based on Equation (11), the flutter stability lobe diagram can be obtained using Floquet theory.

According to the Floquet theorem, the stability boundary of the milling system can be determined, and its judgment criterion is shown in Equation (12).

$$\max \left(\left|\lambda \left({\psi }_1\right)\right|\right)\left\{\begin{array}{c}<1, \mathrm{Stable}\\ =1, \mathrm{Critical}\\ >1, \mathrm{Chatter}\end{array}\right.$$

(12)

2.1.2. Approximating the State Item with the Fourth-Order Hermite Interpolation Method

The fourth-order Hermite interpolating polynomial is used to approximate the state item $X\left(t\right)$ in Equation (3) within the time interval $\left[t_i, t_{i+1}\right]$. The first-order derivative of the response values at any nodes $t=t_i$ and $t=t_{i+1}$ can be expressed as follows:

$$\dot{X}\left(t_i\right)=AX\left(t_i\right)+B\left(t_i\right)\left[X\left(t_i\right)-X\left(t_i-\tau \right)\right]$$

(13)

$$\dot{X}\left(t_{i+1}\right)=AX\left(t_{i+1}\right)+B\left(t_{i+1}\right)\left[X\left(t_{i+1}\right)-X\left(t_{i+1}-\tau \right)\right]$$

(14)

Based on the fourth-order Hermite interpolating polynomial, the state item $x\left(k\tau +\tau -\xi \right)$ can be interpolated and approximated within the time interval $\left[t_i, t_{i+1}\right]$ using the derivative node values $\dot{X}\left(t_i\right)$, $\dot{X}\left(t_{i+1}\right)$ and the state item node values $X\left(t_{i-1}\right)$, $X\left(t_i\right)$, and $X\left(t_{i+1}\right)$ as shown in Equation (3).

$$X\left(t\right)=a_1{X}_{i-1}+b_1{X}_i+c_1{X}_{i+1}+d_1{X}_{i-n}+e_1{X}_{i-n+1}$$

(15)

In the formula, $X_i$ is the abbreviation of $X\left(i\cdot \mathsf{\Delta }t\right)$, and the expressions for $a_1$, $b_1$, $c_1$, $d_1$, and $e_1$ are as follows:

$$\begin{array}{l}{a}_1=\left({\displaystyle \frac{t^4}{4\mathsf{\Delta }{t}^4}}-{\displaystyle \frac{t^3}{2\mathsf{\Delta }{t}^3}}+{\displaystyle \frac{t^2}{4\mathsf{\Delta }{t}^2}}\right)I\\ b_1=I+\left(A+B_i\right)t-{\displaystyle \frac{\left(2I+A\mathsf{\Delta }t+B_i\mathsf{\Delta }t\right)t^2}{\mathsf{\Delta }{t}^2}}-{\displaystyle \frac{\left(A+B_i\right)t^3}{\mathsf{\Delta }{t}^2}}+{\displaystyle \frac{\left(I+A\mathsf{\Delta }t+B_i\mathsf{\Delta }t\right)t^4}{\mathsf{\Delta }{t}^4}}\\ c_i={\displaystyle \frac{7I{t}^2}{4\mathsf{\Delta }{t}^2}}-{\displaystyle \frac{\left(A+B_{i+1}\right)t^2}{2\mathsf{\Delta }t}}+{\displaystyle \frac{I{t}^3}{2\mathsf{\Delta }{t}^3}}+{\displaystyle \frac{\left(2Ah+2B_{i+1}h-5I\right)t^4}{4\mathsf{\Delta }{t}^4}}\\ d_1=-{\displaystyle \frac{B_{i+1}{t}^2}{\mathsf{\Delta }{t}^3}}+{\displaystyle \frac{B_i{t}^3}{\mathsf{\Delta }{t}^2}}+{\displaystyle \frac{B_i{t}^2}{\mathsf{\Delta }t}}-B_{i}t{e}_1={\displaystyle \frac{B_{i+1}{t}^2}{2\mathsf{\Delta }t}}-{\displaystyle \frac{B_{i+1}{t}^4}{2\mathsf{\Delta }{t}^3}}\end{array}$$

(16)

In the formula, $I$ is the identity matrix.

The second-order Newton interpolating polynomial is used to approximate the time-lag item $x\left(k\tau +\tau -\xi -T\right)$ in Equation (3), which is specifically described as follows:

$$x\left(k\tau +\tau -\xi -T\right)=\left({\displaystyle \frac{t^2}{2\mathsf{\Delta }{t}^2}}-{\displaystyle \frac{3t}{2\mathsf{\Delta }t}}\right)X_{i-n}+\left({\displaystyle \frac{2I}{\mathsf{\Delta }t}}-{\displaystyle \frac{t^2}{\mathsf{\Delta }{t}^2}}\right)X_{i-n+1}+\left({\displaystyle \frac{t^2}{2\mathsf{\Delta }{t}^2}}-{\displaystyle \frac{t}{2\mathsf{\Delta }t}}\right)X_{i-n+2}$$

(17)

Linear interpolation is used to approximate the periodic coefficient items $A\left(k\tau +\tau -\xi \right)$ and $B\left(k\tau +\tau -\xi \right)$ in Equation (3), which is specifically described as follows:

$$\begin{array}{l}A\left(k\tau +\tau -\xi \right)={\displaystyle \frac{\mathsf{\Delta }t-t}{\mathsf{\Delta }t}}{A}_i+{\displaystyle \frac{t}{\mathsf{\Delta }t}}{A}_{i+1}\\ B\left(k\tau +\tau -\xi \right)={\displaystyle \frac{\mathsf{\Delta }t-t}{\mathsf{\Delta }t}}{B}_i+{\displaystyle \frac{t}{\mathsf{\Delta }t}}{B}_{i+1}\end{array}$$

(18)

Substituting Equations (15), (17), and (18) into Equation (2), the following is obtained:

$$X_{i+1}=P_i\left[\begin{array}{l}\left(e^{A\mathsf{\Delta }t}+H_{13}{B}_{i+1}+H_{14}{B}_i\right)X_i+\left(H_{15}{B}_{i+1}+H_{16}{B}_i\right)X_{i-1}+\\ \left(H_{17}{B}_{i+1}+H_{18}{B}_i\right)X_{i-n+2}+\left(H_{19}{B}_{i+1}+H_{20}{B}_i\right)X_{i-n+1}+\\ \left(H_{21}{B}_{i+1}+H_{22}{B}_i\right)X_{i-n}\end{array}\right]$$

(19)

From Equation (19), it can be seen that if matrix $P_i$ is non-singular, then Equation (19) can be converted into the following form:

$${\left[X_{i+1}{X}_i{X}_{i-1}\dots X_{i+1-n}\right]}^T=D_k{\left[X_i{X}_i{X}_{i-1}\dots X_{i+1-n}\right]}^T$$

(20)

In the formula, matrix $M_i$ is shown in Equation (21):

$$D_k=\left[\begin{array}{cccccc}{H}_{11}^i& H_{12}^i& \cdots & H_{1, n-1}^i& H_{1, n}^i& H_{1, n+1}^i\\ I& 0& \cdots & 0& 0& 0\\ 0& I& \cdots & 0& 0& 0\\ ⋮& ⋮& & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& I& 0\end{array}\right]$$

(21)

The matrices $H_{11}^i$, $H_{12}^i$, $H_{1,n-1}^i$, $H_{1,n}^i$, and $H_{1,n+1}^i$ in Equation (21) are shown in the following Equation (22):

$$\begin{array}{ll}{M}_{11}^i=P_i\left(e^{A\mathsf{\Delta }t}+H_{13}{B}_{i+1}+H_{14}{B}_i\right)& M_{12}^i=P_i\left(H_{15}{B}_{i+1}+H_{16}{B}_i\right)\\ M_{1, n-1}^i=P_i\left(H_{17}{B}_{i+1}+H_{18}{B}_i\right)& M_{1, n}^i=P_i\left(H_{19}{B}_{i+1}+H_{20}{B}_i\right)\\ M_{1, n+1}^i=P_i\left(H_{21}{B}_{i+1}+H_{22}{B}_i\right)& \end{array}$$

(22)

The stability boundary of the milling system can be determined by the theorem Floquet as shown in Equation (12).

2.2. Simulation Calculation of Third-Order Newton Interpolation Method and Fourth-Order Hermite Interpolation Method

The fully discrete stability prediction method constructed above utilizes the third-order Newton interpolating polynomial (abbreviated as 3N-FDM) and the fourth-order Hermite interpolating polynomial (abbreviated as 4H-FDM) to perform high-precision interpolation approximation of the system’s state item. Here, “N” and “H” serve as identifiers for Newton and Hermite polynomials, respectively, and the digit preceding the symbol explicitly indicates the order of the polynomial.

To thoroughly investigate and determine the optimal order of the improved fully discrete method, different-order prediction methods were used to calculate the convergence and stability lobe diagram of a two-degrees-of-freedom milling dynamics model. Throughout this process, strict adherence to the parameter settings outlined in

Table 1. Tool system parameters.

Parameter	${\mathit{t}}_{\mathit{i}-1}$	${\mathit{t}}_{\mathit{i}}$
Natural Frequency	$f_n$	922 Hz
Damping Ratio	$\zeta$	1.1%
Modal Mass	$m_x$	0.03993 kg
Cutting Force Parameters	$\begin{array}{l}{K}_t\\ K_n\end{array}$	$6\times 10^8\mathrm{N}/{\mathrm{m}}^2$ $2\times 10^8\mathrm{N}/{\mathrm{m}}^2$
Number of Teeth	$N_t$	2

was maintained to ensure the accuracy and reproducibility of the experiments.

To determine the optimal order of the interpolating polynomial for the system state item, a comprehensive comparative analysis of the fully discrete method at different orders is required. To this end, an indicator called the local discrete error $\left|\left|\mu \right|-\left|{\mu }_0\right|\right|$ is introduced, which quantifies the speed of convergence. Here, $\left|\mu \right|$ represents the modulus of the critical eigenvalue of the state transition matrix, while $\left|{\mu }_0\right|$ is the precise value obtained using the second-order FDM (Finite-Difference Method) when m = 200.

When conducting a comparative analysis of the convergence speed of the improved fully discrete method, relevant computational materials were referenced to enhance the accuracy and clarity of the calculation results and accordingly determined the following parameters: the spindle speed was set at $n=5000\mathrm{r}/\mathrm{m}\mathrm{i}\mathrm{n}$, and the radial depth of the cut ratio was selected as $a_p/D=0.1$ (the ratio of the cutting depth $a_p$ to the tool diameter D is 0.1). Given that up-milling is commonly used in the finish machining process of CNC machines, we also followed this principle in our calculations. For the machining of high-temperature alloy materials, cutting depths of 0.2 mm, 0.5 mm, and 1 mm were specified for finish machining, semi-finish machining, and finer stages of machining, respectively. these parameter settings were adopted in this calculation, as detailed in Figure 2, Figure 3 and Figure 4.

Solving Equations (10) and (21) using different-order improved fully discrete methods yields curves illustrating the relationship between the local discrete error $\left|\left|\mu \right|-\left|{\mu }_0\right|\right|$ and the discrete number m. These curves offer a direct comparison of the convergence performance of methods at different orders, facilitating the determination of the optimal interpolating polynomial order.

From the presentation in Figure 2, Figure 3 and Figure 4, the black circles represent the 4H-FDM, the blue solid lines represent the 3N-FDM, and the red stars represent the second-order FDM. It is observed that the 3N-FDM has a significant advantage in convergence speed compared to methods of other orders. All three methods exhibit a gradual decrease in convergence speed as the axial depth of the cut increases, and this change has a particularly notable impact on the initial convergence speed.

In Figure 2, the initial convergence speed of the 3N-FDM is comparable to that of the second-order FDM, both faster than the 4H-FDM. However, when m reaches 40, the convergence speed of the 3N-FDM begins to significantly increase, gradually widening the gap with the second-order FDM. Meanwhile, the convergence speed of the 4H-FDM surpasses the second-order FDM at m = 50 and becomes consistent with that of the 3N-FDM at m = 70.

In Figure 3, the initial convergence speed ranking is second-order FDM, 3N-FDM, and 4H-FDM. As m increases, the 3N-FDM overtakes the second-order FDM at m = 40, while the 4H-FDM matches the second-order FDM at m = 50. When m = 60, the convergence speeds of the 3N-FDM and 4H-FDM are roughly the same.

In Figure 4, the initial convergence speed ranking is second-order FDM, 4H-FDM, and 3N-FDM. However, as m further increases, the 3N-FDM begins to exhibit a faster convergence speed at m = 40, gradually closing the gap with the second-order FDM. Notably, the convergence speeds of the 4H-FDM and the second-order FDM are almost the same at m = 50.

Among the three methods, the second-order FDM exhibits the best stability performance, with its convergence speed less affected by changes in the cutting depth $a_p$. In contrast, although the 3N-FDM experiences fluctuations in the initial convergence speed due to the influence of the cutting depth, it typically achieves the fastest convergence speed at m = 50. The convergence speed of the 4H-FDM, although slower than the 3N-FDM, is relatively less affected by the cutting depth $a_p$. It is concluded that using higher-order polynomial interpolation methods results in a greater computational workload, leading to slower convergence speeds and requiring more computational time.

To further investigate the performance of these methods, the proposed improved fully discrete method was employed to calculate the milling stability lobe diagram. The exact stability limit diagram (SLD) was computed using the second-order FDM with a discrete number m = 200 serving as the reference curve. Subsequently, the differences between the 4H-FDM, 3N-FDM, and first-order SDM at various levels of discretization (m = 10, m = 20, m = 30, and m = 40) were compared with those of the second-order FDM. Specifically, the lines in the comparison are as follows: the blue solid line for the second-order FDM; red dashed line for the first-order SDM; green dashed line for the 3N-FDM; and black dashed line for the 4H-FDM. The relevant results are shown in Figure 5, Figure 6, Figure 7 and Figure 8.

As shown in Figure 5, the discretization number m is 10, and the calculation time for the calculation formula is 34.1 s. The lobe diagram of the 3N-FDM shows the highest similarity and accuracy compared to the second-order FDM. Within the speed range of 5000 r/min to 10,000 r/min, the 3N-FDM performs optimally; within the speed range of 10,000 r/min to 12,000 r/min, the 3N-FDM exhibits a certain error; within the speed range of 12,000 r/min to 13,000 r/min, the calculation accuracy of the 3N-FDM is better than that of the 4H-FDM; and within the speed range of 13,000 r/min to 15,000 r/min, the calculation accuracy of the 4H-FDM is better than that of the 3N-FDM, and subsequently, the calculation accuracies of the 3N-FDM and 4H-FDM are consistent. It is concluded that at a discrete number m = 10, the calculation accuracy of the 3N-FDM is superior to that of the 4H-FDM. The first-order SDM, due to its smaller discrete value, remains largely unchanged, resulting in the lowest accuracy.

As shown in Figure 6, the discretization number m is 20, the calculation time for the calculation formula is 53.8 s, and the 3N-FDM still maintains the highest accuracy. Within the speed range of 10,000 r/min to 15,000 r/min, the accuracy of the 4H-FDM is lower than that of the first-order SDM; in other speed ranges, the two perform basically the same.

As shown in Figure 7, the discretization number m is 30, and the calculation time for the calculation formula is 81.2 s. The lobe diagrams of the 3N-FDM, 4H-FDM, and first-order SDM almost completely overlap with the second-order FDM. Within the speed range of 11,000 r/min to 12,000 r/min, the computational accuracy of the three methods is basically consistent, but they are slightly lower than that of the second-order FDM, although they also reach an ideal level of accuracy.

As shown in Figure 8, the discretization number m is 40, and the calculation time is 128.5 s. The accuracy of the three methods is further improved, and their lobe diagrams almost completely overlap with the second-order FDM, all reaching a very ideal level of accuracy.

Through a comprehensive comparison of the convergence speed and the two-degrees-of-freedom milling stability lobe diagram, among these three methods, the 3N-FDM demonstrates the best performance in approximating the system state item, followed by the 4H-FDM, and the first-order SDM ranks last. It is worth noting that although the 4H-FDM uses a higher interpolation order, its computational accuracy is lower than that of the 3N-FDM. This is mainly because computational errors accumulate during the approximation process of the system state item, leading to a decrease in computational accuracy. Therefore, when choosing an interpolation method, it is not always the case that a higher order is better; instead, multiple factors such as convergence speed, computational accuracy, and computational cost need to be considered comprehensively.

As the discrete number m increases, the computational accuracy does indeed improve, but at the same time, the computation time also increases significantly. However, when comparing the computational accuracy at a discrete number m = 200 with m = 40, it is found that the two are basically consistent. This indicates that a larger discrete number m is not always better because an excessively large value of m can lead to a significant decrease in computational efficiency, with only a slight improvement in accuracy.

Based on the current research situation presented in this article, a reasonable range for the discrete number m can be roughly determined, which is between 35 and 45. Of course, this range is not absolute and may vary depending on the specific research context. In practical applications, an appropriate value of m can be selected based on the precision requirements and computational resources to ensure computational accuracy while maximizing computational efficiency.

In summary, when enhancing the computational accuracy of the fully discrete method, it is not sufficient to merely rely on increasing the interpolation order while neglecting the importance of convergence speed. Research on optimizing the fully discrete method should aim to achieve the dual objectives of high computational accuracy and rapid convergence. Furthermore, during the optimization process of the fully discrete method, many parameters are still set based on experience. The scientific calculation of these parameters will be the key to further improving the computational accuracy of the fully discrete method in the future.

2.3. Verification Test of Milling Stability Using the 3N-FDM

By applying the high-order interpolating polynomial method to approximate the system state item, it is concluded that the 3N-FDM has the optimal performance in approximating the system state item, followed by the 4H-FDM. Using the 3N-FDM fully discrete method and the system parameters listed in Table 1, along with a specific radial cutting depth rate $a_p/D=1$, the lobe diagram is calculated. In this lobe diagram, the upper region of the curve represents the chatter region, while the lower region represents the stable cutting region.

By analyzing the waveform characteristics and the magnitude of the cutting force signals, one can ascertain whether chatter vibration occurs during the machining process. When chatter vibration emerges in the cutting process, the waveform characteristics of the cutting force signals undergo noticeable changes. Chatter vibration results in an augmentation of the amplitude of the cutting force curve, and concurrently, its frequency also rises, causing greater fluctuations in the cutting force signal. Furthermore, chatter vibration gives rise to a substantial increase in the cutting force. This is because chatter vibration, being a detrimental vibration between the tool and the workpiece during the cutting process, alters the contact area and contact pressure between them, thereby inducing an elevation in the cutting force.

To verify the accuracy of the lobe diagram calculated by the 3N-FDM fully discrete method, an experimental validation approach is adopted. Given the restrictions on the spindle speed $n$ due to the experimental conditions, priority is given to testing at lower rotational speeds. To this end, a 3 × 2 test point array containing 6 test points is designed in the left region of the lobe diagram. These test points are named according to the row–column naming convention, sequentially numbered as point 1-1, point 1-2, etc., with the specific test point layout shown in Figure 9.

Cutting experiments were conducted for each of these test points. After the experiments, the occurrence of chatter during the machining process was determined by observing the changes and amplitude fluctuations in the cutting force curves.

Based on the information presented in Figure 9, corresponding cutting depths and rotational speeds have been clearly defined for each test point. After establishing these fundamental processing parameters, further refinement of the feed rate $f$ was achieved through exhaustive research and in-depth consultations with engineers specializing in the cutting of the GH4169 superalloy. Consequently, a detailed processing parameter table was formulated, as outlined in

Table 2. Specific processing parameters for the test experiment.

Number	$\mathbf{Cutting} \mathbf{Depth} {\mathit{a}}_{\mathit{p}}\left(\mathbf{m}\mathbf{m}\right)$	$\mathbf{Spindle} \mathbf{Speed} \mathit{n}(\mathbf{r}/\mathbf{m}\mathbf{i}\mathbf{n})$	$\mathbf{Feed} \mathbf{Rate} \mathit{f}(\mathbf{m}\mathbf{m}/\mathbf{m}\mathbf{i}\mathbf{n})$
1-1	1	7000	70
1-2	1	8250	80
2-1	0.55	7000	110
2-2	0.55	8250	120
3-1	0.2	7000	150
3-2	0.2	8250	160

.

During the cutting experiments, the GH4169 superalloy was selected as the experimental material. The GH4169 superalloy is a typical difficult-to-machine material, prone to chatter during processing. Conducting experiments with this material can more accurately reflect chatter conditions. The specific material properties of the GH4169 superalloy are shown in

Table 3. Performance parameters of GH4169 high-temperature alloy.

Chemical Composition	Ratio	Mechanical Properties	Parameters
Ni	50~55%	Tensile Strength σb (MPa)	1030–1100 (20–25 °C)
Cr	17~21%	Yield Strength σs (MPa)	850–1000 (20–25 °C)
Mo	2.8~3.3%	Rockwell Hardness (HRC)	35–45 (20–25 °C)
Nb	4.75~5.5%	Elongation after Fracture	30% (20–25 °C)
Ti	0.65~1.15%	Fatigue Limit (MPa)	400 (<600 °C)

below.

To accurately measure the cutting forces during the milling process, a triaxial piezoelectric dynamometer model 5073 produced by the Swiss KISTLER company was used, along with a dynamometer stand model 9257 A from KISTLER. To ensure the real-time acquisition and accuracy of the data, an HBM GEN5i handheld high-speed data acquisition system was used, with the sampling frequency set at 6000 Hz. The cutting operation was performed on a DMG103 machine tool, using a φ16R3 mm ring cutter. The final form of the experimental platform is shown in Figure 10.

The cutting tests were conducted based on the machining parameters listed in Table 2. In each set of cutting tests, the consistency of the machining environment was ensured, and the cutting area remained the same. After the cutting tests, the cutting force curves under different cutting parameters were obtained, and the magnitudes of the cutting forces were statistically analyzed.

The specific cutting force curves have been systematically summarized and presented in

Table 4. The cutting force after the completion of the cutting at each point.

	First Column	Second Column
First Line
Second Line
Third Line

. Based on the precise locations of the test points indicated in Figure 9, the cutting force curves are arranged in a sequential and orderly manner. In each accompanying figure, the cutting force variation curves in the X, Y, and Z directions over time or during the process are depicted sequentially from top to bottom. The average cutting forces, obtained through measurement and calculation, have been meticulously recorded in

Table 5. Cutting forces at each point position.

Number	$\mathbf{Cutting} \mathbf{Depth}$ ${\mathit{a}}_{\mathit{p}}\left(\mathbf{m}\mathbf{m}\right)$	$\mathbf{Spindle} \mathbf{Speed}$ $\mathit{n}(\mathbf{r}/\mathbf{m}\mathbf{i}\mathbf{n})$	$\mathbf{Feed} \mathbf{Rate}$ $\mathit{f}(\mathbf{m}\mathbf{m}/\mathbf{m}\mathbf{i}\mathbf{n})$	$\mathbf{X}\text{-}\mathbf{Direction} \mathbf{Cutting} \mathbf{Force}$ ${\mathit{F}}_{\mathit{x}}\left(\mathbf{N}\right)$	$\mathbf{Y}\text{-}\mathbf{Direction} \mathbf{Cutting} \mathbf{Force}$ ${\mathit{F}}_{\mathit{z}}\left(\mathbf{N}\right)$	$\mathbf{Z}\text{-}\mathbf{Direction} \mathbf{Cutting} \mathbf{Force}$ ${\mathit{F}}_{\mathit{z}}\left(\mathbf{N}\right)$
1-1	1	7000	70	615.3	601.8	1713.6
1-2	1	8250	80	690.4	674.2	2085.2
2-1	0.55	7000	110	413.9	386.3	1413.3
2-2	0.55	8250	120	537.8	467.4	1741.1
3-1	0.2	7000	150	280.9	352.1	1257.3
3-2	0.2	8250	160	229.8	235.7	1086.1

. The presentation of these data adheres to the standardized format of Table 2, ensuring the accuracy and readability of the information.

From the data in Table 4, it can be observed that under the processing parameters corresponding to points 1-2 and 2-2, the cutting force curves exhibit significant fluctuations, accompanied by steep peaks and valleys. These characteristics indicate that during machining under these two sets of parameters, obvious tool wear occurred, along with the occurrence of chatter. In contrast, the cutting force curves under other processing parameters show smaller fluctuations, indicating a relatively stable cutting process in a stable cutting state.

Based on the cutting force statistics in Table 5, a cutting force variation curve graph was drawn, as shown specifically in Figure 11. It is found that at points 1-2 and 2-2, abnormal increases in the cutting forces in all three directions occurred. This phenomenon indicates that chatter occurred during machining at these two test points. Chatter not only exacerbates tool wear but also generates greater friction, hindering the smooth progress of the cutting process and leading to a significant increase in cutting forces.

Based on the in-depth analysis of the cutting force curves and their numerical values, it is determined that chatter occurred at points 1-2 and 2-2, while the remaining test points were in a stable cutting state. This conclusion aligns with the stability boundary predicted by the 3N-FDM fully discrete method, thus verifying the accuracy of the 3N-FDM fully discrete method. Therefore, it can be confidently stated that this method is reliable in practical applications and can be directly applied to subsequent research work. The specific judgment results are clearly presented in Figure 12 through graphical representation.

3. Impact of System State Item on Processing Parameters

To further explore the impact of processing parameters on stability and to uncover the underlying relationship between processing parameters and the system state item, research on tool wear and workpiece surface roughness was conducted. During this process, not only were the experimental results tallied, but a multifactorial variance analysis was also performed to ensure the accuracy and reliability of the conclusions.

Subsequently, the lobe diagram calculated using the 3N-FDM fully discrete method was used as a guide to continue the experimental work. Based on the original Figure 9, the number of experimental points was further increased, resulting in a 4 × 4 array that contained a total of 16 test points, as shown in Figure 13. After conducting an extensive literature review and in-depth discussions with engineers specializing in cutting the GH4169 superalloy, specific machining parameters were formulated and are detailed in

Table 6. Specific processing parameters for continuing the experiment.

Number	$\mathbf{Cutting} \mathbf{Depth}$ ${\mathit{a}}_{\mathit{p}}\left(\mathbf{m}\mathbf{m}\right)$	$\mathbf{Spindle} \mathbf{Speed}$ $\mathit{n}(\mathbf{r}/\mathbf{m}\mathbf{i}\mathbf{n})$	$\mathbf{Feed} \mathbf{Rate}$ $\mathit{f}(\mathbf{m}\mathbf{m}/\mathbf{m}\mathbf{i}\mathbf{n})$
1-1	1.4	7000	50
1-2	1.4	8250	55
1-3	1.4	9250	60
1-4	1.4	10,250	65
2-1	1	7000	70
2-2	1	8250	80
2-3	1	9250	90
2-4	1	10,250	100
3-1	0.55	7000	110
3-2	0.55	8250	120
3-3	0.55	9250	130
3-4	0.55	10,250	140
4-1	0.2	7000	150
4-2	0.2	8250	160
4-3	0.2	9250	180
4-4	0.2	10,250	200

.

Cutting tests were conducted based on the machining parameters listed in Table 6 (note: previously completed experiments were not repeated). In each set of cutting tests, the consistency of the machining environment was ensured, and the cutting area remained the same. After the cutting operations, the wear condition of the tool and the roughness of the workpiece surface were inspected.

For a precise evaluation of the tool wear conditions, the German-made ZOLLER genius 3 fully automatic tool inspection machine is employed. This equipment boasts impressive specifications, including an X-axis travel of 600 mm, Y-axis travel of 175 mm, and Z-axis travel of 100 mm, with a repeatability of ±2 μm, a display and positioning accuracy of 1 μm, a concentricity of 2 μm, and an incident light magnification of up to 200 times. During use, the tool to be measured is placed securely in the measuring fixture to ensure proper fixation and alignment. Then, through the computer interface integrated with the tool measuring machine, the appropriate measurement program is selected to execute the measurement, with results displayed instantly. The appearance of this equipment is shown in Figure 14.

To quantify the surface roughness of workpieces, the Chinese-made SMIP-200 surface roughness tester is utilized. This tester features a high resolution of 0.01 μm and a measurement accuracy of (3 + L/200) μm, with a sampling length of 480 mm and an evaluation length of 680 mm. During the measurement process, the stylus needs to gently contact the measured surface and maintain good contact; then, the tester is activated, and the stylus glides linearly at a constant speed across the measured surface while accurately collecting surface roughness data. Upon completion of the measurement, the tester automatically processes the data and displays the surface roughness parameters. The appearance of this tester is shown in Figure 15.

Based on the processing parameters provided in Table 6, the cutting experiments were conducted. During the experimental process, the wear condition of each cutting tool was measured, and a detailed measurement and statistical analysis were carried out on the surface roughness of each workpiece. The summary of the tool wear conditions is presented in

Table 7. Tool wear after cutting at each point position.

	1	2	3	4
1
2
3
4

, while the measurement results of the workpiece surface roughness are recorded separately in

Table 8. Workpiece surfaces after cutting at each point position.

	1	2	3	4
1
2
3
4

and

Table 9. Surface roughness measurement results.

Number	$\mathbf{Cutting} \mathbf{Depth}$ ${\mathit{a}}_{\mathit{p}}\left(\mathbf{m}\mathbf{m}\right)$	$\mathbf{Spindle} \mathbf{Speed}$ $\mathit{n}(\mathbf{r}/\mathbf{m}\mathbf{i}\mathbf{n})$	$\mathbf{Feed} \mathbf{Rate}$ $\mathit{f}(\mathbf{m}\mathbf{m}/\mathbf{m}\mathbf{i}\mathbf{n})$	$\mathbf{Ra} \left(\mathsf{\mu }\mathbf{m}\right)$
1-1	1.4	7000	50	6.4
1-2	1.4	8250	55	7.1
1-3	1.4	9250	60	5.6
1-4	1.4	10,250	65	8.3
2-1	1	7000	70	4.8
2-2	1	8250	80	7.6
2-3	1	9250	90	4.3
2-4	1	10,250	100	6.7
3-1	0.55	7000	110	4.1
3-2	0.55	8250	120	6.2
3-3	0.55	9250	130	4.4
3-4	0.55	10,250	140	6.3
4-1	0.2	7000	150	3.4
4-2	0.2	8250	160	3.2
4-3	0.2	9250	180	3.8
4-4	0.2	10,250	200	3.5

.

After analyzing the inspection results in Table 7, the following phenomena are observed:

At test points 1-1, 1-2, 1-4, 2-2, 2-4, 3-2, and 3-4, severe tool wear occurred, and in some cases, chipping of the tool occurred. Especially at points 1-4 and 2-4, the chipping of the tool was particularly evident, resulting in a significant decrease in the surface quality of the workpiece processed under these point parameters.

At points 1-3, 2-1, and 2-3, the wear on the tool was relatively light, mainly manifested as general localized wear or partial chipping of the cutting edge, and the surface quality of the workpiece processed under these point parameters was at an average level.

At points 3-1, 3-3, 4-1, 4-2, 4-3, and 4-4, the wear on the tool was the least severe, with only localized wear on the cutting edge, and the surface quality of the workpiece processed under these point parameters reached an excellent level.

The surface condition of the workpieces has been summarized in the following Table 8, while the specific measurement results are detailed in Table 9. Based on the data from Table 9, a surface roughness variation curve has been plotted and is shown in Figure 16. These data and observations provide an important basis for the further analysis and optimization of the cutting parameters.

After an in-depth analysis of the workpiece surface conditions presented in Table 8 and Table 9, the following conclusions are drawn:

At test points such as 1-1, 1-2, 1-4, and 2-2, the surface roughness exceeded 6.3 μm, categorizing them as rough surfaces. These workpiece surfaces exhibited prominent “fish scale marks,” which were caused by severe wear of the cutting tool at these locations. This led to increased friction between the workpiece and the tool, leaving obvious traces after chatter occurred. The emergence of these marks significantly increased the surface roughness of the workpiece.

At points such as 2-4, 3-2, and 3-4, the surface roughness was predominantly above 6.3 μm, indicating rough surfaces with clearly visible and unevenly distributed cutting marks. These were accompanied by variations in the cutting depths, suggesting that chatter occurred during the cutting process at these points and that the wear of the cutting tool was also severe.

At points such as 1-3, 2-1, and 2-3, the surface roughness ranged between 3.2 μm and 6.3 μm, categorizing them as moderately rough surfaces. Although the cutting marks were evident, they were relatively evenly distributed, indicating that the cutting process was relatively stable at these points and that the wear of the cutting tool was also relatively moderate.

At points such as 3-1, 3-3, 4-1, 4-2, 4-3, and 4-4, the surface roughness was approximately 3.2 μm, categorizing them as medium rough surfaces. The cutting marks were almost imperceptible and extremely densely distributed, further proving that the cutting process was very stable at these points and that the degree of tool wear was relatively light.

In Figure 16, the X-axis represents the experimental sequence, and the Y-axis represents the surface roughness. The graph clearly shows an overall downward trend in surface roughness. Combining the specific parameters in Figure 13 and Table 6, it is concluded that the processing stability gradually improved as the cutting depth $a_p$ decreased.

The experimental results regarding the degree of tool wear and workpiece surface roughness were consistent. Chatter occurred during processing at points such as 1-1, 1-2, 1-4, 2-2, 2-4, 3-2, and 3-4; mild chatter occurred at points such as 1-3, 2-1, and 2-3; while stable processing was achieved at points such as 3-1, 3-3, 4-1, 4-2, 4-3, and 4-4. This confirms the accuracy of the inferences about chatter at each point, as shown in the specific statistical results in Figure 17. The processing conditions at each point coincided with the predicted stability boundaries of the lobe diagram, thus verifying the accuracy of the cutting experiments. The cutting experiment data can be directly applied in subsequent research.

To delve into the combined influence of multiple parameters on surface roughness, based on the statistical results in Table 9, three-dimensional surface graphs illustrating the impact of various parameters on surface roughness were meticulously drawn. These graphs not only intuitively display the coupling relationships between the parameters but also make the analysis results clearer and easier to understand through sequential labeling.

In this collection, each three-dimensional surface graph corresponds to a specific combination of parameters. By observing and comparing these graphs, the degree and trend of the influence of different parameters on surface roughness can be determined.

In Figure 18, the X-axis represents the cutting depth $a_p$, the Y-axis represents the spindle speed $n$, and Z represents the surface roughness. It can be clearly seen from the graph that the cutting depth $a_p$ has a significant impact on surface roughness: as the cutting depth $a_p$ decreases, the surface roughness also shows a downward trend. In contrast, the impact of the spindle speed $n$ on surface roughness is less significant; when the spindle speed $n$ parameter changes, the change in surface roughness is relatively small.

In Figure 19, the X-axis still represents the cutting depth $a_p$, the Y-axis becomes the feed rate $f$, and Z represents the surface roughness. In this graph, both the cutting depth $a_p$ and feed rate $f$ have a significant impact on surface roughness: a decrease in cutting depth $a_p$ leads to a reduction in surface roughness, while an increase in feed rate $f$ reduces surface roughness. This indicates that by appropriately reducing the cutting depth $a_p$ and increasing the feed rate $f$, a lower surface roughness can be achieved.

In Figure 20, the X-axis represents the spindle speed $n$, the Y-axis represents the feed rate $f$, and Z represents the surface roughness. It can be seen from the graph that the feed rate $f$ has a more significant impact on surface roughness compared to rotational speed: as the spindle speed $n$ and feed rate $f$ increase, the surface roughness shows a downward trend.

To more accurately determine the significance of the impact of each parameter on the surface roughness, based on the statistical results in Table 9, the single-factor variance analysis method was selected. The variance analysis was performed using SPSS 27 software, and the specific analysis results are shown in Figure 21, and the statistical results are shown in

Table 10. Statistical results of analysis of variance.

	F	${\mathbf{eta}}^2$	p	Significance
Cutting Depth $a_p$	3.667	0.423	0.044	Significant
Spindle Speed $n$	1.94	0.280	0.052	Insignificant
Feed Rate $f$	7.085	0.673	0.011	Significant

.

The analysis of the statistical results in Table 10 led to the following conclusions:

The F-value for the cutting depth $a_p$ is 3.667, the p-value is 0.044, and the ${\mathrm{eta}}^2$-value is 0.423. Since the p-value is less than the significance level of 0.05, it is considered that the cutting depth $a_p$ has a significant impact on surface roughness. Additionally, due to the effect size being greater than 0.3, it is believed that this impact has a substantial effect.

The F-value for the spindle speed $n$ is 1.94, the p-value is 0.052, and the ${\mathrm{eta}}^2$-value is 0.28. Since the p-value is greater than the significance level of 0.05, it is considered that the spindle speed $n$ does not have a significant impact on surface roughness. However, its effect size is greater than 0.15, indicating that it still has a certain medium effect.

The F-value for the feed rate $f$ is 7.085, the p-value is 0.011, and the ${\mathrm{eta}}^2$-value is 0.673. Since the p-value is less than the significance level of 0.05, it is considered that the feed rate $f$ has a significant impact on surface roughness, and this impact has a very large effect.

In summary, the cutting depth $a_p$ and feed rate $f$ have a significant impact on surface roughness, while the spindle speed $n$ does not have a significant impact. Therefore, it can be inferred that during cutting processing, the cutting depth $a_p$ and feed rate $f$ have a significant impact on cutting stability, while the impact of the spindle speed $n$ is relatively small.

Further comparing the impact of the cutting depth $a_p$ and feed rate $f$, it is found that the p-value and effect size of the cutting depth $a_p$ are both smaller than those of the feed rate $f$. Considering that the system state item (such as the displacement and trajectory of the tool) is a comprehensive indicator reflecting the instantaneous status of the system, it can be inferred that changes in the system state item have the greatest impact on the feed rate $f$, followed by the cutting depth $a_p$, and finally the spindle speed $n$.

Based on these conclusions, it is determined that optimizing the system state item can more efficiently explore the comprehensive impact of processing parameters on processing quality, production efficiency, and tool wear. This discovery not only clarifies the advantages brought by optimizing the system state item but also provides valuable guidance for subsequent research and actual production.

4. Conclusions

This study tightly integrates theory with practice, successfully highlighting the pivotal role of optimizing system state terms in the research of process parameter optimization. The high-order interpolation method was adopted to optimize system state terms, and through comparisons of computational accuracy and convergence speed, it was ultimately found that the 3N-FDM exhibited particularly outstanding performance.

Furthermore, the specific value of the discretization number m was thoroughly discussed. The research results indicated that although increasing the value of m can enhance computational accuracy to a certain extent, it also correspondingly increases computation time. Therefore, in practical applications, it is recommended to set the value of m within the range of 35 to 45. Of course, this range is not absolute and may need to be adjusted appropriately according to specific research circumstances. During the optimization process of the fully discrete method, many parameters are still set based on experience, and the scientific calculation of these parameters will be the key to further improving the computational accuracy of the fully discrete method in the future.

By conducting a detailed analysis of the waveform characteristics and values of cutting force signals, the accuracy of the 3N-FDM was successfully verified. On this basis, further cutting experiments were conducted, and through a one-way ANOVA, the influencing factors of tool wear conditions and surface roughness of the workpiece were thoroughly investigated. The research results showed that the cutting depth $a_p$ and feed rate $f$ have significant impacts on surface roughness, while the influence of the spindle speed $n$ is relatively smaller. At the same time, it was also inferred that changes in system state terms have the most significant impact on the feed rate $f$, followed by the cutting depth $a_p$, and finally the spindle speed $n$.

In summary, optimizing system state terms can more effectively explore the comprehensive impacts of processing parameters on processing quality, production efficiency, and tool wear. This conclusion not only deepens the understanding of parameter optimization in the machining process but also provides valuable guidance and reference for the selection of process parameters in actual production.