Output-Only Time-Varying Modal Parameter Identification Method Based on the TARMAX Model for the Milling of a Thin-Walled Workpiece

Abstract

The process parameters chosen for high-performance machining in the milling of a thin-walled workpiece are determined by a stability prediction model, which needs accurate modal parameters of the machining system. However, the in-process modal parameters are different from the offline modal parameters and are difficult to precisely obtain due to material removal. To address this problem, an accurate time-dependent autoregressive moving average with an exogenous input (TARMAX) method is proposed for the identification of the modal parameters in the milling of a thin-walled workpiece. In this process, a TARMAX model considering external force excitation is constructed to characterize the actual condition in the milling of a thin-walled workpiece. Then, recursive method and sliding window recursive method are used to identify TARMAX model parameters under time-varying cutting conditions. Subsequently, a three-degree of freedom (3-DOF) time-varying structure numerical model under theoretical milling forces and white-noise excitation is established, and the computational results show that the predicted natural frequencies using the proposed method are in close agreement with the simulated values. Finally, several experiments are designed and carried out to validate the effectiveness of the proposed method. The experimental results show that the predicted accuracy of the proposed method using actual cutting forces is 95.68%. Good agreement has been drawn in the numerical simulation and machining experiments. Our further research objectives will focus on the prediction of the damping ratios, modal stiffness, and modal mass.

1. Introduction

Thin-walled workpieces such as aero-engine blades, casings, and impellers are widely used in the aerospace industry [1] because of their light weight, high strength, and heavy load-bearing characteristics [2,3]. However, machining chatter and dimensional error caused by the problem of low rigidity occur in the milling of a thin-walled workpiece, which can affect the surface quality, machining productivity, and tool life [4,5,6]. In the milling of thin-walled workpieces, the methods of milling stability prediction are extensively investigated for avoiding chatter [7,8,9]. Generally, accurate stability lobe diagrams calculated by solving a time period delay differential equation are critical for milling stability [10,11,12]. In this process, the modal parameters of the pressing system are important factors affecting accurate stability lobe diagrams. In particular, the modal parameters are time-varying due to material removal and the tool–workpiece contact position in the milling of a thin-walled workpiece [12,13,14,15], which will lead to low cutting process parameters in milling and will reduce the machining productivity [10]. Therefore, it is necessary to accurately identify the time-varying modal parameters of thin-walled workpieces in milling in order to improve their efficiency and quality.

For extracting the modal parameters, the common method that is used is the hammer test. However, this method cannot be adopted to obtain modal parameters at in-process machining [16,17]. To deal with this problem, operational modal analysis (OMA) is used to identify the dynamic modal parameters in the milling of a thin-walled workpiece. For the efforts of identifying the modal parameters in milling, OMA can be divided into two categories: frequency domain modal identification and time domain modal identification. For frequency domain modal identification, frequency response functions or power density functions are used to obtain time-varying modal parameters, which are mainly applied in large-scale bridge and building structure health monitor fields. Nevertheless, in milling, excitation signals mainly contain strong harmonic excitation and white-noise excitation due to spindle rotation and material removal [18], which may result in the harmonic components being mistaken for structural modes in the OMA frequency domain identification method. Liang et al. [19] eliminated the response generated by harmonic excitation by the properties of white-noise excitation signals. Kiss et al. [20] employed a comb filter to eliminate the response generated by harmonic excitation. Yuan et al. [21] used Kalman filters, which effectively attenuate spindle frequency and its harmonics. However, if the harmonic fundamental frequency cannot be determined before using these methods, even a small error will be significantly amplified in the high-order harmonic components and cause filter failure. For this, Liu et al. [18] proposed a modified least square method to fit harmonics with multiple fundamental frequencies and to extract the modal parameters of the workpiece–tool system. Weijtjens et al. [22] estimated the structural dynamics of the system by using transfer functions that could be unaffected by periodic harmonic excitations. Wan et al. [23] expressed the power spectral density matrix as a spectral decomposition form with modal parameters and an extract damping ratio using inverse Fourier transform, which ignored the natural frequency identification. These methods eliminated harmonic components in machining signals and identified the modal parameters, but they relied on time domain signals to decrease the conversion errors.

For the time domain modal identification, the modal parameters could be directly identified in the time domain and avoid signal errors compared with the frequency domain method. Li et al. [24] extracted machine tool modal parameters through the stochastic subspace identification (SSI) method. Then, the effectiveness of the SSI method was validated by Yan et al. [25]. However, the identification accuracy and the time-varying tracking ability of the SSI method need to be improved. Subsequently, Burney et al. [26] employed a time-dependent autoregressive moving average (TARMA) model to estimate the machine tool modal parameters. In addition, Kim et al. [27] applied TARMA to the milling operations. Then, Zaghbani et al. [28] adopted the TARMA approach to estimate the modal parameters in real-time during machining. Ma et al. [29] proposed a kernelized TARMA model that extended on the TARMA approach, and the computational efficiency was improved. However, these methods only identified modal parameters through white-noise signals. In the actual machining, Kang et al. [30] obtained the actual modal parameters by removing the harmonic components based on the statistical characteristics of the response, but the amount of calculations were too large, and the signal produced errors due to Fourier transform. Zhuo et al. [31] adopted the TARMA method to identify the modal parameters in the milling of a thin-walled workpiece, and false modes eliminated the convergence properties of the stability lobe diagram.

Through the above analysis, it was found that a number of methods for modal parameter identification are used in a wide field of processing. However, for time-varying processing conditions, especially the milling of thin-walled workpieces, the identification of time-varying modal parameters relevant to the prediction of the system stability is less well studied. To overcome this problem, a modified TARMAX modal parameters identification method considering the actual cutting force excitation was proposed in this paper. Taking the dynamic milling force and environment noise as excitation signals, the TARMAX modal parameter identification algorithm was derived based on “frozen time”, and the relationship between time-varying autoregressive coefficients and modal parameters was established. Subsequently, the white-noise excitation signal and the time-varying autoregressive coefficients were estimated by the least square method, and then, the time-varying modal parameters were determined. This paper is organized as follows. In Section 2, the TARMAX model based on milling force excitation was established, the recursive estimation of the TARMAX model parameters was identified, and recursive estimation of the time-varying coefficient sliding windows was determined. Then, in Section 3, a 3-DOF time-varying machining system numerical simulation model was established and the modal parameter identifications were obtained and compared under different noise and window numbers. Subsequently, the effectiveness and feasibility of the proposed method were validated by several machining experiments in Section 4. Finally, some conclusions are drawn in Section 5.

