Operational Modal Analysis of CNC Machine Tools Based on Flank-Milled Surface Topography and Cepstrum

Abstract

Conducting research on the dynamics of machine tools can prevent chatter during high-speed operation and reduce machine tool vibration, which is of significance in enhancing production efficiency. As one of the commonly used methods for studying dynamic characteristics, operational modal analysis is more closely aligned with the actual working state of mechanical structures compared to experimental modal analysis. Consequently, it has attracted widespread attention in the field of CNC machine tool dynamic characteristics research. However, in the current operational modal analysis of CNC machine tools, discrepancies between the excitation methods and the actual working state, along with unreasonable vibration response signal acquisition, affect the accuracy of modal parameter identification. With the development of specimen-based machine tool performance testing methods, the practice of identifying machine tool characteristics based on machining results has provided a new approach to enhance the accuracy of CNC machine tool operational modal analysis. Existing research has shown that vibration significantly influences surface topography in flank milling. Therefore, a novel operational modal analysis method is proposed for the CNC machine tool based on flank-milled surface topography. First, the actual vibration displacement of the tooltip during flank milling is obtained by extracting vibration signals from surface topography, which enhances the accuracy of machine tool operational modal analysis from both the aspects of the excitation method and signal acquisition. A modified window function based on compensation pulses is proposed based on the quefrency domain characteristics of the vibration signals, which enables accurate extraction of system transfer function components even when the high-frequency periodic excitation of the machine tool causes overlap between the system transfer function components and the excitation components. Experimental results demonstrate that the proposed method can obtain accurate operational modal parameters for CNC machine tools.

1. Introduction

With the development of cutting-edge technologies such as aerospace, the demand for high-precision and complex components in industries is increasing steadily. Consequently, as one piece of key manufacturing equipment for precision components, CNC machine tools are also required to have higher production efficiency. For CNC machine tools, high production efficiency implies high feed rates, high cutting speeds, significant load variations, or frequent substantial accelerations and decelerations. Under such operating conditions, the structural vibration of the machine tool becomes significant [1,2]. Vibration often affects the surface quality of parts and can even lead to structural damage to the machine tool, resulting in safety incidents. Therefore, machine tool vibration is considered a major obstacle limiting further improvement in production efficiency. Dynamic characteristics of machine tools determine the severity of vibration, making it a focal point of research in the field of CNC machine tools [3,4].

Modal analysis is one of the common methods to study dynamic characteristics, with its core focus on identifying modal parameters. Depending on the state of the identification object, modal analysis of CNC machine tools can be classified into experimental Modal Analysis (EMA) and operational modal analysis (OMA). EMA refers to applying artificial excitation to a stationary machine tool and recording both the excitation signal applied as input and the vibration response signal as output using sensors. Based on these recorded signals, frequency response function curves are computed, and ultimately, the modal parameters of the machine tool are obtained by fitting the curves. Since its first proposal in the 1970s, EMA has undergone half a century of development. With its outstanding ability to resolve modal orders and accurate parameter identification, it has become one of the most commonly used modal analysis methods [5,6,7]. However, as a complex mechanical system composed of multiple components, CNC machine tools exhibit differences in dynamic characteristics between operational and static states due to various influencing factors such as contact conditions between components, component temperatures, and component centroid positions [8,9]. Therefore, to obtain CNC machine tool system modes that are more practically relevant, OMA, primarily applied in civil engineering, has been introduced into the study of CNC machine tool dynamic characteristics. OMA of CNC machine tools involves using the excitation of the measured machine tool structure in its working state and identifying the modal parameters of the machine tool based solely on the collected vibration signals of the machine tool system [10,11,12].

However, unlike OMA of civil engineering structures, the excitation in machine tool OMA often includes periodic components with extremely high energy, therefore affecting the effective identification of modal parameters. To address this challenge, scholars have conducted extensive research with the aim of eliminating the influence of periodic components. Current solutions can be categorized into two main approaches: artificially designing excitation conditions and signal processing. However, the approach of artificially designing excitation conditions contradicts the original purpose of OMA, whereas signal processing methods aimed at eliminating periodic components are more valuable in terms of research [13,14,15]. Methods for processing vibration signals can be classified into time-domain methods and frequency-domain methods based on the analytical perspective. However, these classical methods are all premised on the presence of wideband white noise in the system excitation, which is often difficult to meet in practice, especially for machine tools requiring high-speed spindle rotation [16,17,18]. Although many studies have attempted to address this issue, the algorithms often tend to be overly complex and lack generalization capability [19,20,21]. In recent years, cepstrum analysis, originally used for speech signal processing and seismic echo detection, has emerged as a new research hotspot in the field of OMA owing to its superior deconvolution capabilities. In 1996, Gao et al. successfully reconstructed the transfer function of a single input with multiple output (SIMO) system using cepstrum analysis, marking the first implementation of cepstrum-based OMA [22]. Subsequently, Hanson et al. expanded the applicability of cepstrum-based OMA from SIMO systems to multiple input with multiple output (MIMO) systems. They successfully extracted the system transfer function under the influence of periodic excitation [23]. Smith et al. improved the identification efficiency of cepstrum-based OMA and reduced the complexity of the algorithm through optimization of the fitting method [24]. In recent years, with the increasingly widespread application of cepstrum analysis in modal analysis, cepstrum-based OMA has also been gradually refined [25,26]. Currently, cepstrum-based OMA is only applied in fault detection for simple rotating mechanical structures such as bearings and gears, with limited application in complex and large-scale mechanical structures like CNC machine tools [27,28]. Due to the relatively simple structures of the application objects and the clear distinction between periodic excitation and system transfer function distributions in the quefrency domain, existing cepstrum-based OMA can separate different components using only traditional window functions or notch filters. However, for complex mechanical structures like CNC machine tools primarily involved in precision machining tasks, dynamic characteristics of the system are not only more intricate but also subject to high-frequency periodic excitation due to high spindle speed. In such cases, the overlap between periodic excitation and system transfer functions in the quefrency domain is highly likely, therefore affecting the accuracy of modal parameter identification.

Additionally, the vibration signal in OMA is usually collected by placing sensors on the mechanical structure. However, compared to other simple rotating machinery, the environment for collecting vibration signals from operating machine tools is more complex. The working space of the machine tool is filled with splashing cutting fluid, chips, and mist during machining. In such harsh environments, it is not feasible to obtain accurate vibration signals through sensors placed on the machine structure. While it is possible to place sensors on the exterior casing of the machine tool to collect vibration signals, research on machine tool dynamic characteristics often focuses on the scenario where the tooltip serves as the system output point. However, attempts to improve the collection environment by designing the machining process or altering cutting parameters to reduce chips, mist, and heat may distort the working state. With standardized specimens widely used in machine tool performance tests, there has been growing academic interest in identifying machine tool errors and characteristics based on machining results which also provides a new approach to obtaining vibration signals of the machine tool [29,30,31]. Extracting vibration signals from surface topography offers advantages for the OMA of the machine tool. Not only does it eliminate the need for manually adjusting the machine tool state, but the vibration signal collection process is also less affected by environmental factors. As a result, the identification results can be more aligned with the actual working state of the machine tool. It is widely accepted that machine tool vibration does indeed affect the surface topography of workpieces, and adjusting cutting parameters can improve the surface quality of workpieces. Zhang et al. derived a formula for surface topography displacement considering vibration based on the geometric principles of point contact cutting. Through an analysis based on the formula, they concluded that the influence of horizontal vibration on surface morphology is smaller than that of vertical vibration, while both can be reduced by adjusting cutting parameters [32]. Cai et al. utilized end milling surfaces as an example and employed finite element analysis to investigate the relationship between two-dimensional tool paths and three-dimensional surface topography [33]. They established a surface topography model considering system vibration and tool deformation after multiple cutting passes and then illustrated that vibration has a direct impact on the surface topography phase while tool runout errors indirectly affect the surface topography amplitude. While Zhang et al. considered the vibration amplitude and cutting parameters when analyzing the influence mechanism between cutting-edge trajectory and surface topography, the model they established did not account for the coupling effect between vibration and cutting parameters [34]. Currently, there is an amount of research on pairwise relationships, but there is still a lack of systematic research on vibration, cutting parameters and surface topography, resulting in unclear mechanisms among the three. This also often leads to unreasonable cutting parameter optimization for surface quality. Research on extracting vibration signals from surface topography can further elucidate the relationship between vibration and surface topography, thus refining the principles of cutting parameter optimization for surface quality. The optimized principles can help avoid selecting overly conservative feed per tooth to ensure surface quality, which may lead to inefficient machining.

Aimed at obtaining a more realistic representation of the dynamic characteristics of the machine tool under working conditions, this paper proposes a CNC machine tool operational modal analysis method to enhance the identification accuracy of machine tool operational modal parameters. The method involves extracting vibration signals of the machine tool from the surface topography of flank-milled workpieces and obtaining modal parameters through cepstrum analysis. The research systematically considers the interaction relationship between vibration, cutting parameters, and surface topography, providing theoretical guidance for optimizing cutting parameters aimed at surface quality. Furthermore, the research enhances the extraction accuracy of system transfer functions when periodic excitation and system transfer functions overlap in the quefrency domain using a modified window function. A surface topography generation model for flank milling considering vibration was first established in Section 2. In Section 3, the feasibility of extracting vibration signals from surface topography was illustrated, and a specific vibration signal acquisition method based on tooth passing frequency (TPF) was proposed. The necessary conditions for implementing this method were analyzed, and thus, the influence mechanism among vibration, cutting parameters, and surface topography was further illustrated. Section 4 focused on cepstrum-based OMA, emphasizing the characteristics of vibration signals in the quefrency domain during flank milling and the cepstrum editing process based on the characteristics. Finally, in Section 5, a detailed introduction to the experimental validation process and results was provided.

2. Surface Topography Generation Model for Flank Milling

The surface topography generation model proposed in this section is aimed at the flank milling of developable ruled surfaces. The flank milling process of a workpiece with an end mill cutter is illustrated in Figure 1. The vibration of the machine tool is reflected onto the surface topography of the machined workpiece through its influence on the trajectory of the cutting edge. Therefore, the first step is to establish a trajectory model of the cutting edge in flank milling considering vibration. A tool coordinate system (TCS) O_T-X_TY_TZ_T is established at the center of the cutter tip plane. The Z_T axis aligns with the cutter axis, while the X_T and Y_T axes lie on the cutter tip plane and are mutually perpendicular.

Each cutting edge is discretized into $N$ cutting-edge elements according to equal length. For the $j$-th ($j=\mathrm{1,2},3\dots ,N_t$, where $N_t$ is the number of cutting edges of the tool) cutting edge, the $i$-th cutting-edge element $TE\left(i,j\right)$ has its homogeneous coordinates in the TCS as follows:

$$\left[\begin{array}{c}{{}_E{}^{T}x}_{ij}\\ {{}_E{}^{T}y}_{ij}\\ \begin{array}{c}{{}_E{}^{T}z}_{ij}\\ 1\end{array}\end{array}\right]={}_E{}^T\mathbf{T}\cdot \left[\begin{array}{c}\mathit{Rcos}{\theta }_{ij}\\ \mathit{Rsin}{\theta }_{ij}\\ \begin{array}{c}R{\theta }_{ij}/\mathit{tan}\beta \\ 1\end{array}\end{array}\right],{}_E{}^T\mathbf{T}=\left[\begin{array}{ccc}\mathit{cos}{\varphi }_j& -\mathit{sin}{\varphi }_j& \begin{array}{cc}0& 0\end{array}\\ \mathit{sin}{\varphi }_j& \mathit{cos}{\varphi }_j& \begin{array}{cc}0& 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}1\\ 0\end{array}& \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\right]$$

(1)

where ${}_E{}^T\mathbf{T}$ is the homogeneous coordinate transformation matrix of each cutting-edge element to the TCS, $R$ is the radius of the cutter. ${\varphi }_j={\varphi }_0+2\pi \left(j-1\right)/N_t$ represents the distribution angle of $TE\left(1,j\right)$ in the TCS, where ${\varphi }_0$ stand for the initial phase and ${\theta }_{ij}$ is the helix lag angle of $TE\left(1,j\right)$. $\beta$ is the helix angle of the cutter.

During the milling process, the cutter moves at a constant feed rate $V_F$ parallel to the surface to be machined while the tool rotates at a constant speed $\omega$ about its axis. A workpiece coordinate system (WCS) O_w–X_wY_wZ_w is established on the workpiece, where the Z_w axis is parallel to the Z_T axis of the TCS, the X_w axis is perpendicular to the surface of the workpiece, and the Y_w axis is parallel to the feed rate $V_F$. Except for machining thin-walled parts, the dynamic stiffness of the workpiece is often larger, resulting in smaller vibration compared to the tool under cutting forces. Therefore, in flank milling, the main consideration is the tool vibration affecting the surface topography of the machined workpiece. The clamping method for the tool during flank milling ensures enough stiffness in the axial direction. Hence, only vibrations on the tool side of the machine tool along the X_w direction and Y_w direction with possible multiple modes are considered. Thus, the homogeneous coordinates of the cutting-edge element $TE\left(i,j\right)$ in the WCS are:

$$\left[\begin{array}{c}{{}_E{}^{W}x}_{ij}\\ {{}_E{}^{W}y}_{ij}\\ \begin{array}{c}{{}_E{}^{W}z}_{ij}\\ 1\end{array}\end{array}\right]={}_T{}^W\mathbf{T}\cdot {}_E{}^T\mathbf{T}\cdot \left[\begin{array}{c}\mathit{Rcos}{\theta }_{ij}\\ \mathit{Rsin}{\theta }_{ij}\\ \begin{array}{c}R{\theta }_{ij}/\mathit{tan}\beta \\ 1\end{array}\end{array}\right],{}_T{}^W\mathbf{T}=\left[\begin{array}{ccc}\mathit{cos}\omega t& -\mathit{sin}\omega t& \begin{array}{cc}0& V_x\left(t\right)\end{array}\\ \mathit{sin}\omega t& \mathit{cos}\omega t& \begin{array}{cc}0& V_{F}t\end{array}{+V}_y\left(t\right)\\ \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}1\\ 0\end{array}& \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\right]$$

(2)

where ${}_T{}^W\mathbf{T}$ represents the homogeneous coordinate transformation matrix of the cutting-edge element from TCS to WCS. $V_x\left(t\right)$ and $V_y\left(t\right)$, respectively, denote the vibration displacements of the tool along the X_w and Y_w directions at time $t$.

Similar to the discrete modeling of tool cutting-edge trajectory described above, the surface topography of the workpiece can also be viewed as composed of $N$ layers of different height two-dimensional profiles. The surface topography profile of each layer is formed by the cutting-edge element located at the same height. When the tool structure, cutting parameters, and vibration displacement in Equation (2) are specified, the trajectory of each cutting-edge element can be determined according to this equation. Ultimately, the obtained lowest envelope line is regarded as the two-dimensional surface topography of the workpiece at that height, as illustrated in Figure 2. It is worth noting that the figure magnifies a portion of the actual topography to illustrate the principle of surface topography generation clearly.