2. TARMAX Model Modal Identification Algorithm

TARMA methods are computationally small and can capture the time-varying properties of structures in real-time. In contrast with the identification methods in the frequency domain, by continuously introducing new data, the TARMA methods only need to modify the model parameters estimated in the previous step and can quickly model the system at the current moment [32]. This section is divided into the TARMAX model considering the actual milling force excitation, TARMAX model parameter recursive estimation, and sliding window recursive estimation.

2.1. TARMAX Model Considering Actual Milling Force Excitation

From the probability statistical chart, it can be found that the dynamic milling force signals were a superposition of the harmonic components and Gaussian white noise [18]. Therefore, the excitation exerted on the workpiece could be considered as a superposition of the theoretical milling force and Gaussian white noise. Subsequently, the time-varying dynamic equation of the machining system was expressed by an n-dimensional equation:

$$M(t)\ddot{x}(t)+C(t)\dot{x}(t)+K(t)x(t)=u(t)$$

(1)

where $M(t)\in {\mathrm{R}}^{n\times n}$, $C(t)\in {\mathrm{R}}^{n\times n}$, and $K(t)\in {\mathrm{R}}^{n\times n}$ are the mass, damping, and stiffness matrices of the machining system, respectively. t is continuous time. $\ddot{x}(t)\in {\mathrm{R}}^{n\ast 1}$, $\dot{x}(t)\in {\mathrm{R}}^{n\ast 1}$, and $x(t)\in {\mathrm{R}}^{n\ast 1}$ are acceleration, velocity, and displacement vectors, respectively. $u(t)\in {\mathrm{R}}^{n\ast 1}$ is the input excitation signals, which contain theoretical milling force $F(t)\in {\mathrm{R}}^{n\ast 1}$ and the environmental excitation $e(t)\in {\mathrm{R}}^{n\ast 1}$, namely, u(t) is defined as

$$u(t)=F(t)+e(t) e(t)\sim NID(0,{\mathbf{\sigma }}^2(t))$$

(2)

where e(t) is the environmental excitation, which can be treated as an uncorrelated white-noise excitation signal when the expectation and variance of the environmental excitation are 0 and σ²(t), respectively.

By continuously “freezing” the mass, damping, and stiffness matrices [33], the time-varying dynamic equation of Equation (1) can be transformed into TARMAX [32], which is given by

$$x[t]+{\displaystyle \sum _{i=1}^{N_a}{a}_i[t]x[t-i]}={\displaystyle \sum _{j=0}^{N_b}{b}_j[t]u[t-j]}={\displaystyle \sum _{j=0}^{N_b}{b}_{\mathrm{j}}[t](F[t-j]+e[t-j])}$$

(3)

where x[t] is the discrete displacement response; u[t] is the discrete excitation signal; F[t] is the discrete milling force; e[t] is the white noise sequence; a_i[t] and b_i[t] are the time-varying autoregressive coefficient and external excitation matrix, respectively; and N_a and N_b are the order of a_i[t] and b_i[t], respectively. Then, introducing backward shift operator z, and

$$z^{i}x[t]=x[t-i]z^i$$

(4)

Substituting Equation (4) into Equation (3), Equation (3) can be rewritten in the state-space form

$$A[z,t]x[t]=B[z,t]u[t]$$

(5)

where