3. Method for Collecting Vibration Signals from Surface Topography

The mechanism of machine tool vibration on flank-milled surface topography indicates that the milled surface indeed contains vibration displacement perpendicular to the feed direction. In the discretization modeling method adopted in this paper, the influence of machine tool vibration on the trajectory of cutting-edge elements at different heights can be considered the same. Therefore, studying the influence of vibration signals on a cutting-edge element at any height as an example will lead to conclusions that are equally applicable to two-dimensional surface topography profiles at different heights.

The two-dimensional surface topography profile of the milled workpiece can be considered to be a series of ideal structural elements connected in sequence with $\stackrel{⏜}{P_1{P}_2}$ in Figure 3 being one such element. The projection length of $\stackrel{⏜}{P_1{P}_2}$ along the feed direction is equal to the feed per tooth $F_t$. P₁ represents the entry point of the cutting-edge element $TE\left(i,j\right)$ to machining the current ideal structural element, while P₂ represents the entry point of $TE\left(i,j+1\right)$ to machining the latter ideal structural element. Both P₁ and P₂ have the same X_w coordinate in the WCS. O₁ represents the center point position of the tooltip when $TE\left(i,j\right)$ reaches the entry point P₁ at time $t_1$. O₂ represents the center point position of the tooltip when $TE\left(i,j\right)$ reaches the lowest point of the profile contour at time $t_2$. O₃ represents the center point position of the tooltip when $TE\left(i,j\right)$ reaches the entry point P₂ at time $t_3$. O₄ represents the center point position of $TE\left(i,j+1\right)$ when it reaches the entry point P₂. α denotes the angle of rotation of $TE\left(i,j\right)$ between times $t_2$ and $t_3$.

According to the principle of flank milling, the following relationships can be obtained:

$$\begin{gathered} t_3-t_1=2\left(t_3-t_2\right) \\ {t}_3-t_2=\frac{\alpha }{\omega } \end{gathered}$$

(3)

Moreover, there are the following geometric relationships:

$$R\cdot \mathit{sin}\alpha +V_F\cdot \left(t_3-t_2\right)=\frac{1}{2}{F}_t$$

(4)

Let $T_0=t_3-t_1$ and denoting $T_0$ as the effective cutting time per tooth. Substituting Equation (3) and $F_t=2\pi V_F/N_t\omega$ into Equation (4):

$$R\cdot \mathit{sin}\frac{1}{2}\omega T_0+\frac{1}{2}{V}_F{T}_0=\frac{\pi V_F}{\omega N_t}$$

(5)

For the function $f\left(\alpha \right)=\alpha -\mathit{sin}\alpha$, it holds true that $\forall \alpha \in \left(0,+\infty \right),f^{\prime }\left(\alpha \right)\ge 0$. Moreover, in the actual machining process of precision workpieces, $\alpha$ usually does not exceed $\pi /2$. Therefore, there are:

$$\forall T_0\in \left(0,\frac{\pi }{\omega }\right),\exists C\in \left(0,\frac{\pi }{2}-1\right), \frac{1}{2}\omega R{T}_0+\frac{1}{2}{V}_F{T}_0-C=\frac{\pi V_F}{\omega N_t}$$

(6)

The calculation formula for $T_0$ can be obtained when the relevant cutting parameters are given:

$$T_0=\frac{2\pi }{N_t}\cdot \frac{1}{\frac{R}{V_F}\cdot {\omega }^2+\omega }+\frac{2C}{R\omega +V_F}$$

(7)

As shown in Figure 4, According to the principle of signal acquisition, there are theoretically two possible ways to retain vibration signals on the surface topography:

• Effective cutting time per tooth T₀ is longer than the period of any component contained in the vibration signal. This ensures that the vibration signal is completely preserved on only one ideal structural element of the surface topography. In this paper, this sampling method is also referred to as the vibration signal sampling method based on the cutting period per tooth.
• Effective cutting time per tooth T₀ is much shorter than the period of any component contained in the vibration signal. Under this condition, the ideal structural element of the surface topography formed by each tooth element is extremely short compared to the length of the vibration signal. Thus, it can be regarded as a single sampling of the vibration signal. In this paper, this sampling method is also referred to as the vibration signal sampling method based on the TPF.

The above two sampling methods have different requirements for cutting parameters. For vibration signal sampling based on the cutting period per tooth, the necessary condition is that for $T_0\in \left(0,\pi /\omega \right)$, it holds true for $\forall C\in \left(0,\pi /2-1\right)$ that:

$$T_0\ge \mathit{max}\left\{T\right\}$$

(8)

where $\left\{T\right\}$ represents the period set of vibration signal components. To ensure that Equation (8) holds true for any $T_0\in \left(0,\pi /\omega \right)$, it is equivalent to:

$$\frac{2\pi }{N_t}\cdot \frac{1}{\frac{R}{V_F}\cdot {\omega }^2+\omega }\ge \mathit{max}\left\{T\right\}$$

(9)

Equation (9) indicates that vibration signals within any frequency band can be theoretically accommodated by appropriately designing cutting parameters, tool tooth number, and tool radius, ensuring their complete representation in the ideal structural element of the surface topography. However, in practical machining scenarios, machine tool vibrations often occur within the low-frequency range. According to Equation (9), this situation necessitates minimizing the spindle speed $\omega$, maximizing the feed rate $V_F$, and minimizing the machine tool radius. The milling process under such special conditions has significant discrepancies from the general milling process. For most conventional machine tools, ensuring normal milling under these specific conditions requires machine tool modifications and custom tooling, which are unavoidable. Therefore, vibration signal sampling based on the cutting period per tooth is generally impractical and lacks practical value in typical scenarios, contradicting the research objectives of this paper.

Compared to vibration signal sampling based on the cutting period per tooth, vibration signal sampling based on the cutting frequency per tooth imposes much looser requirements on cutting and tool parameters. The essential condition for retaining vibration signals on the surface profile through this method is that:

$$T_0\ll \mathit{min}\left\{T\right\},f_t\ge 2\mathit{max}\left\{f\right\}$$

(10)

where $f_t=N_t\omega /2\pi$ represents TPF, and $\left\{f\right\}$ denotes the frequency set of vibration signal components. During the vibration signal sampling process, the surface topography of the workpiece is measured and collected with multiple sampling points distributed along each ideal structural element. In other words, the actual sampling process involves multiple samplings of the vibration signal within every effective cutting time per tooth $T_0$. Therefore, $T_0\ll \mathit{min}\left\{T\right\}$ is no longer necessary. It is only essential to ensure that $f_t$ is sufficiently large to prevent distortion of the vibration signal. According to the sampling theorem, taking a frequency proportion factor of 5, i.e., $f_t$ is at least five times that of the highest frequency of the vibration signal components. Thus, the revised necessary condition for the validity of the vibration signal sampling method based on TPF can be expressed as:

$$\frac{N_t\omega }{10\pi }\ge \mathit{max}\left\{f\right\}$$

(11)

Reference [34] indicates that vibrations perpendicular to the feed direction may cause certain structural elements of the surface topography that should be reflected on the surface profile to be masked by other structural elements, resulting in a reduction in the actual number of vibration signal samplings based on TPF. This ultimately affects the integrity of the obtained vibration signal. Figure 5, taking a milling cutter with two teeth as an example, illustrates the incomplete retention of signals on the surface profile due to vibrations perpendicular to the feed direction when the cutting parameters are unreasonable. The blue lines represent the trajectory of $TE\left(i,1\right)$, while the red lines represent the trajectory of $TE\left(i,2\right)$, and the black lines represent the final formed surface profile. The trajectory of $TE\left(i,j\right)$ within the N-th spindle rotation cycle is defined as the N-th cutting curve of $TE\left(i,j\right)$. The solid lines of the trajectory represent the cutting curves corresponding to the final formed surface profile, while the dashed lines represent the cutting curves that did not contribute to the formation of the final surface profile. When $F_t$ is small, due to the influence of vibration, the N + 1-th cutting curve of $TE\left(i,1\right)$, the N-th and N+1-th cutting curve of $TE\left(i,2\right)$ do not contribute to the formation of the final surface profile. Consequently, the corresponding vibration when the cutting-edge elements passing through these curves are not retained on the final surface profile. For the vibration signal sampling method based on the TPF, this implies that the sampling frequency will be lower than TPF, therefore affecting the integrity of the retained vibration signals in the surface profile.

To achieve a sampling frequency equal to the TPF, it is necessary to adjust the cutting parameters according to the vibration conditions, ensuring that milling is performed at a reasonable $F_t$. Figure 6 illustrates the critical state, where the N-th cutting curve of $TE\left(i,j\right)$ intersects with the N-th cutting curve of $TE\left(i,j+1\right)$ at point P on the N-th cutting curve of $TE\left(i,j+1\right)$, which is the lowest point in height from the workpiece surface. Additionally, the vibration displacement of the tool reaches a maximum when $TE\left(i,j+1\right)$ passes through point P. In this state, the distance between O₂ and O₄ in the feed direction is $d$. Only when $F_t$ is greater than $d$ can the N-th cutting curve of $TE\left(i,j+1\right)$ remain unaffected by the adjacent two cutting curves and contribute to forming the final surface profile. Geometrically, this can be expressed as:

$$\begin{gathered} d=V_F\left(t_3-t_2\right)+\sqrt{R^2-(R+h_1-h_2{)}^2} \\ {h}_1=V_x\left(t_3\right)-V_x\left(t_2\right) \\ {h}_2=V_x\left(t_4\right)-V_x\left(t_2\right)={\int }_{t_2}^{\frac{1}{f_t}+t_2}{V_x}^{\prime }\left(t\right)dt \end{gathered}$$

(12)

where $t_2,t_3,t_4$ correspond to the moments when the tooltip center is located at O₂, O₃, and O₄, respectively. $V_x\left(t_2\right), V_x\left(t_3\right),V_x\left(t_4\right)$ represent the vibration displacements of the tooltip center along the X_w direction at these respective moments. Equation (12) demonstrates that the time-varying vibration displacement results in a continuously changing $d$ between every two cutting curves. Thus, only when $F_t$ exceeds the maximum value of $d$, can any cutting curve remain unaffected by the adjacent two cutting curves and effectively contribute to forming the final surface profile. Geometric analysis reveals that $d$ attains its maximum value when $h_2$ reaches its maximum, and the magnitude of $h_2$ is solely determined by $t_2$ and the temporal variation pattern of the vibration signal.

However, characteristics of vibration signals, such as ${V_x}^{\prime }\left(t\right)$ are unknown in practice until the operational modal analysis of the machine tool is completed. Therefore, adjusting $F_t$ based on the temporal variation pattern of the vibration signal before modal analysis is unattainable. However, $\mathit{max}[V_x(t)]-\mathit{min}[V_x(t)]$ can be estimated before modal analysis by measuring the maximum peak-to-valley value on the milled workpiece surface. Additionally, for any $t_2$, the following inequality always holds:

$${\int }_{t_2}^{\frac{1}{f_t}+t_2}{V_x}^{\prime }\left(t\right)dt\le \mathit{max}[V_x(t)]-\mathit{min}[V_x(t)]$$

(13)

The range of $F_t$ can be determined based on Equation (13) with an unknown temporal variation pattern of the vibration signal, therefore ensuring the integrity of the vibration signal retained in the surface topography. Let $Z=\mathit{max}[V_x(t)]-\mathit{min}[V_x(t)]$ and substituting Equation (13) into Equation (12), the range of $F_t$ can be obtained as:

$$F_t\ge V_F\left(t_3-t_2\right)+\sqrt{R^2-(R+h_1-Z{)}^2}$$

(14)

$T_0$ is relatively small compared to $1/f_t$ in practical milling and the influence of vibration on $T_0$ is limited. Therefore, it can be assumed that:

$$\begin{gathered} T_0\approx 2\left(t_3-t_2\right) \\ {h}_1=\frac{T_0{f}_{t}Z}{2} \end{gathered}$$

(15)

Substituting Equation (15) into Equation (14) and rearranging:

$$F_t\ge \frac{1}{2}{V}_F{T}_0+\sqrt{\left(2-T_0{f}_t\right)ZR-\left(1-\frac{T_0{f}_t}{2}\right)Z^2}$$

(16)

When the tool and cutting parameters suit for the vibration signal sampling method based on the TPF, there is always $R\gg Z$. Thus, Equation (16) can be further scaled to derive a more stringent range for $F_t$:

$$F_t\ge \frac{1}{2}{V}_F{T}_0+\sqrt{\left(2-T_0{f}_t\right)ZR}$$

(17)

It can be observed from the derivation process that the lower limit provided by Equation (17) is higher than the theoretical value. However, such a deviation is allowed for guiding the setting of $F_t$. Moreover, compared to the lower limit calculation method given by Equation (12), Equation (17) exhibits significantly higher feasibility in engineering application practices.

The cutting parameter tuning process for the vibration signal sampling method based on the TPF is illustrated in Figure 7. The tool’s structural parameters, such as the number of teeth, are considered predetermined. First, determine the frequency range of the vibration signal that is intended to be preserved in the surface topography. Next, based on the number of teeth of the tool and the determined frequency range, select the spindle speed $\omega$ according to Equation (11). Then, assign an initial value to the feed rate $V_F$ with set $\omega$ and solve for the effective cutting time per tooth $T_0$ by Equation (5). The fourth step is to estimate the maximum peak-to-valley amplitude $Z$ based on former machining experience and substitute it along with the calculated $T_0$. into Equation (17) to determine whether the current $F_t$ meets the requirements. Finally, iterate to determine the appropriate feed rate $V_F$.

After clarifying the cutting parameter requirements for vibration signal sampling based on the TPF, the feasibility of this method was validated through simulation. The tool radius and number of teeth are set to $R=10 \left(\mathrm{m}\mathrm{m}\right)$ and $N_t=2$, respectively, and these values remain constant across the following three sets of simulations. A sinusoidal vibration signal with a preset frequency of 20 Hz and an amplitude of 0.01 mm is defined as $V_x\left(t\right)=0.01sin40\pi t$. According to the process shown in Figure 7, the spindle speed is set to $\omega =300\pi (\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s})$, the maximum peak-to-valley amplitude of the vibration signal is $Z=0.02 (\mathrm{m}\mathrm{m})$, and the feed rate is set to $V_F=200\left(\mathrm{m}\mathrm{m}/\mathrm{s}\right)$. The surface topography simulation results of a part of the workpiece at a certain height after flank milling and the analysis in the frequency domain are shown in Figure 8a and Figure 8b, respectively. After performing the Fourier transform on the surface topography, the amplitude at the frequency of 20 Hz is found to be 0.0098 mm. The simulation results indicate that the exact frequency and amplitude of the vibration signal can be extracted from the surface topography. Therefore, the vibration signal sampling based on the TPF is verified feasible.