$$\left\{\begin{array}{l}A[z,t]=I+{\displaystyle \sum _{i=1}^{N_a}{a}_i[t]z^i}\\ B[z,t]={\displaystyle \sum _{j=0}^{N_b}{b}_j[t]z^j}\end{array}\right.$$

(6)

Subsequently, the time-varying transfer function is obtained from Equation (5) as

$$H[z,t]=A{[z,t]}^{-1}\circ B[z,t]$$

(7)

Defining $z=e^{-\mathrm{j}\omega \Delta t}$, the “frozen time” transfer function of TARMAX is expressed as follows [34]

$$H[\omega ,t]=A{[e^{-\mathrm{j}\omega \Delta t},t]}^{-1}\circ B[e^{-\mathrm{j}\omega \Delta t},t]\approx \frac{B[e^{-\mathrm{j}\omega \Delta t},t]}{A[e^{-\mathrm{j}\omega \Delta t},t]}$$

(8)

where j is an imaginary number and △t denotes the sampling period. Considering Equation (8), in every moment t, assume that $\left|A[e^{-\mathrm{j}\omega \Delta t},t]\right|=0$, then, the “frozen time” transfer function poles p_r (r = 1, …, N_a) are determined, so the natural frequency of the continuous system Equation (1) can be calculated by [35]

$$f_r[t]=\frac{\left|\ln p_r[t]\right|}{2\mathsf{\pi }\cdot \Delta t}$$

(9)

For $\left|A[e^{-\mathrm{j}\omega \Delta t},t]\right|=0$, i.e.,

$$\left|A[e^{-\mathrm{j}\omega \Delta t},t]\right|=\left|I+{\displaystyle \sum _{i=1}^{N_a}{a}_i[t]z^i}\right|=0$$

(10)

To avoid a non-linear system of Equation (10), the roots of Equation (10) can be transformed into a generalized eigenvalue problem [35], so Equation (10) can be rewritten as

$$\left(D\right[t]-p_r[t]I)V_r[t]=0$$

(11)

where V_r[t] is an eigenvector, D[t] is the eigenvalues matrix, which is defined as

$$D[t]=\left[\begin{array}{cccc}0& I& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & I\\ -a_{N_a}& -a_{N_a-1}& \dots & -a_1\end{array}\right]$$

(12)

Therefore, if the time-varying autoregressive coefficients ɑ_i[t] can be accurately estimated, and the natural frequency $f_r[t]$ of the time-varying dynamic system can be solved by combining Equations (9) and (11).

2.2. TARMAX Model Parameters Recursive Estimation

Based on the milling force model [36], the theoretical milling force can be predicted. However, according to Equation (3), it can be seen that ɑ_i[t], b_i[t], and e[t] are unknown, and it is time-consuming to directly determine the ɑ_i[t] matrix based on the TARMAX model. Therefore, in this section, the least square method is used to improve the calculation efficiency.

2.2.1. Theoretical Milling Force Model

The natural frequencies in the x and y directions are lower than the natural frequencies in the z direction due to the clamping method and the machining mode. As the material is removed in milling, the workpiece is prone to chatter under the dynamic milling forces in the x and y directions. Then, in this paper, milling forces in the x and y directions are mainly considered and milling forces in the z direction are ignored. Therefore, in the milling of a thin-walled workpiece, the time-varying dynamic milling system can be expressed by an n-dimensional linear time periodic system, which is shown in Figure 1.

According to Figure 1, the dynamic milling force F[t] is represented as follows [36]

$$F\left[t\right]=a{K}_c\left(x\left[t\right]-x\left[t-T\right]\right)+a{f}_0$$

(13)

where a is the depth of cut, T denotes the time delay period, f₀ expresses the static force and K_c is the milling force coefficients. In addition, f₀ and K_c are defined as

$$\left\{\begin{array}{l}{K}_c={\displaystyle \sum _{j=1}^{n_1}g\left({\varphi }_j\left[t\right]\right)\left[\begin{array}{cc}{k}_{xx}& k_{xy}\\ k_{yx}& k_{yy}\end{array}\right]}\\ k_{xx}=-k_{tc}\sin \left({\varphi }_j\left[t\right]\right)\cos \left({\varphi }_j\left[t\right]\right)-k_{rc}\sin {\left({\varphi }_j\left[t\right]\right)}^2\\ k_{xy}=-k_{tc}\cos {\left({\varphi }_j\left[t\right]\right)}^2-k_{rc}\sin \left({\varphi }_j\left[t\right]\right)\cos \left({\varphi }_j\left[t\right]\right)\\ k_{yx}=k_{tc}\sin {\left({\varphi }_j\left[t\right]\right)}^2-k_{rc}\sin \left({\varphi }_j\left[t\right]\right)\cos \left({\varphi }_j\left[t\right]\right)\\ k_{yy}=k_{tc}\sin \left({\varphi }_j\left[t\right]\right)\cos \left({\varphi }_j\left[t\right]\right)-k_{rc}\cos {\left({\varphi }_j\left[t\right]\right)}^2\end{array}\right.$$

(14)

$$f_0={\displaystyle \sum _{j=1}^{n_1}g\left({\varphi }_j\left[t\right]\right)\left[\begin{array}{l}{f}_t\left(\begin{array}{l}-k_{tc}\sin \left({\varphi }_j\left[t\right]\right)\cos \left({\varphi }_j\left[t\right]\right)-k_{rc}\sin {\left({\varphi }_j\left[t\right]\right)}^2\\ k_{tc}\sin {\left({\varphi }_j\left[t\right]\right)}^2-k_{rc}\sin \left({\varphi }_j\left[t\right]\right)\cos \left({\varphi }_j\left[t\right]\right)\end{array}\right)\\ +\left(\begin{array}{l}-k_{te}\cos \left({\varphi }_j\left[t\right]\right)-k_{re}\sin ({\varphi }_j(t))\\ k_{te}\sin \left({\varphi }_j\left[t\right]\right)-k_{re}\cos \left({\varphi }_j\left[t\right]\right)\end{array}\right)\left({\varphi }_j\left[t\right]\right)\end{array}\right]}$$

(15)

where n₁ is number of teeth; f_t is feed per tooth in machining; k_tc, k_rc, k_te, and k_re are the milling force coefficients; g(ϕ_j[t]) is the cut-in function; and ϕ_j[t] is the angular position of the j-th tooth; then, ϕ_j[t] and g(ϕ_j[t]) are

$${\varphi }_j\left[t\right]=\left(2\mathsf{\pi }\Omega /60\right)t+2\mathsf{\pi }(j-1)/N$$

(16)

$$g\left({\varphi }_j\left[t\right]\right)=\left\{\begin{array}{l}1 \mathrm{if} {\varphi }_{st}<{\varphi }_j\left[t\right]<{\varphi }_{ex}\\ 0 \mathrm{otherwise}\end{array}\right.$$

(17)

where Ω is the spindle speed, and ϕ_st and ϕ_ex are the entry and exit angles on the j-th cutter tooth. For down milling, ϕ_st = arccos(2H-1) and ϕ_ex = π, and for up milling, ϕ_st = 0 and ϕ_ex = arccos(1-2H). H is the radial immersion ratio.

2.2.2. Estimating Environmental Excitation ê₁[t]

Substituting Equation (2) into Equation (5) yields the following equation

$$\begin{array}{l}A[z,t]x[t]=B[z,t](F[t]+e[t])\\ \iff G[z,t]x[t]=F[t]+e[t]\end{array}$$

(18)

where

$$G[z,t]=B{[z,t]}^{-1}\circ A[z,t]$$

(19)

Theoretically, the excitation can be represented by an infinite order inverse function model [30], namely, Equation (18) can be rewritten as

$${\displaystyle \sum _{i=0}^{\infty }{g}_i[t]}x[t-i]=F[t]+e[t]$$

(20)

To estimate e[t], the N_g items of the inverse function model are taken as

$$\begin{array}{l}{\displaystyle \sum _{i=0}^{N_g}{g}_i[t]}x[t-i]=F[t]+e[t]\\ \iff \left[g_0[t]g_1[t] \cdots g_{N_g}[t]\right]\left[\begin{array}{c}x[t]\\ x[t-1]\\ ⋮\\ x[t-N_g]\end{array}\right]=F[t]+e[t]\end{array}$$

(21)

with

$$\left\{\begin{array}{l}G{[t]}^{\mathrm{T}}=\left[g_0[t]g_1[t] \cdots g_{N_g}[t]\right]\\ \psi {[t]}^{\mathrm{T}}=\left[x{[t-1]}^{\mathrm{T}}x{[t-2]}^{\mathrm{T}}\cdots x{[t-N_g]}^{\mathrm{T}}\right]\end{array}\right.$$

(22)

Then, Substituting Equation (22) into Equation (21), Equation (21) can be reduced to

$$G^T[t]\psi [t]=F[t]+e[t]$$

(23)

where G^T[t] and e[t] are estimated using the least square method, and the estimation function of the least square method is given by

$$\min \left[\frac{1}{2}{\displaystyle \sum _{\tau =1}^t{‖G^{\mathrm{T}}[t]\psi [\tau ]-F[\tau ]‖}^2}\right]$$

(24)

where $‖•‖$ is the Euclidean norm. For this optimal problem, the solution of Equation (24) is calculated by

$${\widehat{G}}_1[t]={({\widehat{\psi }}_1[t]{\widehat{\psi }}_1{[t]}^{\mathrm{T}})}^{-1}{\widehat{\psi }}_1[t]{\widehat{F}}_1[t]$$

(25)

where

$$\left\{\begin{array}{l}{\widehat{\psi }}_1[t]=\left[\psi [1] \psi [2]\cdots \psi [t]\right] \\ {\widehat{F}}_1[t]={\left[F[1] F[2]\cdots F[t]\right]}^{\mathrm{T}}\end{array}\right.$$

(26)

Subsequently, substituting Equation (25) into Equation (23), the environmental excitation ê₁[t] can be estimated by

$${\widehat{e}}_1[t]={\widehat{G}}_1{[t]}^{\mathrm{T}}\psi [t]-F[t]$$

(27)

With the increase in iterative steps, the computational complexity of ${({\widehat{\psi }}_1[t]{\widehat{\psi }}_1{[t]}^{\mathrm{T}})}^{-1}$ progressively enlarges. To improve the computational efficiency of Equation (25), define the following

$$P_1[t]={({\widehat{\psi }}_1[t]{\widehat{\psi }}_1{[t]}^{\mathrm{T}})}^{-1}$$

(28)

Next, the following equation can be received

$$\begin{array}{l}{P}_1{[t]}^{-1}={\widehat{\psi }}_1[t]{\widehat{\psi }}_1{[t]}^{\mathrm{T}}={\widehat{\psi }}_1[t-1]{\widehat{\psi }}_1{[t-1]}^{\mathrm{T}}+\psi [t]\psi {[t]}^{\mathrm{T}}\\ =P_1{[t-1]}^{-1}+\psi [t]\psi {[t]}^{\mathrm{T}}\end{array}$$

(29)

Based on the matrix inverse theorem [34], P₁[t] can be rewritten as

$$P_1[t]=\left(P_1[t-1]-\frac{P_1[t-1]\psi [t]\psi {[t]}^{\mathrm{T}}{P}_1[t-1]}{I+\psi {[t]}^{\mathrm{T}}{P}_1[t-1]\psi [t]}\right)$$

(30)

Substituting Equation (30) into Equation (25), we can obtain

$$\begin{array}{l}{\widehat{G}}_1[t]=\left(P_1[t-1]-\frac{P_1[t-1]\psi [t]\psi {[t]}^{\mathrm{T}}{P}_1[t-1]}{I+\psi {[t]}^{\mathrm{T}}{P}_1[t-1]\psi [t]}\right)\left({\widehat{\psi }}_1[t-1]{\widehat{F}}_1[t-1]+\psi [t]F{[t]}^{\mathrm{T}}\right)\\ ={\widehat{G}}_1[t-1]+\frac{P_1[t-1]\psi [t]\left(F{[t]}^{\mathrm{T}}-\psi {[t]}^{\mathrm{T}}{\widehat{G}}_1[t-1]\right)}{I+\psi {[t]}^{\mathrm{T}}{P}_1[t-1]\psi [t]}\\ ={\widehat{G}}_1[t-1]+K_1[t]\end{array}$$

(31)

where

$$K_1[t]=\frac{P_1[t-1]\psi [t]\left(F{[t]}^{\mathrm{T}}-\psi {[t]}^{\mathrm{T}}{\widehat{G}}_1[t-1]\right)}{I+\psi {[t]}^{\mathrm{T}}{P}_1[t-1]\psi [t]}$$

(32)

From Equation (31), it is found that the relationship between Ĝ₁[t-1] and Ĝ₁[t] can be established. Therefore, the calculation volume of Equations (31) and (27) is reduced by iterative calculations.

2.2.3. Estimating Coefficient Matrices ŵ₁[t]

Substituting Equation (27) into Equation (3), the time-varying TARMAX dynamic equation can be expressed as

$$\begin{array}{l}x[t]+{\displaystyle \sum _{i=1}^{N_a}{a}_i[t]x[t-i]}={\displaystyle \sum _{j=0}^{N_b}{b}_j[t]({\widehat{e}}_1[t-j]+F[t-j])}\\ \iff x[t]=w^{\mathrm{T}}[t]\phi [t]\end{array}$$

(33)

where

$$\left\{\begin{array}{l}w[t]={\left[a_1[t]\cdots a_{N_a}[t]b_0[t]\cdots b_{N_b}[t]\right]}^{\mathrm{T}}\\ \phi {[t]}^{\mathrm{T}}=\left[-x{[t-1]}^{\mathrm{T}}\cdots -x{[t-N_a]}^{\mathrm{T}}F{[t]}^{\mathrm{T}}+{\widehat{e}}_1{[t]}^{\mathrm{T}}\cdots F{[t-N_b]}^{\mathrm{T}}+{\widehat{e}}_1{[t-N_b]}^{\mathrm{T}}\right]\end{array}\right.$$

(34)

Then, the coefficient matrices ŵ₁[t] can be estimated using the least square method, and

$${\widehat{w}}_1[t]={({\widehat{\phi }}_1[t]{\widehat{\phi }}_1{[t]}^{\mathrm{T}})}^{-1}{\widehat{\phi }}_1[t]{\widehat{X}}_1[t]$$

(35)

where

$$\left\{\begin{array}{l}{\widehat{\phi }}_1[t]=\left[\phi [1] \phi [2]\cdots \phi [t]\right]\\ {\widehat{X}}_1[t]={\left[x[1] x[2]\cdots x[t]\right]}^{\mathrm{T}}\end{array}\right.$$

(36)

Subsequently, combining Equations (34) and (35), the time-varying coefficient matrices a_i[t] and b_i[t] in the TARMAX model are obtained. Using the same simplified method in Equation (25) for Equation (35), defines the following

$$P_2[t]={({\widehat{\phi }}_1[t]{\widehat{\phi }}_1{[t]}^{\mathrm{T}})}^{-1}$$

(37)

Simplified in the same way as Equation (31), Equation (35) can be reduced to

$${\widehat{w}}_1[t]={\widehat{w}}_1[t-1]+K_2[t]$$

(38)

where

$$\begin{array}{l}{P}_2[t]=\left(P_2[t-1]-\frac{P_2[t-1]\phi [t]\phi {[t]}^{\mathrm{T}}{P}_2[t-1]}{I+\phi {[t]}^{\mathrm{T}}{P}_2[t-1]\phi [t]}\right)\\ K_2[t]=\frac{P_2[t-1]\phi [t]\left(x{[t]}^{\mathrm{T}}-\phi {[t]}^{\mathrm{T}}{\widehat{w}}_1[t-1]]\right)}{I+\phi {[t]}^{\mathrm{T}}{P}_2[t-1]\phi [t]}\end{array}$$

(39)

2.3. TARMAX Model Parameters Sliding Windows Recursive Estimation

According to Equations (26) and (36), it can be found that the dimensions of ${\widehat{\psi }}_1[t]$ and ${\widehat{\phi }}_1[t]$ gradually increase with time, and the amount of calculations for e[t] and a_i[t] will multiply. Therefore, to simplify the calculation process, a sliding window N is introduced. When t > N, the estimation of e[t] and a_i[t] only depends on the latest data. Then, e[t] and a_i[t] are calculated using the following method.

2.3.1. Estimating Environmental Excitation ê₂[t]

To improve the calculation efficiency, the latest data N are applied for the least square estimation function. Therefore, the estimation function of the sliding window least square method is used for ê₂[t] and is defined as

$$\min \left[\frac{1}{2}{\displaystyle \sum _{\tau =t-N+1}^t{‖G^{\mathrm{T}}[t]\psi [\tau ]-F[\tau ]‖}^2}\right]$$

(40)

For Equation (40), the solution is given by

$${\widehat{G}}_2[t]={({\widehat{\psi }}_2[t]{\widehat{\psi }}_2{[t]}^{\mathrm{T}})}^{-1}{\widehat{\psi }}_2[t]{\widehat{F}}_2[t]$$

(41)

where

$$\left\{\begin{array}{l}{\widehat{\psi }}_2[t]=\left[\psi [t-N+1] \psi [t-N+2]\cdots \psi [t]\right]\\ {\widehat{F}}_2[t]={\left[F[t-N+1] F[t-N+2]\cdots F[t]\right]}^{\mathrm{T}}\end{array}\right.$$

(42)

Subsequently, considering Equations (23) and (41), the environmental excitation ê₂[t] can be expressed by

$${\widehat{e}}_2[t]={\widehat{G}}_2{[t]}^{\mathrm{T}}\psi [t]-F[t]$$

(43)

To reduce the computational effort of Equation (41), use the following

$$\left\{\begin{array}{l}{P}_3{[t]}^{-1}={\widehat{\psi }}_2[t]{\widehat{\psi }}_2{[t]}^{\mathrm{T}}=P_3{[t-1]}^{-1}+U_1[t]U_1{[t]}^{\mathrm{T}}\\ U_1[t]=[\psi [t]j\psi [t-N]]\end{array}\right.$$

(44)

Similarly, considering matrix inverse theorem [34], P₃[t] can be converted to

$$P_3[t]=\left(P_3[t-1]-\frac{P_3[t-1]U_1[t]U_1{[t]}^{\mathrm{T}}{P}_3[t-1]}{I+U_1{[t]}^{\mathrm{T}}{P}_3[t-1]U_1[t]}\right)$$

(45)

Then, substituting Equation (45) into Equation (41) yields

$$\begin{array}{l}{\widehat{G}}_2[t]=\left(P_3[t-1]-\frac{P_3[t-1]U_1[t]U_1{[t]}^{\mathrm{T}}{P}_3[t-1]}{I+U_1{[t]}^{\mathrm{T}}{P}_3[t-1]U_1[t]}\right)\left({\widehat{\psi }}_2[t-1]{\widehat{F}}_2[t-1]+U_1[t]\overline{F}{[t]}^{\mathrm{T}}\right)\\ ={\widehat{G}}_2[t-1]+K_3[t]\end{array}$$

(46)

where

$$\left\{\begin{array}{l}{K}_3[t]=\frac{P_3[t-1]U_1[t]\left(\overline{F}{[t]}^{\mathrm{T}}-U_1{[t]}^{\mathrm{T}}{\widehat{G}}_2[t-1]\right)}{I+U_1{[t]}^{\mathrm{T}}{P}_3[t-1]U_1[t]}\\ \overline{F}[t]=\left[F[t] \mathrm{j}F[t-N]\right]\end{array}\right.$$

(47)

2.3.2. Estimating Coefficient Matrices ŵ₂[t]

Based on Equations (3) and (43), the estimation method in this section is similar to that of ŵ₁[t]. Therefore, the parameter ŵ₂[t] can be estimated using the least square method, and

$${\widehat{w}}_2[t]={({\widehat{\phi }}_2[t]{\widehat{\phi }}_2{[t]}^{\mathrm{T}})}^{-1}{\widehat{\phi }}_2[t]{\widehat{X}}_2[t]$$

(48)

where

$$\left\{\begin{array}{l}{\widehat{\phi }}_2[t]=\left[\phi [t-N+1] \phi [t-N+2]\cdots \phi [t]\right]\\ {\widehat{X}}_2[t]={\left[x[t-N+1] x[t-N+2]\cdots x[t]\right]}^{\mathrm{T}}\end{array}\right.$$

(49)

Defined

$$\left\{\begin{array}{l}{P}_4{[t]}^{-1}={\widehat{\phi }}_2[t]{\widehat{\phi }}_2{[t]}^{\mathrm{T}}=P_4{[t-1]}^{-1}+U_2[t]U_2{[t]}^{\mathrm{T}}\\ U_2[t]=\left[\phi [t] \mathrm{j}\phi [t-N]\right]\end{array}\right.$$

(50)

Similarly, Equation (48) can be reduced to

$${\widehat{w}}_2[t]={\widehat{w}}_2[t-1]+K_4[t]$$

(51)

where

$$\left\{\begin{array}{l}{P}_4[t]=\left(P_4[t-1]-\frac{P_4[t-1]U_2[t]U_2{[t]}^{\mathrm{T}}{P}_4[t-1]}{I+U_2{[t]}^{\mathrm{T}}{P}_4[t-1]U_2[t]}\right)\\ K_4[t]=\frac{P_4[t-1]U_2[t]\left(\overline{X}{[t]}^{\mathrm{T}}-U_2{[t]}^{\mathrm{T}}{\widehat{w}}_2[t-1]]\right)}{I+U_2{[t]}^{\mathrm{T}}{P}_4[t-1]U_2[t]}\\ \overline{X}[t]=\left[x[t] \mathrm{j}x[t-N]\right]\end{array}\right.$$

(52)

In addition, during calculation, false modes may appear, which may lead to the natural frequency estimation deviation. Therefore, to avoid false modes, a judging function J(·) is defined as

$$J({\widehat{f}}_r[t])=\left\{\begin{array}{l}{\widehat{f}}_r[t-1] \mathrm{if} \left|\left({\widehat{f}}_r[t]-{\widehat{f}}_r[t-1]\right)/{\widehat{f}}_r[t]\right|>\zeta \\ {\widehat{f}}_r[t] \mathrm{otherwise}\end{array}\right.$$

(53)

where ζ is the threshold value. If the relative error between the natural frequency at moment t and that at the moment t-1 is greater than the threshold value, then, the output natural frequency at moment t will be replaced by that at moment t-1.

From above calculation, when t < N, the environment excitation ê₁[t] can be determined by Equations (27) and (31), and the time-varying coefficient ŵ₁[t] can be obtained by Equation (38). Nevertheless, when t > N, the environment excitation ê₂[t] be determined by Equations (43) and (46), and the time-varying coefficient ŵ₂[t] can be obtained by Equation (51). In addition, the time-varying autoregressive coefficient a_i[t] can be extracted from ŵ₁[t] or ŵ₂[t] by Equation (34). Finally, the system modal parameters can be calculated by Equation (9), and Equations (11) and (53). Then, the flow chart of the proposed algorithm in milling is shown in Figure 2.

3. Numerical Simulation Verification

In this section, several numerical simulation experiments were used to verify the effectiveness of the proposed method. As the natural frequency of the first few orders of a thin-walled part has a large influence on the milling stability, the same three-degree-of-freedom structure as shown in the literature [18,37,38] was used to calculate the response of a thin-walled part subjected to a milling force excitation, which is shown in Figure 3. Then, the parameters were selected for numerical simulation and are shown in

Table 1. Simulation parameters for the 3-DOF time-varying structure [39].

Parameters	Values
Mass, stiffness and damping	m1 = 3 − 0.3t, m2 = 3 − 0.1t, m3 = 2 k1 = 5 × 107 − 107t, k2 = 3 × 107 − 5 × 106t, k3 = 107 − 106t, c1 = 200, c2 = 100, c3 = 50,
Variances of white noise signal	${\sigma }_1^2$ = 4, ${\sigma }_2^2$ = 2, ${\sigma }_3^2$= 1
Process parameters	Ω = 1500, m = 40, ϕex = 0, ϕts = π, n1 = 2, a = 2 × 10 − 4, ft = 10 − 3
Milling force coefficients	ktc1 = 6 × 107, krc1 = 2 × 107, kte1 = 3 × 103, kre1 = 103, Ktc2 = ktc3 = 8 × 107, kte2 = kte3 = 6 × 103, Krc2 = krc3 = 3 × 107, kre2 = kre3 = 4 × 103

. To analyze the system dynamic characteristics, three theoretical cutting forces were calculated, and 30 sets of white-noise signals were generated for simulation using the Monte Carlo method according to the parameters in Table 1. Subsequently, based on the “frozen time” method, the simulated time-varying frequency response functions and the natural frequencies of the 3-DOF time-varying structure were obtained, which are shown in Figure 4 and Figure 5, respectively.

From Figure 3 and Table 1, it can be seen that the system is excited by the independent milling force signals (theoretical milling force signal and Gaussian white noise signal). From Figure 4, the frequency response functions are obtained under excited conditions, but the natural frequency values of the system are approximately equal. In addition, from Figure 5, the first, the second, and the third natural frequencies are extracted, and we can see that the natural frequencies gradually reduced with time.

Subsequently, for the numerical simulation system, the response under a theoretical cutting force and Gaussian white noise was obtained in order for identifying the system modal parameters. In this process, 30 sets of displacement responses were obtained, then, to closely simulate the actual cutting condition, white noise signals with signal-to-noise ratios (SNR) of 15 dB, 20 dB, and 30 dB were added to contaminate the obtained displacement response signals; a set of displacement response signals with SNR = 15 dB is shown in Figure 6. Then, the cutting displacement responses contaminated with 15 dB, 20 dB, and 30 dB were used to calculate and determine the modal parameters of the simulated system.

To reduce the effect of initialization errors on the calculation results, the initial 500 samples (0.2 s) were discarded. Subsequently, the TARMAX model for the sliding window lengths of N = 50, N = 100, and N = 200 was used to estimate the system modal parameters, which were compared with the benchmark theoretical frequency and are shown in Figure 7.

From Figure 7, it was found that the estimated natural frequency presented good tracking results. As the SNR decreased, many scattered points deviated from the theoretical value. Under the constraint of Equation (53), the estimated natural frequencies fluctuated within the allowable error range, which validated the effectiveness of the proposed method. Then, when the sliding window length was relatively small, the natural frequencies were well estimated compared with the theoretical value, but fluctuated obviously. When the sliding window length was large, all of the estimated natural frequencies agreed well with the theoretical value, and fewer false modals appeared. In addition, it could be seen that the second-order natural frequency was not adequately tracked at SNR = 15 due to the large sliding window length, and the existence of some errors.

To further validate the proposed method, the mean absolute error (MAE) was introduced and is defined as

$$\mathrm{MAE}=\frac{1}{R}{\displaystyle \sum _{i=1}^R\frac{1}{L}{\displaystyle \sum _{t=1}^L\left|f[t]-{\widehat{f}}_i[t]\right|.}}$$

(54)

where f[t] is the theoretical natural frequency, ${\widehat{f}}_i[t]$ denotes the estimated natural frequency in the i-th Monte Carlo experiment, and L is the data length. Based on this, the different mean absolute error (MAE) was calculated and is shown in

Table 2. Modal parameter identification errors under different external excitations.

Noise	Mean Absolute Error/Hz
	N = 50	N = 100	N = 200
SNR = 30	18.9807	16.3046	10.2242
SNR = 20	22.2530	20.2492	16.0250
SNR = 15	26.1829	22.6237	19.0390

. From Table 2, with the increase in sliding window length, the prediction accuracy of the proposed method was obviously improved. Especially, when N = 200, the estimated accuracy was higher at the different SNR, and the MAE of the proposed method was also significantly smaller than that of the other sliding window lengths, which could indicate that the noise pollution was robust. Therefore, the proposed method could be widely used to identify the modal parameters under the data noise pollution, and be applied to practical conditions.

4. Experimental Validation and Discussion

In order to verify the effectiveness of the proposed method in the milling of thin-walled plates, several experiments were conducted. In the experiments, the materials of the plate and cutter were TC4 and cemented carbide, respectively, and the material parameters are shown in

Table 3. Workpiece and tool parameters.

Workpiece	Density	Poisson Ratio	Young’s Modulus	Materials
	4.6 g/cm3	0.34	108 GPa	TC4
Cutter	Diameter	Number of teeth	Spiral angle	Length
	12 mm	2	30°	75 mm

. The size of the plate used in the experiments was 80 × 40 × 3 mm. In addition, all of the experiments were carried out on a three-axis machine center (VMC-850E), the dynamic cutting forces in the milling were measured using a Kistler9257B, and the modal parameters were measured using a model hammer (500 N), an acquisition instrument DH5981, and acceleration sensors (ref. sensitivity 10.25 mV/g). The machining and modal test setups are shown in Figure 8.

4.1. Milling Force Coefficients Identification

To validate the effectiveness of the proposed method, the predicted cutting forces were needed. Therefore, the cutting force coefficients needed to be calibrated. Then, five slotting tests were carried out with a spindle speed of 1000 r/min and an axial depth of cut of 0.5 mm, while the feed rates were 40 mm/min, 80 mm/min, 120 mm/min, 160 mm/min, and 200 mm/min. Then, the milling force coefficient and the average milling force were obtained [36,39].

$$\left\{\begin{array}{l}{\overline{F}}_x=-\frac{na}{4}{k}_{rc}{f}_t-\frac{na}{\mathsf{\pi }}{k}_{re}\\ {\overline{F}}_y=\frac{na}{4}{k}_{tc}{f}_t+\frac{na}{\mathsf{\pi }}{k}_{te}\\ {\overline{F}}_z=\frac{na}{\mathsf{\pi }}{k}_{ac}{f}_t+\frac{na}{2}{k}_{ae}\end{array}\right.$$

(55)

Forces in the x and y directions were mainly considered in the model. In addition, the actually tested milling forces in the z-direction were small and easily contaminated by noise. Therefore, the cutting force coefficient in the z direction was not identified. Subsequently, the cutting force coefficients in the x and y directions were identified as follows k_tc = 1120.8 N/mm², k_rc = 2285.6 N/mm², k_te = 9.16 N/mm, and k_re = 13.21 N/mm.

4.2. Modal Parameters Evolution Considering Material Removal

To analyze the evolution of the modal parameters caused by the material removal, different impact experiments were conducted after five milling stages. In the experiments, the impact tests were conducted to assess the workpiece modal parameters before and after each machining stage. The experimental setups contained the acquisition apparatus DH5981, the accelerometer (5000 g, sensitivity 10.25 mV/g, and mass 5 g), and the modal hammer (500 N, sensitivity 1 mV/N). Then, the dynamic responses were changing at different impact positions and machining stages. Therefore, taking into account the fixture constraints, three points (point 1, point 2, and point 3) on a thin plate were selected for the tests, which are shown in Figure 9a, where it can be seen that point 2 is located in the middle of the edge of the plate, and point 1 and point 3 are symmetrical with respect to point 2. Subsequently, the frequency response functions on different impact measured points can be obtained under different machining stages, and it was found that the frequency response functions were almost the same for point 1 and point 3. However, considering the position of point 1 and point 3 on the plate, an unstable vibration signal caused by low rigidity would be generated. Therefore, point 2 was chosen for the measured response. Then, the sensor position and the machining region were determined and are shown in Figure 9b,c.

Next, five groups of milling experiments were conducted and the material removal patterns are shown in Figure 10. The selected cutting parameters were a spindle speed of 2000 rpm, axial depth of cut of 5 mm, radial depth of cut of 0.2 mm, feed rate of 120 mm/min, down milling, and the vibration responses in different cutting stages were obtained and are shown in Figure 11. After every cut, each impact test after milling was carried out, and the frequency response functions were obtained, which are shown in Figure 12. Then, from Figure 11, it can be seen that the acceleration response of the workpiece was unstable in the cut-in and cut-out stages due to the structural properties and rigidity. When the tool was fully cut into a workpiece, the coupled effect of workpiece–tool was weak, so the vibration was stable. As the number of cutting layers increased, the acceleration response was progressively reduced. The possible reasons for this are as follows: In the first cut, chatter occurred, which caused a larger acceleration response. However, as the material was removed, the modal parameters were time-varying in the milling and used the same machining parameters for the next milling, the chatter gradually decreased, which led to a low acceleration response. Subsequently, from Figure 12, we can find that with the workpiece material removal, the modal parameters of the machining system were significantly reduced from 525 Hz, 501 Hz, 504 Hz, 482 Hz, 465 Hz, to 453 Hz. Therefore, it is necessary to investigate the time-varying modal parameters for improving the stability.

4.3. Model Validation and Discussion

To identify the modal parameters of a machining system, a displacement response should be used for the proposed method. Therefore, a program by MATLAB was applied to convert the vibration acceleration responses measured in the milling to the vibration displacement responses, and the process is as follows.

Without a loss of generality, the acceleration signal $\ddot{x}\left(t\right)$ in the frequency domain can be expressed as follows

$$\ddot{x}\left(t\right)=A{e}^{\mathrm{j}\omega t}$$

(56)

where A is the coefficient corresponding to $\ddot{x}\left(t\right)$. Based on the inverse Fourier transform theory, assuming that the initial velocity and displacement of the system are 0, the displacement signal x[t] can be obtained using Equation (56).

$$x[t]={\displaystyle \sum _{k=0}^{N-1}-\frac{1}{{(2\mathsf{\pi }k\mathsf{\Delta }f)}^2}H[k]x[k{]\mathrm{e}}^{2\mathrm{j}\mathsf{\pi }k\mathrm{r}/N}}$$

(57)

where $X\left[k\right]=-A/{\omega }_{\mathrm{k}}^2$, ${\omega }_k=\mathrm{k}\Delta f$, H[k] denotes the cut-off frequency, which is defined as

$$H[k]=\left\{\begin{array}{l}1 (f_d\le k\Delta f\le f_u)\\ 0 (\mathrm{otherwise})\end{array}\right.$$

(58)

Based on Equations (56)–(58), the displacement response of the machining system after five cuts was calculated and is shown in Figure 13.

From Figure 11 and Figure 13, the vibration responses obviously fluctuated in the cut-in and cut-out stages, with the reason for this being that the dynamic cutting force was higher and the rigidity of the plate was lower in the cut-in and cut-out stages, respectively. For example, at the beginning and end of the fourth group of experiments, there was a significant change in the displacement signal. Therefore, to avoid anomalous data interference, the relatively stable response data over a 16 s period (from 3.3 s to 19.3 s) was chosen to identify the system modal parameters. Considering the actual milling process, the first order natural frequencies were closely associated with the system machining stability, which was validated by many references. In addition, as the modal parameters during machining were not available in real time, it was assumed in the paper that the natural frequency of the workpiece varied linearly during each experiment. Finally, the first order natural frequencies were calculated by TARMAX (N = 200) in different states, and the predicted results are shown in Figure 14.

From Figure 14, it was found that the predicted natural frequencies using the proposed method agreed well with the experimental values. The experimental average of the five tests was 513, 502.5, 493, 473.5, and 459 Hz and the predicted average of the five tests was 516.16, 506.24, 496.45, 479.31, and 462.71 Hz. The predicted values were all higher than the experimental values due to the coupling benefits of the tool and the workpiece during machining, which is consistent with the experimental phenomena in the literature [15]. In addition, several scatter points offset the baseline under the machining start and end stage. The reason for this is that the milling was at the beginning or end of cut at this time, and the vibration displacement response obviously changed, which caused a large deviation between the predicted values and the experimental values.

To further demonstrate the validity of the proposed method, the mean absolute error (MAE) and mean relative error (MRE) were introduced, and are defined as

$$\left\{\begin{array}{l}\mathrm{MAE}=\frac{1}{L}{\displaystyle \sum _{t=1}^L\left|f[t]-\widehat{f}[t]\right|}\\ \mathrm{MRE}=\frac{1}{L}{\displaystyle \sum _{t=1}^L\frac{\left|f[t]-\widehat{f}[t]\right|}{f[t]}}\end{array}\right.$$

(59)

where f[t] is the theoretical natural frequency, $\widehat{f}[t]$ denotes the estimated natural frequency, and L is the data length. Based on this, MAE and MRE were calculated and shown in Figure 15.

From Figure 15, the average absolute error of the TARMAX method was around 20 Hz and the relative error was around 4.3%, and the prediction results were relatively stable. The third group experiments had the highest prediction accuracy because smoother displacement signals were chosen. Next, there were some non-smooth signals in the displacement signals in the other sets of experiments. As the number of milling experiments increased, the relative error of the TARMAX method showed an upward trend. The possible reasons for this are as follows. As the number of experiments increased, the tool wore out and the dynamic milling forces changed, thus causing fluctuations in the vibration response, which is consistent with the observation that the displacement signals in the fourth and fifth experiment included some non-stationary components. In addition, the calculation time and complexity of the TARMAX method are mainly relevant to the sliding window length, which makes the proposed method handle vibration responses online better. Overall, the TARMAX method can commendably predict the modal parameters in the milling of a thin-walled workpiece.

5. Conclusions

In the milling of a thin-walled workpiece, the dynamic characteristics of the thin-walled workpiece-fixture system are time-varying due to material removal. To investigate the time-varying characteristics of the modal parameters in milling, an accurate TARMAX identification method under external excitation was proposed. The proposed method can rapidly and precisely identify the modal parameters in the milling of thin-walled workpieces. Thus, the contributions of the paper are as follows.

(1)

The TARMAX modal parameter identification method under external excitation was proposed. In this process, the dynamic modal parameters are estimated through the measured vibration response in milling using the recursive estimation method or sliding windows recursive estimation.

(2)

A 3-DOF time-varying structure model under milling forces and white noise excitation was established. The predicted natural frequencies using the proposed method were in close agreement with the simulated values.

(3)

The proposed model is validated by the numerical simulation and machining experiments. Moreover, the prediction accuracy of the TARMAX identification method was 95.68%, and good agreement was drawn in the numerical simulation and machining experiments. Our further research objectives will focus on the prediction of damping ratios, modal stiffness, and modal mass.