It is worth mentioning that the impact of the TPF on the surface topography is also reflected in the frequency-domain analysis results. However, this does not affect the accuracy of the extracted vibration signal frequency-domain characteristics. Since the TPF acts as the actual sampling frequency for the vibration signal, a low-pass filter with a bandwidth of $f_t/2$ will subsequently be applied to the surface topography.

Further simulations were conducted to verify the impact of TPF on the vibration signal retained in the surface topography. Keeping the remaining parameters constant, the vibration signal frequency is altered to obtain the relationship between TPF and the acquired vibration signal frequency. Simulations were conducted under the conditions of vibration signal frequencies of 80 Hz, 120 Hz, and 160 Hz, respectively. The results are shown in Figure 9. When the vibration frequency is 80 Hz, TPF is less than 5 times the vibration frequency, resulting in a decrease in the amplitude of the vibration signal to 0.0083 mm. Moreover, the signal waveform becomes distorted due to the reduced number of sampling points within the vibration cycle, leading to noticeable components in the analysis results at other frequencies. When the vibration frequency is 120 Hz, the accuracy of the vibration signal obtained from the surface topography further decreases. The amplitude of the obtained signal at this point is only 0.0066 mm, with further increased components at other frequencies. When the vibration frequency is 160 Hz, TPF is less than twice the vibration frequency, which does not meet the requirements of the sampling theorem. The frequency-domain analysis results of the surface topography show frequency aliasing occurring at this point. The fundamental peak reflecting the overall frequency-domain characteristics of the signal appears at 140 Hz. The vibration signal acquired from the surface topography is severely distorted, with significant deviations in both frequency and amplitude.

The final part of the simulation experiment was conducted to verify the influence of feed rate per tooth $F_t$ on the vibration signal retained in the surface topography. Assuming the vibration is with a frequency of 20 Hz and an amplitude of 0.01 mm, expressed as $V_x\left(t\right)=0.01sin40\pi t$. The spindle speed $\omega$ remains constant. Since the vibration signal is assumed to be known, the accurate theoretical range of $F_t$ can be calculated according to Equation (12), yielding $F_t\ge 0.1312 \mathrm{m}\mathrm{m}$. Because TPF is already determined, the range of $F_t$ can be converted into a range of feed rates as $V_F\ge 2361 \mathrm{m}\mathrm{m}/\mathrm{m}\mathrm{i}\mathrm{n}$. Six Simulations were conducted at $V_F=1500, 2000, 2500, 3000, 4000, 5000,$ and the results of the six simulations are shown in Figure 10. The blue curve fitted based on the simulation results represents the amplitude accuracy changing with feed speed. The results show that the accuracy of the vibration amplitude extracted from the surface topography decreases sharply when the feed speed does not meet the required range and only when $V_F\ge 2361 \mathrm{m}\mathrm{m}/\mathrm{m}\mathrm{i}\mathrm{n},$ the accuracy of the vibration amplitude reaches 90%.

The feasibility and sampling accuracy of the proposed vibration signal sampling method based on the TPF was verified through three sets of simulations conducted under appropriate cutting parameters. The principles and methods discussed in this section for extracting vibration signals from surface topography are indispensable for achieving the OMA of machine tools based on flank-milled surfaces. The investigations ensure that the vibration response signals remain undistorted which serve as the primary subjects of analysis in subsequent analysis. It is a crucial prerequisite for ensuring the accuracy of modal parameter identification that the vibration signals remain undistorted.

4. Cepstrum-Based OMA

4.1. Characteristics of Vibration Signals in Quefrency Domain

Vibration signals $y\left(t\right)$ in the time domain can be expressed as the convolution of excitation $e\left(t\right)$ and system transfer function $h\left(t\right)$ under the assumption of linearity and time-invariance of the machine tool system:

$$y(t)=e(t)\ast h(t)$$

(18)

Cepstrum is defined as the Fourier inverse transform of the natural logarithm of the spectrum. The definition of complex cepstrum $C_c\left(t\right)$ is as follows:

$$C_c(t)={\Gamma }^{-1}\left\{\mathit{ln}\left[\Gamma (f(t))\right]\right\}$$

(19)

where $\Gamma \left(\cdot \right)$ and ${\Gamma }^{-1}\left(\cdot \right)$ denote the Fourier transform and inverse transform, respectively. $f\left(t\right)$ represents the signal in the time domain. In this definition, the independent variable of the complex cepstrum is named quefrency and has the same dimension as time $t$. To distinguish between the quefrency domain and the time domain, variable $\tau$ will be used to represent quefrency. Substituting Equation (18) into Equation (19), the complex cepstrum of the vibration signal is given by:

$$\begin{array}{ll}\hfill C_c(\tau )& ={\Gamma }^{-1}\left\{\ln [E(f)\cdot H(f)]\right\}\\ & ={\Gamma }^{-1}\left\{\ln [E(f)]\right\}+{\Gamma }^{-1}\left\{\ln [H(f)]\right\}\\ & =C_{c,e}(\tau )+C_{c,h}(\tau )\end{array}$$

(20)

where $E\left(f\right)$ and $H\left(f\right)$ are the Fourier transforms of the excitation $e\left(t\right)$ and the system transfer function $h\left(t\right)$, respectively. $C_{c,e}\left(\tau \right)$ and $C_{c,h}\left(\tau \right)$ are the complex cepstrum of the excitation and the system transfer function, respectively. Equation (20) demonstrates the most significant feature of the cepstrum: it transforms the convolution in the time domain or the multiplication in the frequency domain into a linear addition in the quefrency domain. This transformation facilitates the removal of the excitation component from the vibration signal, therefore retaining only the transfer function component that contains the modal information.

Rewriting the Fourier transform pairs in exponential form yields:

$$\begin{array}{ll}\hfill C_c(\tau )& ={\Gamma }^{-1}\left\{\ln [E(f)\cdot H(f)]\right\}\\ & ={\Gamma }^{-1}\left\{\ln [A_e(f)\cdot e^{j{\varphi }_e(f)}]\right\}+{\Gamma }^{-1}\left\{\ln [A_h(f)\cdot e^{j{\varphi }_h(f)}]\right\}\\ & ={\Gamma }^{-1}\left\{\ln [A_e(f)\cdot A_h(f)]\right\}+{\Gamma }^{-1}[j{\varphi }_e(f)+j{\varphi }_h(f)]\end{array}$$

(21)

where $A_e\left(f\right)$ and $A_h\left(f\right)$ represent the amplitude spectrum of the excitation and the system transfer function in the frequency domain, respectively, while ${\varphi }_e\left(f\right)$ and ${\varphi }_h\left(f\right)$ represent their corresponding phase spectrum. Equation (21) analyses the composition of the complex cepstrum from another perspective, starting from the amplitude and phase spectrum in the frequency domain. In cepstrum analysis, the real cepstrum $C_r\left(\tau \right)={\Gamma }^{-1}\left\{\mathit{ln}\left[A\left(f\right)\right]\right\}$ is defined to represent the characteristics of the amplitude spectrum in the quefrency domain. According to Equation (21), directly calculating the complex cepstrum requires the signal to satisfy phase unwrapping. However, the commonly encountered stationary signals in practical applications do not satisfy this requirement. Therefore, the complex cepstrum is usually obtained indirectly by first calculating the real cepstrum [26]. The relationship between the complex cepstrum and the real cepstrum can be expressed as:

$$C_c\left(\tau \right)=\left\{\begin{array}{cc}\begin{array}{c}2C_r\left(\tau \right)\\ C_r\left(\tau \right)\\ 0\end{array}& \begin{array}{c}\tau >0\\ \tau =0\\ \tau <0\end{array}\end{array}\right.$$

(22)

Therefore, the analysis of signal characteristics in the quefrency domain and subsequent cepstrum editing mainly focus on the real cepstrum. The real cepstrum $C_{r,y}\left(\tau \right)$ corresponding to the vibration $y\left(t\right)$ of the machine tool system can be expressed as:

$$\begin{array}{ll}\hfill C_{r,y}(\tau )& ={\Gamma }^{-1}\left\{\ln [A_y(f)]\right\}\\ & ={\Gamma }^{-1}\left\{\ln [A_e(f)\cdot A_h(f)]\right\}\\ & =C_{r,e}(\tau )+C_{r,h}(\tau )\end{array}$$

(23)

where $C_{r,e}\left(\tau \right)$ and $C_{r,h}\left(\tau \right)$ represent the real cepstrum of the excitation and the system transfer function, respectively. As a real cepstrum formed by linear superposition, the distribution characteristics of $C_{r,y}\left(\tau \right)$ can be obtained by separately analyzing $C_{r,e}\left(\tau \right)$ and $C_{r,h}\left(\tau \right)$.

For the transfer function of a linear time-invariant multi-degree-of-freedom system, its real cepstrum $C_{r,h}\left(\tau \right)$ is composed of the superposition of the real cepstrum $C_{r,sDoF}\left(\tau \right)$ of multiple single-degree-of-freedom systems, which can be analytically expressed. Thus, $C_{r,h}\left(\tau \right)$ can be written as:

$$\begin{gathered} C_{r,h}(\tau )=\sum C_{r,sDoF}(\tau ) \\ {C}_{r,sDoF}(\tau )=\frac{p^n}{n}+\frac{{\widehat{p}}^n}{n}=\frac{2e^{-\sigma n\Delta \tau }}{n}\mathit{cos}[n\mathit{arg}(p)] \end{gathered}$$

(24)

Here $p$ and $\widehat{p}$ represent a pair of conjugate poles, $n$ is the number of sampling points, $\sigma$ is the damping coefficient,$\mathit{arg}(p)$ represents the phase angle of $p$, and $\Delta \tau$ is the sampling time interval. The distribution of $C_{r,h}(\tau )$ in the quefrency domain exhibits two main characteristics. First, it demonstrates a low-pass filtering property, where non-zero values are predominantly located at lower quefrencies. Second, its waveform depicts damped oscillatory decay with its absolute magnitude eventually converging to zero. These distribution features provide theoretical support for the removal of corresponding excitation components through cepstrum editing.

In flank milling the cutting force acting as excitation comprises two components: periodic and stochastic. The periodic component, as the main body of excitation, is generated by the regular shearing and compressing between the cutting edge and the workpiece material based on the TPF. It is characterized by higher energy and exhibits a sawtooth waveform. The stochastic component, acting as noise, is caused by uncertain factors such as tool wear and randomly distributed hard points in the workpiece material. It is characterized by lower energy and can be considered to be a zero-mean stationary stochastic process. As clarified in Section 3, flank milling must be carried out under conditions of a relatively large feed rate per tooth to realize vibration sampling. Hence, it can be assumed that the energy of the periodic component is much greater than that of the stochastic component at this point. Under this premise, the excitation $e(t)$ can be rewritten as:

$$e(t)=\sum _{i=-\infty }^{+\infty }\frac{A}{ji\pi }\mathit{cos}(i\pi )e^{2ji\pi f_{t}t}$$

(25)

where $j^2=-1$, $A$ represents the amplitude of the sawtooth wave, and $f_t$ is the TPF. The amplitude spectrum of the excitation $A_e(f)$ can be obtained from the Fourier transform of the periodic signal. Considering that the negative half of the spectrum does not have physical significance, only the positive half spectrum $A_e(f{)}^{+}$ is studied:

$$A_e(f{)}^{+}=\sum _{i=1}^{+\infty }\frac{2A}{i}\cdot \delta (f-i{f}_t)$$

(26)

$\delta (f-i{f}_t)$ represents the impulse at $f=i{f}_t$. Equation (26) shows that the sawtooth wave excitation, as a typical periodic signal, also has discrete characteristics in its amplitude spectrum. Moreover, there are significant differences in the amplitudes at the fundamental frequency and each harmonic frequency, and the spectral lines appear steep. Taking the natural logarithm of the amplitude spectrum to obtain $\mathit{ln}[A_e(f{)}^{+}]$:

$$\mathit{ln}[A_e(f{)}^{+}]=\left\{\begin{array}{cc}\begin{array}{c}\mathit{ln}A-\mathit{ln}i+\mathit{ln}2\\ a,a\to -\infty \end{array}& \begin{array}{c}f=i{f}_t,i\in N^{+}\\ f\ne i{f}_t,i\in N^{+}\end{array}\end{array}\right.$$

(27)

The logarithmic operation reduces the differences in amplitude at the fundamental frequency and each harmonic frequency, making the spectral envelope smoother. Additionally, Equation (27) indicates that under conditions of large cutting forces where the amplitude A is large, it can be assumed for a finite series expansion of the sawtooth wave that the amplitudes of the discrete spectrum in the natural logarithm of the amplitude spectrum are approximately equal to $\mathit{ln}A$. The real cepstrum of the excitation $C_{r,e}(\tau {)}^{+}$ can be calculated through the Fourier inverse transform of the natural logarithm of the amplitude spectrum:

$$\begin{array}{ll}\hfill C_{r,e}{(\tau )}^{+}& ={\Gamma }^{-1}\{\ln [A_e{(f)}^{+}]\}\\ & =0.5a\delta (\tau )+{\displaystyle \sum _{i=1}^m\frac{\ln A-\ln i+\ln 2-a}{2\pi }{e}^{j2\pi i{f}_t\tau }}\end{array}$$

(28)

where $m$ represents the number of terms in the finite series expansion of the Fourier series. Equation (28) demonstrates two distribution characteristics of the real cepstrum of the sawtooth wave excitation. First, it has significant components at zero quefrency, and second, the spectral lines at non-zero quefrencies are discretely periodic, with non-zero amplitudes only near $\tau =k/f_t,k\in N^{+}$. The distribution characteristics of the real cepstrum stem from the similarity between the inverse Fourier transform and the Fourier transform. They can both be regarded as spectral analyses, leading to the interpretation of the real cepstrum as a spectrum in the quefrency-domain of the frequency spectrum. The natural logarithm of the magnitude spectrum $\mathit{ln}[A_e(f{)}^{+}]$ exhibits a constant component due to the logarithmic calculation, which is also the reason for the first distribution characteristic. Additionally, since the spectral lines in the frequency spectrum are periodically discrete, the spectral lines in the real cepstrum are also periodically discrete. From the above analytical perspective, it can be considered that the spectral line distribution of periodic signals in the quefrency domain is the same as that of $C_{r,e}(\tau {)}^{+}$.

To further illustrate the differences in distribution characteristics between $C_{r,e}\left(\tau \right)$ and $C_{r,h}\left(\tau \right)$, a vibration response of a third-order modal system under sawtooth wave excitation was obtained through simulation. The real cepstrum $C_{r,y}\left(\tau \right)$ of the vibration, as well as the excitation of real cepstrum $C_{r,e}\left(\tau \right)$ and the real cepstrum of the system transfer function $C_{r,h}\left(\tau \right)$ were calculated. The results are shown in Figure 11. Results visually demonstrate significant differences in the distribution of the real cepstrum between the excitation and the system transfer function. Moreover, due to the deconvolution property of the cepstrum, the vibration response curve is synthesized by the superposition of the first two cepstrum. Hence, separating the components corresponding to the periodic excitation from the real cepstrum of the time-domain vibration response does not require complex mathematical processing.

4.2. Modal Parameters Identification Based on Cepstrum Editing

4.2.1. Identification Process

The widely used modal parameter identification process based on cepstrum editing includes four main steps: cepstrum calculating, cepstrum editing, spectrum reconstruction, and rational fractional polynomial fitting [35]. Combining the characteristics of the vibration signal sampling method based on the TPF and sawtooth wave excitation in the flank milling proposed in this paper, the modal parameters identification process based on surface topography is depicted in Figure 12. In this process, the step of editing the cepstrum by liftering or other means based on the differences in the distribution characteristics of each component in the cepstrum determines the integrity and accuracy of the modal information. This step is the core of the identification process. Therefore, the study of the modal parameter identification process based on cepstrum editing should focus on cepstrum editing methods. Other relatively mature steps in the identification process, such as spectrum reconstruction and rational fractional polynomial fitting can be found in the reference and are not reiterated here.

4.2.2. Cepstrum Editing Method for Extracting Transfer Function

Common techniques used for cepstrum editing to extract transfer function components include notch liftering and window functions. The characteristic of these editing techniques is to preserve the spectral lines in the low quefrency range that may contain modal information while editing only the spectral lines in the high quefrency range. Equation (28) indicates that these editing techniques are appropriate when the frequency of the periodic excitation is low. However, as the frequency of the periodic excitation increases, the spectral lines corresponding to the periodic excitation component in the cepstrum gradually shift to the left. Eventually, they may overlap with the spectral lines in the low quefrency range that contain modal information, therefore affecting the identification accuracy. In the milling process, TPF is often relatively high, while the modal damping of the machine tool system is low. Therefore, the situation described above is highly likely to occur. In addition, Equation (28) also indicates the presence of interference from periodic excitation at the quefrency of 0, which should also be eliminated during cepstrum editing.

First, the influence of random disturbances on identification accuracy should be eliminated. Flank milling is conducted under certain axial cutting depths. Therefore, in Section 2, a discretized modeling approach was adopted to divide the milled surface topography into $N$ layers. The $N$ sets of vibration data retained in topography obtained through a single flank milling exhibit temporal synchronicity. This implies that leveraging this advantage of flank milling enables the acquisition of multiple sets of vibration data at a lower experimental cost, therefore providing favorable conditions for employing the time-domain synchronous averaging method to eliminate random noise interference. Denoting the cepstrum of vibration obtained by transforming the two-dimensional topography data at different heights as $C_{r,y_i}(\tau )$ and then there is:

$$C_{r,\stackrel{-}{y}}(\tau )=\frac{1}{N}\sum _{i=1}^N{C}_{r,y_i}(\tau )$$

(29)

where $C_{r,\stackrel{-}{y}}(\tau )$ represents the cepstrum of the vibration after removing random disturbances.

After removing random disturbances by time-domain synchronous averaging, the next step is cepstrum editing to eliminate the periodic excitation component from the cepstrum. Considering that the periodic excitation component may overlap with the transfer function component, a modified window function is proposed to remove the periodic excitation component. Define an isolation window $\tau \in (0,s{\tau }_0],s\in N^{+}$ where the length coefficient $s$ is determined based on prior experience design principles and manually set. However, when setting the parameters, the boundaries of the window function need to satisfy:

$$C_{r,h}(s{\tau }_0)=0,{\tau }_0=\frac{1}{f_t}$$

(30)

Equation (28) indicates that the amplitude cepstrum of the periodic excitation at $\tau =k{\tau }_0,k\in N^{+}$ should theoretically be equal. Therefore, the amplitude $A_{c,e}$ of the periodic excitation at $\tau =k{\tau }_0,k\in N^{+}$ can be calculated based on the portion of $C_{r,\stackrel{-}{y}}(\tau )$ outside the isolation window:

$$A_{c,e}=\frac{1}{m-s}\sum _{k=s+1}^m{C}_{r,\stackrel{-}{y}}(k{\tau }_0)$$

(31)

where $m$ represents the integer multiples of ${\tau }_0$ within the range covered by $C_{r,\stackrel{-}{y}}(\tau )$. It should be noted that, theoretically, $m$ tends to infinity, but in practical applications, it is limited by the sampling frequency. Therefore, a constant is used here as an approximation.

Once the isolation window is set, the cepstrum $C_{r,\stackrel{-}{y}}(\tau )$ can be divided into three regions for editing to remove the periodic excitation components. First, eliminate the influence of the periodic excitation at zero quefrency as indicated by Equation (28). It is worth noting that in practical data acquisition and processing, Discrete Fourier Transform (DFT) or Fast Fourier Transform (FFT) is commonly used. It results in discrepancies between the actual value to be subtracted during the cepstrum editing process and the theoretical value. However, this falls within the scope of data processing and will not be further elaborated here. For the part of the cepstrum within the isolation window, compensation pulses with an amplitude of $-A_{c,e}$ and a width of $1/A_{c,e}$, are added at $\tau =k{\tau }_0, k=\mathrm{1,2},3,\cdots ,s$ to remove the influence of periodic excitation within the window. Considering possible interference from other vibration sources in the actual milling process on the cepstrum, the part outside the isolation window is zeroed out. The edited cepstrum $C_{r,y-edited}(\tau )$ is thus obtained:

$$C_{r,y-edited}(\tau )=\left\{\begin{array}{cc}\begin{array}{c}{C}_{r,\stackrel{-}{y}}(\tau )-0.5a\\ C_{r,\stackrel{-}{y}}(\tau )-A_{c,e}\cdot \delta (\tau -k{\tau }_0)\\ 0\end{array}& \begin{array}{c}\tau =0\\ 0<\tau \le s{\tau }_0,k=\mathrm{1,2},3,\cdots ,s\\ \tau >s{\tau }_0\end{array}\end{array}\right.$$

(32)

After extracting the components corresponding to the system transfer function through cepstral editing, subsequent modal identification is completed through spectrum reconstruction and rational fraction polynomial fitting.

The cepstrum editing procedure is shown in Figure 13. The proposed modified window function for cepstrum editing is based on flank milling. This method leverages the characteristics of flank milling where the periodic excitation is relatively fixed in form, predictable in period, and has a large excitation amplitude. Additionally, multiple sets of data can be obtained from a single milling operation, reducing experimental costs and eliminating random interference. By utilizing a clear model of the excitation in the quefrency domain, the method accurately removes periodic excitation components. Compared to traditional window functions, the proposed window function has minimal impact on the low quefrency part of the cepstrum containing modal information. Therefore, in subsequent modal parameter calculations, there is no need to subtract additional damping introduced by the window function.

5. Experiment Verification

5.1. Verification of Proposed Cepstrum-Based OMA

Experiments were conducted on a certain AC swivel-type five-axis milling machining center to validate the proposed cepstrum-based OMA based on flank-milled surface topography. Considering the potential impact of tool wear, an aluminum alloy with a hardness significantly lower than that of the tool was selected as the workpiece for the flank milling experiments. Thus, it can be assumed that the change in tool wear during the experiments is minimal and negligible. Consequently, the impact of tool wear on the vibration signals extracted from the surface topography can be considered to be an addition of a constant offset. Although tool wear is not uniform along the tool, it does not vary with time during the experiments. Therefore, its impact can be removed using high-pass filtering.

After installing a five-teeth end mill with a diameter of 16 mm on the spindle of the machine tool, hammer impact tests (HIT) were conducted at the tooltip point to measure the experimental modal parameters of the machine tool, which serves as the control group for subsequent experiments. The HIT instrumentations are shown in Figure 14a. The hammer impact point and the accelerometer were located on opposite sides of the tooltip. Machine tool modal parameters of the first three orders were measured through the experiment, as shown in Table 1.

The proposed modal parameters identification method was first verified. The process of flank milling the workpiece is shown in Figure 14b. The spindle speed $\omega$ and the feed rate $V_F$ were set to be 200π (rad/s) and 150 (mm/s), respectively. The plane with a length of 150 mm was milled under the conditions of an axial depth of 2 mm and a radial depth of 1 mm. The surface topography of the milled surface was measured using a white light interferometer, as shown in Figure 15. Due to the workspace size of the white light interferometer, the milled surface of 150 mm in length was divided into multiple sections for individual measurements and then assembled. As a result, the final measured length is slightly less than 150 mm.

According to the cepstrum editing method proposed in Section 4, two-dimensional topography profile data at 10 different positions along the tool axis are initially selected to eliminate random disturbances. Considering the relative magnitudes between TPF and the natural frequencies, the isolation window length coefficient $s$ was set to 50. The magnitude of $A_{c,e}$ obtained by time-domain synchronous averaging was 0.106. A comparison before and after cepstrum editing is shown in Figure 16. In Figure 16a, the amplitudes of the periodic excitation component are not identical for two main reasons. First, to prevent tool wear or breakage, the cutting parameters are set in such a way that the magnitude of the cutting force does not satisfy the approximation condition discussed based on Equation (27), leading to deviation in the amplitudes of the periodic excitation component in low quefrency range. Second, there are other sources of non-random disturbances.

After cepstrum editing, the frequency spectrum was reconstructed. The initial amplitude spectrum data obtained from HIT, the frequency response function (FRF) fitted from initial data, and the FRF reconstructed by the proposed modal identification method are shown in Figure 17. The modal parameters identification results obtained by rational fractional polynomial fitting are presented in Table 2. The sampling time was calculated to be 1 (s) based on the milled surface length and the feed rate. Therefore, according to the sampling theorem, the precision achievable for the identification of natural frequencies was 1 (Hz). Compared to the results of HIT, the proposed modal parameters identification method yields slightly increased natural frequencies for all orders. This is attributed to the variation in the dynamic stiffness of the machine tool during milling. Therefore, the proposed method can be considered accurate and effective for the identification of natural frequencies. As for the identification results of modal damping ratios for each order, significant differences exist between the proposed method and HIT. This is because the machine tool system experiences changes in contact states between components and preload forces due to the cutting force.

Conclusions can be drawn from the validation experiments that the proposed cepstrum editing method effectively reduces random disturbances and significantly diminishes periodic excitation components. The modal parameter identification method based on the surface topography obtained was verified.

5.2. Influence of Cutting Parameters on Vibration Signal Sampling

After verifying the effectiveness of the proposed method, the reconstructed FRF obtained by the proposed method was taken as the vibration signal to be extracted from the surface topography to verify the influence of cutting parameters on vibration signal sampling. After obtaining the highest frequency of vibration components through HIT, it can be obtained from Equations (11) and (17) that when $Z=0.005$ (mm), the necessary conditions for complete acquisition of the vibration signal are $\omega \ge 189\pi$ (rad/s) and $F_t\ge 0.285$ (mm).

First, validate the influence of TPF on the vibration signal sampling results. By adjusting the feed rate $V_F$ to maintain a constant feed rate per tooth $F_t=0.3$ (mm). Vibration signal sampling and reconstruction of FRF were carried out at spindle speeds $\omega$ of 220π, 160π, 120π, 80π, and 40π. The FRF curves reconstructed under different spindle speeds are shown in Figure 18. When $\omega =220\pi$, TPF is greater than five times the highest vibration frequency. Hence, it can be considered that the amplitude spectrum of the acquired signal is accurate at this spindle speed. However, when $\omega =160\pi , 120\pi , 80\pi$, TPF decreases successively while it still satisfies the sampling theorem that the sampling frequency is greater than twice the vibration frequency. Therefore, although the accuracy of the amplitude decreases gradually, the natural frequencies indicated by the FRF remain accurate. When $\omega =40\pi$, TPF is less than twice the maximum vibration frequency, which does not satisfy the sampling theorem. The results show that each mode was superimposed, leading to severe distortion of the obtained FRF. Significant deviations exist in both the indicated natural frequencies and amplitudes.

After verifying the influence of TPF on the proposed vibration signal sampling method, the influence of feed rate per tooth on the proposed method was further discussed through experiments. Maintaining a constant spindle speed $\omega$ of 220π (rad/s), experiments were conducted with feed rate per tooth $F_t$ of 0.333, 0.267, 0.2, 0.133, and 0.067, respectively. FRF curves reconstructed under different $F_t$ are shown in Figure 19. The influence of $F_t$ on the vibration signal sampling method mainly manifests in the amplitude. As $F_t$ decreases, the surface roughness decreases gradually, resulting in a decrease in the amplitude of FRF across the entire frequency range, especially noticeable near the natural frequencies of each mode. Experimental results indicate that when $F_t$ is less than 0.285, there is a significant decrease in amplitude accuracy, leading to deviations in the identified parameters.

The two sets of validation experiments, respectively, explored the influence of TPF and feed rate per tooth on the vibration signal sampling based on TPF, confirming the necessary conditions stated in Section 3. It is worth noting that the variations in cutting parameters result in changes in the working state of the machine tool, leading to changes in operational modal parameters. However, the influence is assumed to be negligibly small compared to the influence on vibration signals sampling, which results in an obvious deviation of amplitude during experiments. This assumption may have some impact on the assessment of identification accuracy, but it does not affect the exploration of the influence of cutting parameters on vibration signals sampling based on TPF.

6. Conclusions

Aiming to address discrepancies between the excitation method and the actual operating conditions and the unreasonable collection of vibration signals in current OMA for CNC machine tools, a method is proposed in this paper for identifying the operational modal parameters of machine tools based on the milled surface topography and cepstrum analysis. This method not only ensures the alignment between the excitation method and the actual operating conditions by extracting vibration signals from the surface topography but also captures the intended displacements at the tooltip. The proposed method offers low experimental costs and does not require operators to possess additional knowledge of experimental instrument operation, making it easier to promote and apply widely. Moreover, the identification results closely resemble the actual state of the machine tool system. Combining surface topography with identification results also provides a more intuitive reference for optimizing surface quality through the adjustment of cutting parameters in subsequent processes.

The highlights of this study in terms of theory and methodology can be summarized as follows:

(a)

In comparison to numerous existing pairwise relationship studies, this study systematically considers the interaction relationship between modal parameters, vibration, and surface topography. It makes the proposed modal analysis method a valuable tool for future systematic studies on the interaction mechanisms among these factors.

(b)

A theoretical framework for vibration signal sampling based on tooth passing frequency and flank-milled surface was proposed. The framework also provides theoretical guidance for optimizing cutting parameters in flank milling aimed at workpiece surface quality enhancement.

(c)

Under the premise of known periodic excitation and estimable excitation amplitude, a cepstrum editing method based on a modified window function is proposed. The cepstrum analysis method adopted in this paper provides a new perspective for modeling and deconstructing surface topography in machining.

Furthermore, limitations were also identified during the research process of the proposed method. As the natural frequency of the mode to be identified increases and the amplitude under unit cutting force grows, it is necessary to use more stringent cutting parameters to directly obtain accurate and complete vibration signals from surface topography. This poses a significant challenge to machine tool performance. Further research may focus on vibration signal reconstruction algorithms to explore how to expand the identifiable frequency range under more relaxed cutting parameter requirements.