Systematic study of the prediction methods for machined surface topography and form error during milling process with flat-end cutter

Abstract

An enhanced time-domain simulation method of cutter/workpiece dynamic engagement during milling process is proposed in this article, which comprehensively considers the effect of multi-order modal characteristics of cutter system and cutter runout including offset and inclination. Based on the cutter dynamic displacement response, this article further presents the machined surface reconstruction algorithm and evaluation method for the form error. This research systematically studies and compares the calculation accuracy between the proposed method and other previous three kinds of methods. The effectiveness of the proposed method has been verified by a series of milling experiments successfully. By comparing with the other three methods, the proposed method shows a high calculation performance, especially under the milling condition with a large axial depth of cut and low damping or stiffness of cutter system. Besides, the results indicate that the form error has a strong dependence characteristic on the milling parameters, particularly on spindle speed. Additionally, cutter runout would easily cause over cut phenomenon on the machined surface and seriously deteriorate the surface roughness.

Introduction

With the continuous improvement of the part quality requirement, how to obtain a high precision machining has become an important challenge for the manufacturing process. In order to improve the machining quality and save the manufacturing cost, virtual simulation technology ¹ of machining process has been gradually applied to study the relationship between the machining parameters and part quality, and it further helps manufacturer to choose reasonable process parameters. At present, as one of the main methods of metal cutting, the machining quality in milling process has been studied by a lot of researchers.

During milling process, the tool/workpiece engagement process is the superimposed result of the macroscopic relative motion between the cutter and workpiece and the micro-amplitude vibration induced by milling force, which leads to the actual machined surface deviating from its desired position and thus producing the form error. In order to study the form error, it is necessary to reconstruct the machined surface based on the vibration response of cutter and workpiece through theoretical or experiment method. For the form error during milling process, there are mainly four kinds of different methods. The first method is the cutter/workpiece geometric engagement method, which ignores the system dynamic excitation and response caused by milling force and has been mainly used to analyze the surface roughness. Based on the geometric relative motion between cutter and workpiece, Gao et al. ² presented a numerical simulation method for machined surface topography during milling process with flat- and ball-end cutter. Arizmendi et al.^3,4 studied the effect of cutter runout on surface roughness in peripheral milling, especially analyzed the position and width of heterogeneity band and the characteristic of roughness distribution along cutter axis direction. Corral et al. ⁵ discussed the influence of feed, eccentricity and helix angle on machined surface topography in side milling processes. Later, Corral et al. ⁶ applied the same method to study the surface reconstruction algorithm with ball-end cutter, and emphatically researched the effect of feed per tooth and radial depth of cut on surface roughness, and gave out the strategy for improving the machining efficiency.

The second one is the static deflection method, which considers the static excitation effect by milling force. In early studies, Kline et al. ⁷ analyzed the influence of cutter static deflection on machining accuracy during milling process based on the mechanistic cutting force model. DeVor and Sutherland ⁸ further improved this method, in which the feedback effect of cutter deflection on cutting load is considered. Altintas and Lee ⁹ applied the cutter static deflection method and concluded that the variable-pitch helical end mills can increase the machining accuracy of machined surface by 10%–20%. Law and Geddam ¹⁰ studied the form error in straight line milling and corner milling during rough and finishing operations. Yun et al. ¹¹ constructed the three-dimensional (3D) surface error map by obtaining the cutting edge position at which the flutes pass over the machined surface. Ryu et al.^12,13 built a form error calculation model and discussed the variation of error topography under different milling parameters and milling types. Dedai and Rao ¹⁴ presented a classification method of form error profile corresponding to milling parameters, which was further to be used for error prediction during curved surface milling. During ball-end milling of inclined convex and concave surfaces, Gök et al. ¹⁵ calculated the tool deflections depending on cutting force acting on the cantilever beam model of cutter. For workpiece static deflection, Ratchev et al. ¹⁶ reported on a new simulation methodology of static deflection of thin-wall parts during machining based on finite element method (FEM), which was further used for compensation of poor geometric accuracy. Bolar and Joshi ¹⁷ proposed a 3D workpiece static deflection method, which provided very useful insights into the complex physical interaction of cutter and workpiece.

Different from the first two methods, the third one is the harmonic balance method based response directly solving method, which considers the dynamic response of cutter and workpiece induced by milling force. Hence, the advantage of this method is that the dynamic excitation characteristics of milling force with harmonic components depended on spindle speed are considered. In this method, it is necessary to have a prior-knowledge about the milling force. Montgomery and Altintas ¹⁸ early proposed the simulation and analysis method for machining surface topography based on milling dynamics. Then, the influence of cutter eccentricity and system dynamic response on the surface error morphology is further analyzed in their studies. ¹⁹ Ismail et al. ²⁰ put forward a generation model of machined surface, in which the cutter runout, vibration and wear were comprehensively considered. Schmitz and colleagues^21,22 gave out a closed-form solution for cutter displacement response and pointed out that the cutter runout would affect the roughness and surface location error. Yuan et al. ²³ analyzed the surface topography and surface location error based on the harmonic balance method, in which the influence of feed per tooth, runout and number of cutter teeth on topography was discussed. Denkena et al. ²⁴ reconstructed the cutter actual motion and extracted the form surface morphology based on the measured milling force data and the frequency coupling dynamic model of tool system. Besides, Jiang et al., ²⁵ Arizmendi et al. ²⁶ and Moreau and Costes ²⁷ researched the machined surface and form error through the measured cutter vibration response.

The last one is the time-domain simulation method based on cutter/workpiece dynamic engagement process. As the cutter/workpiece engagement process exist the force–displacement coupling effect, it is necessary to consider this kind of coupling effect for obtaining high simulation accuracy. With the aid of time-domain simulation, Zhang et al. ²⁸ presented an integrated model of feed drive system with milling process and pointed out that the displacement fluctuation has an important effect on the contour roughness. Surmann and Enk ²⁹ and Surmann and Biermann ³⁰ predicted the transient behavior of cutter vibration along changing engagement conditions, which can allow the feedback effect of cutter displacement into the cutter/workpiece engagement and thus on chip form, and then constructed a geometric model of machined surface with the knowledge of tooth vibration trajectory. Schmitz ³¹ focused on the time-domain simulation method and applied the surface location error (SLE) to evaluate the machining quality.

The above literature views show that a great improvement about form error during milling process has been made, but what should be pointed out is that there also exist some disadvantages in these methods. First, all the four kinds of methods can be used to study the machined surface topography and form error, but the application ranges of these methods are not given out, especially the comparison of superiority among them is ignored. Besides, as a flexible system, the cutter system, being not a single-degree-of-freedom (DOF) dynamics model, includes multi-order modal characteristics in fact. Especially under the milling condition with a large axial depth of cut, the multi-order modal characteristics will become prominent, which should be considered. Finally, cutter runout is an important factor to affect the machining quality, but its influence cannot be taken into account effectively. Hence, this article presents an enhanced time-domain simulation method for cutter/workpiece engagement process with consideration of the multi-DOF cutter dynamics model and cutter runout and systematically compares the calculation accuracy between the proposed method and other previous methods. The contents of this article are organized as following. In section “Simulation method of machined surface topography during milling process”, this article gives out the simulation method of machined surface topography during milling process, which includes the proposed method and other three kinds of methods, and then presents the reconstruction algorithm of machined surface. In section “Simulation and experiment verification”, a series of milling tests are conducted to verify the effectiveness of the proposed method. In section “The influence of different milling conditions on form error”, the influences of milling parameters, cutter runout, milling types, tooth helix angle, damping and stiffness of cutter system and multi-DOF cutter dynamics model on form error are discussed, respectively.

Simulation method of machined surface topography during milling process

Definition of cutter runout

Figure 1 is the parametric description model of cutter runout. ${\mathbf{X}}_c{\mathbf{Y}}_c{\mathbf{Z}}_c$ is the cutter geometry coordinate system, in which the tooth edge shape is represented. ${\mathbf{X}}_r{\mathbf{Y}}_r{\mathbf{Z}}_r$ is the cutter rotation coordinate system, and the existence of cutter runout will lead to the geometry center of cutter deviating from its rotation center. Generally, cutter runout includes two cases: offset and inclination, where cutter offset can be described by the offset value $\rho$ and offset angle $\lambda$ , and cutter inclination can be described by the inclination value $\gamma$ and inclination angle $\eta$ . In order to analyze the tool/workpiece engagement process, the cutter will be discretized into $N_a$ cutting disk elements along the cutter axis direction with equal thickness within the axial depth of cut a_p, and the thickness of each disk element is db.

Figure 1.

General definition of cutter runout.

In the cutter geometry coordinate system, the geometry shape of tooth cutting edge of the jth cutting disk element on the ith tooth is expressed as

$$\begin{array}{l}\left[\begin{array}{c}{x}_c(i,j)\\ y_c(i,j)\\ z_c(i,j)\end{array}\right]=\left[\begin{array}{c}-R_0\sin [2\pi (i-1)/N_t+(j-1){\beta }_0]\\ R_0\cos [2\pi (i-1)/N_t+(j-1){\beta }_0]\\ (j-1)db\end{array}\right]\end{array}$$

(1)

where $R_0$ is the nominal cutter radius. $N_t$ is the tooth number. ${\beta }_0=\mathrm{db}\tan \alpha /R_0$ is the lag angle of each disk element due to the helix angle, in which $\alpha$ is the helix angle.

Then, in the cutter rotation coordinate system, the tooth cutting edge can be derived as

$$\begin{array}{l}\left[\begin{array}{c}{x}_r(i,j)\\ y_r(i,j)\\ z_r(i,j)\end{array}\right]=\left[\begin{array}{ccc}{\cos }^2\eta (1-\cos \gamma )+\cos \gamma & -\sin \eta \cos \eta (1-\cos \gamma )& -\sin \eta \sin \gamma \\ -\sin \eta \cos \eta (1-\cos \gamma )& {\sin }^2\eta (1-\cos \gamma )+\cos \gamma & -\cos \eta \sin \gamma \\ \sin \eta \sin \gamma & \cos \eta \sin \gamma & \cos \gamma \end{array}\right]\\ \left[\begin{array}{c}{x}_c+\rho \sin \lambda \\ y_c+\rho \cos \lambda \\ z_c-L\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ L\end{array}\right]\end{array}$$

(2)

where $L$ is the cutter overhang length. Finally, as a result of cutter runout, the actual cutting radius is given by

$$R_{i,j}=\sqrt{x_r{(i,j)}^2+y_r{(i,j)}^2}$$

(3)

Solving cutter displacement response and modeling tooth trajectory

This subsection mainly proposes a time-domain simulation method for cutter/workpiece dynamic engagement process with consideration of multi-DOF modal characteristics of cutter system and cutter runout, solves the cutter displacement and builds up the tooth trajectory model. Then, other three kinds of modeling methods for the tooth trajectory are also given briefly.

Method 1: cutter/workpiece dynamic engagement–based time-domain simulation method

In the proposed method, as shown in Figure 2(a), the milling process assumes that the cutter is flexible and the workpiece is rigid. The machined surface will deviate from the desired surface resulting from the excitation of milling force on cutter system. In order to analyze the cutter vibration, the cutter system is assumed as a multi-DOF dynamics system with a series of lumped mass nodes shown as the green nodes in Figure 2(a), and the cutter/workpiece engagement process of the jth disk element on section A in Figure 2(a) is further detailed in Figure 2(b). Here, $\mathrm{XOY}$ is the global coordinate system. ${\mathrm{X}}_{\mathrm{sc}}{\mathrm{O}}_{\mathrm{sc}}{\mathrm{Y}}_{\mathrm{sc}}$ is the cutter static coordinate system, which moves along the feed direction without cutter vibration. ${\mathrm{X}}_{\mathrm{dc},j}{\mathrm{O}}_{\mathrm{dc},j}{\mathrm{Y}}_{\mathrm{dc},j}$ is the cutter dynamic coordinate system of the jth cutting disk element, which is fixed at the cutter rotation center.

Figure 2.

(a) Milling process and (b) cutter/workpiece dynamic engagement.

In physical space, the dynamic equation of the flexible cutter system is

$$M_c{\stackrel{··}{x}}_c+C_c{\stackrel{.}{x}}_c+K_c{x}_c=F_c^N$$

(4)

where $x_c=[O_{\mathrm{dc},1,x},\dots ,O_{\mathrm{dc},N_s,x},O_{\mathrm{dc},1,y},\dots ,O_{\mathrm{dc},N_s,y}{]}^T$ is the displacement vector of cutter nodes, in which $N_s$ is the number of nodes. ${\mathbf{M}}_c$ is the mass matrix. ${\mathbf{C}}_c$ is the damping matrix. ${\mathbf{K}}_c$ is the stiffness matrix. $F_c^N=[F_{1,x},\dots ,F_{N_s,x},F_{1,y},\dots ,F_{N_s,y}{]}^T$ is the force vector acting on cutter nodes.

Based on modal analysis theory, in modal space, the dynamic equation of the flexible cutter system can be given by

$${\mathbf{I}}_c^m{\stackrel{··}{\zeta }}_c+{\mathbf{C}}_c^m{\stackrel{·}{\zeta }}_c+{\mathbf{K}}_c^m{\zeta }_c={\mathbf{U}}_c^T{F}_c^N$$

(5)

where ${\zeta }_c=[O_{\mathrm{dc},1,x}^m,\dots ,O_{\mathrm{dc},N_s,x}^m,O_{\mathrm{dc},1,y}^m,\dots ,O_{\mathrm{dc},N_s,y}^m{]}^T$ is the displacement vector of cutter in modal space. ${\mathbf{I}}_c^m={\mathbf{U}}_c^T{\mathbf{M}}_c{\mathbf{U}}_c$ is the mass-normalized matrix, inwhich ${\mathbf{U}}_c$ is the modal shape matrix and satisfies $x_c={\mathbf{U}}_c{\zeta }_c$ . ${\mathbf{C}}_c^m={\mathbf{U}}_c^T{\mathbf{C}}_c{\mathbf{U}}_c=\mathrm{diag}[2{\xi }_{x1}{\omega }_{\mathrm{nx}1},\cdots ,2{\xi }_{x{N}_s}{\omega }_{\mathrm{nx}{N}_s},2{\xi }_{y1}{\omega }_{\mathrm{ny}1},\cdots ,2{\xi }_{y{N}_s}{\omega }_{\mathrm{ny}{N}_s}{]}^T$ is the damping matrix in modal space, in which $\xi$ and ${\omega }_n$ are the damping ratio and natural frequency, respectively. ${\mathbf{K}}_c^m={\mathbf{U}}_c^T{\mathbf{K}}_c{\mathbf{U}}_c=\mathrm{diag}[{{\omega }_{\mathrm{nx}1}}^2,\dots ,{{\omega }_{\mathrm{nx}{N}_s}}^2,{{\omega }_{\mathrm{ny}1}}^2,\dots ,{{\omega }_{\mathrm{ny}{N}_s}}^2{]}^T$ is the stiffness matrix in modal space.

During the cutter/workpiece engagement process, there exists a coupling effect between milling force and cutter displacement, and the instantaneous uncut chip thickness is the most important parameter for presenting this kind of coupling effect. In the proposed time-domain simulation method, first, the continuous time is discretized into a time sequence with an interval of $\mathrm{\Delta }t$ . Then, in current time interval, the instantaneous uncut chip thickness can be calculated based on the input of previous cutter displacement and milling force. Finally, the current cutter displacement can be obtained through the dynamic equation of the flexible cutter system. Repeating the above iteration process, we can acquire the cutter displacement response during the total machining time. At current time $t_k$ , the general tooth trajectory equation of cutting point $\mathbf{S}$ in the global coordinate system can be expressed as

$$\left[\begin{array}{c}{S}_{x,i,j,k}\\ S_{y,i,j,k}\end{array}\right]=\left[\begin{array}{c}{O}_{sc,x}(t_k)\\ O_{sc,y}(t_k)\end{array}\right]+\left[\begin{array}{c}{O}_{dc,j,x}(t_k)\\ O_{dc,j,y}(t_k)\end{array}\right]+\left[\begin{array}{c}{R}_{i,j}\sin ({\varphi }_{i,j}(t_k))\\ R_{i,j}\cos ({\varphi }_{i,j}(t_k))\end{array}\right]$$

(6)

where $[O_{\mathrm{sc},x}{O}_{\mathrm{sc},y}{]}^T$ is the cutter static displacement at time $t_k$ , in which the component in X-direction is equal to $f_v{t}_k$ and $f_v$ is the feed speed. $[O_{\mathrm{dc},j,x}{O}_{\mathrm{dc},j,y}{]}^T$ is the cutter dynamic displacement of the jth disk element at time $t_k$ . As the number of disk elements is larger than that of DOF, the cutter dynamic displacement vector can be obtained by interpolating based on the displacement vector of cutter $x_c$ . ${\varphi }_{i,j}(t_k)=\pi n{t}_k/30-{\gamma }_{i,j}-(j-1){\beta }_0$ is the rotation angle of the jth disk element on the ith tooth, in which $n$ is the spindle speed, and ${\gamma }_{i,j}$ is the actual included angle between the ith and first tooth.

At current time, the instantaneous uncut chip thickness $h(i,j,t_k)$ is equal to the distance between point $\mathbf{S}$ and $\mathbf{T}$ . Due to the cutter runout, the point $\mathbf{T}$ maybe on the machined surface left behind by previous m_ith tooth. Besides, point $\mathbf{T}$ is not only on straight line $\overline{{\mathbf{O}}_{\mathrm{dc}}\mathbf{S}}$ but also on the straight line $\overline{S_{i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}}{S}_{i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}+1}}$ , where $\mathrm{\Delta }{k}_{m_i}$ is the sampling number of interval between the current tooth and previous m_ith tooth with the same rotation angle, and it can be obtained by iterative searching. Then, point $\mathbf{T}$ can be acquired by solving the variables $u,v(0⩽u,v⩽1)$ in the following equation

$$\{\begin{array}{c}{S}_{x,i,j,k}+u({\mathbf{O}}_{\mathrm{dc},x}-S_{x,i,j,k})=S_{x,i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}}+u(S_{x,i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}+1}-S_{x,i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}})\hfill \\ \multicolumn{1}{c}{S_{y,i,j,k}+v({\mathbf{O}}_{\mathrm{dc},y}-S_{y,i,j,k})=S_{y,i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}}+v(S_{y,i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}+1}-S_{y,i-m_i,j,k-\mathrm{\Delta }{k}_{m_i}})(m_i=1,2,\dots ,N_t)}\end{array}$$

(7)

Furthermore, calculating the distance between point $\mathbf{S}$ and point $\mathbf{T}$ , the possible instantaneous uncut chip thickness can be obtained as $h_{i,j}(t_k,m_i)=\overline{\mathrm{ST}}$ . For the cutting point $\mathbf{S}$ at current time, the actual instantaneous uncut chip thickness is the maximum between the minimum of all possible instantaneous uncut chip thickness and zero, which can be mathematically expressed by

$$\begin{array}{l}{h}_{i,j}(t_k)=\max \left(\min \left(h_{i,j}(t_k,m_i)\right),0\right)W({\varphi }_{i,j}(t_k)),\\ m_i=1,2,\mathrm{\dots },N_t\end{array}$$

(8)

where $W$ is the window function. Considering the dynamic engagement process, $W$ of down milling and up milling can be derived as

$$\begin{array}{l}\mathrm{Down\; milling}:\\ W({\varphi }_{i,j}(t_k))=\{\begin{array}{cc}1& S_{y,i,j,k}-O_{sc,y}⩽-R_0+a_e\\ 0& S_{y,i,j,k}-O_{sc,y}>-R_0+a_e\end{array}\end{array}$$

(9)

$$\begin{array}{l}\mathrm{Up\; milling}:\\ W({\varphi }_{i,j}(t_k))=\{\begin{array}{cc}1& S_{y,i,j,k}-O_{sc,y}⩾R_0-a_e\\ 0& S_{y,i,j,k}-O_{sc,y}<R_0-a_e\end{array}\end{array}$$

(10)

Combining the mechanistic cutting force model and the calculated instantaneous uncut chip thickness, we can further achieve the milling force of each cutting disk element of cutter at current time

$$\left[\begin{array}{c}{F}_{c,j,x}(t_k)\hfill \\ \multicolumn{1}{c}{F_{c,j,y}(t_k)}\end{array}\right]=\sum _i\mathrm{Rot}[\left\{\begin{array}{c}{k}_{\mathrm{ts}}\hfill \\ \multicolumn{1}{c}{k_{\mathrm{rs}}}\end{array}\right\}h(i,j,t_k)+\left\{\begin{array}{c}{k}_{\mathrm{tp}}\hfill \\ \multicolumn{1}{c}{k_{\mathrm{rp}}}\end{array}\right\}]\mathrm{db}$$

(11)

where $k_{\mathrm{qs}}$ and $k_{\mathrm{qp}}(q=t,r)$ are the shearing and plowing specific cutting force coefficient in tangential and radial direction, respectively. $\mathrm{Rot}$ is the rotation matrix. At this moment, the milling force acting on cutter is a kind of distributed force, and we should transform it into concentrated force for each DOF nodes of cutter system. The concentrated force acting on the sth node of cutter system in X and Y directions can be expressed as

$$F_{c,s,q}^N(t_k)=\sum _{j=L_s}^{j=U_s}{F}_{c,j,q}(t_k)(q=x,y)$$

(12)

where $L_s$ and $U_s$ are the upper and lower limit number of disk elements for the sth node of cutter system. Let $Y=[x_c,{\stackrel{·}{x}}_c{]}^T$ , the state space function of cutter dynamic system is

$$\stackrel{·}{Y}=\mathbf{A}Y+\mathbf{B}{F}_c^N$$

(13)

where $\mathbf{A}=\left[\begin{array}{cc}0\hfill & {\mathbf{I}}_c\hfill \\ \multicolumn{1}{c}{-{\mathbf{M}}_c^{-1}{\mathbf{K}}_c}& -{\mathbf{M}}_c^{-1}{\mathbf{C}}_c\hfill \end{array}\right]$ and $\mathbf{B}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{{\mathbf{M}}_c^{-1}}\end{array}\right]$ .

With the aid of the response $Y(t_k)$ at current time $t_k$ , we can further calculate the system response $Y(t_{k+1})$ at next time $t_{k+1}$ using the following equation

$$Y(t_{k+1})=e^{\mathbf{A}\mathrm{\Delta }t}Y(t_k)+\underset{t_k}{\overset{t_{k+1}}{\int }}{e}^{\mathbf{A}(t_{k+1}-\tau )}\mathbf{B}{F}_c^N(t_k)d\tau$$

(14)

After obtaining the cutter displacement response at time $t_{k+1}$ , we can set it as the initial condition for the next time interval. Repeating the above iteration process within the total machining time, the cutter displacement response during the milling process can be obtained and then the tooth trajectory can be also calculated by equation (6).

Method 2: harmonic balance method based response directly solving method

Unlike the method 1, method 2 calculates the cutter displacement response using the harmonic balance method with the aid of pre-calculated milling force and cutter dynamics. In this method, the dynamic engagement process between cutter and workpiece is ignored.

Considering the cutter runout and tooth trochoid trajectory, the possible instantaneous uncut chip thickness during cutter/workpiece geometric engagement can be shown as equation (15), which is detailed in Zhang et al. ³²

$$\begin{array}{l}{h}_{i,j}(t,m_i)=\sin ({\varphi }_{i,j}(t))f_v{\tau }_{i,j}(t,m_i)+R_{i,j}\\ -R_{i-m_i,j}\cos \left(\frac{2\pi n{\tau }_{i,j}(t,m_i)}{60}-{\gamma }_{i,j}+{\gamma }_{i-m_i,j}\right)\end{array}$$

(15)

where ${\tau }_{i,j}(t,m_i)$ is the time delay term.

Combining equations (8)–(12), we can also obtain the milling force. In the harmonic balance method, the milling force acting on DOF nodes of cutter is expressed as a Fourier series expansion

$$F_{c,s,q}^N=\sum _{k=-\infty }^{\infty }{A}_k{e}^{\mathrm{ik}(\omega t+{\theta }_k)}(q=x,y)$$

(16)

where ${\mathbf{A}}_k=[A_{k,1,x},\dots ,A_{k,N_s,x},A_{k,1,y},\dots ,A_{k,N_s,y}{]}^T$ and ${\theta }_k=[{\theta }_{k,1,x},\dots ,{\theta }_{k,N_s,x},{\theta }_{k,1,y},\dots ,{\theta }_{k,N_s,y}{]}^T$ . Then, the cutter displacement response is

$$x_c=\sum _{k=-\infty }^{\infty }{\mathbf{B}}_k{e}^{\mathrm{ik}\omega t}$$

(17)

where ${\mathbf{B}}_k=[B_{k,1,x},\dots ,B_{k,N_s,x},B_{k,1,y},\dots ,B_{k,N_s,y}{]}^T$

Substituting equations (16) and (17) into equation (4), we can further obtain

$$\begin{array}{l}{\displaystyle \sum _{k=-\infty }^{\infty }\left(-{(k\omega )}^2{M}_c+ik\omega C_c+K_c\right)B_k}\\ ={\displaystyle \sum _{k=-\infty }^{\infty }{A}_k{e}^{{\theta }_k}}(q=x,y)\end{array}$$

(18)

Then, ${\mathbf{B}}_k$ can be calculated by

$${\mathbf{B}}_k=\mathbf{H}(k\omega ){\mathbf{A}}_k{e}^{{\theta }_k}(q=x,y)$$

(19)

where $\mathbf{H}(k\omega )=(-(k\omega {)}^2{\mathbf{M}}_c+\mathrm{ik}\omega {\mathbf{C}}_c+{\mathbf{K}}_c{)}^{-1}$ is the frequency response function of cutter dynamic system.

Finally, the cutter displacement response can be acquired by substituting equation (19) into equation (17), and the tooth trajectory can be also calculated with the aid of equation (6).

Method 3: cutter static deflection method

In this method, the cutter displacement is only calculated based on its static deflection, and the dynamic response of cutter system is not considered in this model. Here, the calculation procedure of milling force is the same with method 2. Then, the static deflection of cutter can be shown as equation (20), and the tooth trajectory can be also obtained

$$x_c(t_k)={{\mathbf{K}}_c}^{-1}{F}_c^N(t_k)$$

(20)

Method 4: cutter/workpiece geometric engagement method

For method 4, the dynamic excitation of milling force for cutter system is not considered, and the cutter/workpiece engagement process is only based on the geometric constraint between them. At this moment, the dynamic response of cutter system is zero. Hence, let $O_{\mathrm{dc},s,q}=0(q=x,y)$ in equation (6), we can directly obtain the tooth trajectory.

Machined surface reconstruction algorithm

The machined surface is the result of the intersection between tooth trajectory surface and the workpiece. As shown in Figure 3(a), the actual machined surface $S_{u,w}^m$ of a milling process with a four-teeth cutter is a connected result of all tooth trajectory surfaces, where $u$ and $w$ are the sampling numbers in X and Z directions. As the existence of cutter runout, the tooth trajectory surfaces for all teeth are not the same with each other. Besides, as a result of dynamic excitation by milling force, the actual milling machined surface will deviate from the desired surface $S_{u,w}^d$ , and the deviation in the normal direction is the form error.

Figure 3.

(a) Actual machined surface and desired surface and (b) machined surface reconstruction.

In order to analyze the form error, the actual machined surface should be described in the global coordinate system, and the sampling point on it should also be corresponded to the regular points on the desired surface $S_{u*,w*}^m$ . Taking the point $S_{u*,w*}^d$ on the desired surface for example, its corresponding point on the machined surface is $S_{u*,w*}^m$ , which coincides with the point $S_{i*,j*,k*}$ on tooth trajectory surface and is on the normal line of desired surface through the point $S_{u*,w*}^d$ . Here, the coordinate component in Z direction of point $S_{u*,w*}^m$ and point $S_{u*,w*}^d$ is the same with each other, that is, $S_{z,u*,w*}^m=S_{z,u*,w*}^d$ and $j*=S_{z,u*,w*}^m/\mathrm{db}$ . Hence, it is only needed to solve the coordinate components in X and Y directions of point $S_{u*,w*}^m$ , that is to say, we should only find out the variables $i*$ and $k*$ .

As shown in Figure 3(b), the point $S_{i*,j*,k*}$ is not only on the machined surface but also on the tooth trajectory surface. For point $S_{u*,w*}^d$ , each tooth trajectory has a corresponding point $S_{i,j,k_i}(i=1,2,\dots ,N_t)$ , but only the point with the minimum distance from the ideal point is the point $S_{i*,j*,k*}$ on the machined surface. The solving procedure for the variables $i*$ and $k*$ of point $S_{i*,j*,k*}$ is given as

• Step 1. For a given point $S_{u*,w*}^d$ , extract two period data of tooth trajectory and storage them in $S_{i,j,k}^0$

$$S_{i,j,k}^0=S_{i,j,k}(k_p⩽k⩽k_p+2T/\mathrm{\Delta }t)$$

(21)

where $k_p$ is the sampling number, which satisfies $S_{x,1,j,k_p}^0⩽S_{x,u*,w*}^d⩽S_{x,1,j,k_p+2T/\mathrm{\Delta }t}^0$ .

• Step 2. Determine the point $S_{i,j,k_i}^0$ for each tooth, the distance of which and point $S_{u*,w*}^d$ has a minimum component in X-direction

$$min|S_{x,i,j,k}^0-S_{x,u*,w*}^d|$$

(22)

• Step 3. As shown in equation (23), let the point with smallest coordinate in Y direction be the point $S_{i*,j*,k*}$ and then $i*$ and $k*$ can be determined. Then, assign the X-direction coordinate of point $S_{i*,j*,k*}$ to the point $S_{u*,w*}^m$ as equation (24)

$$S_{y,u*,w*}^m=S_{y,i*,j*,k*}=\min S_{y,i,j,k_i}(i=1,2,\mathrm{\dots },N_t)$$

(23)

$$S_{x,u*,w*}^m=S_{x,i*,j*,k*}$$

(24)

• Step 4. By repeating above steps (1–3), find out corresponding points on the machined surface for all given points on the desired surface, and finally, the actual machined surface can be reconstructed.

Form error evaluation

After obtaining the machined surface $S_{u,w}^m$ , this article applies the following three indexes to evaluate the form error

$$\begin{array}{l}\mathrm{Average\; form\; error}:\\ \mathrm{AFE}=\frac{1}{UW}{\displaystyle \sum _{u=1}^U{\displaystyle \sum _{w=1}^W(S_{u,w}^m-S_{u,w}^d)}}\end{array}$$

(25)

$$\begin{array}{l}\mathrm{Mean\; absolute\; value\; of\; form\; error}:\\ \mathrm{MAE}=\frac{1}{UW}{\displaystyle \sum _{u=1}^U{\displaystyle \sum _{w=1}^W\left|S_{u,w}^m-\mathrm{AFE}\right|}}\end{array}$$

(26)

$$\begin{array}{l}\mathrm{Roughness\; along\; feed\; direction\; on}w*\mathrm{th\; contour}:\\ \mathrm{Ra}=\frac{1}{\mathrm{U}}\sum _{\mathrm{u}=1}^{\mathrm{U}}|{\mathrm{S}}_{\mathrm{u},\mathrm{w}*}^{\mathrm{m}}-{\overline{\mathrm{S}}}_{\mathrm{u},\mathrm{w}*}^{\mathrm{m}}|\end{array}$$

(27)

Simulation and experiment verification

Experiment setup

In this section, a series of milling tests are conducted to verify the effectiveness of the proposed method. As shown in Figure 4, a machining center and a three teeth carbide tool with a diameter of 16 mm, a helix angle of 45° and an overhang length of 80 mm are used in the milling tests. The workpiece material is aluminum alloy 7050. For the milling condition with a large axial depth of cut, in order to accurately reflect the dynamic characteristics of the cutter system and calculate the displacement response of the cutter, the cutter system is assumed as a multi-DOF flexible system. The experiment modal analysis with least mean square (LMS) is used to identify the modal parameters of cutter system. As shown in Figure 4, there exist five nodes in X and Y directions, and the DOF nodes are located at the axial position of 0, 10, 20, 30 and 40 mm from the bottom of cutter. Table 1 contains the identified results of modal parameters including natural frequency, damping ratio and mass-normalized shape. For comparison, this article also assumes the cutter system as a single-DOF dynamic system. At this moment, the modal parameters are identified using the dominant modal of the origin frequency response function (FRF) of the first node. During the milling process, 12 milling tests with its parameters are given in Table 2, in which the cutter runout parameters and the specific cutting force coefficients calibrated based on the thin-plate milling experiment ³³ are also determined. After the milling, the machined surface contour is measured by a contacting aspheric surface contour measurer, the measurement resolution and speed of which are 0.125 μm and 1 mm/s, respectively.

Figure 4.

Milling experiment setup.

Table 1.

Modal parameters of cutter system.

Modal parameters of multi-DOF dynamic model
Direction	Mode no.	Natural frequency (Hz)	Damping ratio	Mass-normalized mode shape $(1/\sqrt{\mathrm{kg}})$
X	1	622.1	0.0682	0.7951, 0.7334, 0.7566, 0.7287, 0.6461
	2	820.3	0.0412	0.6980, 0.6839, 0.5498, 0.6039, 0.4022
	3	1823.6	0.0286	3.5068, 1.7726, 1.8035, 1.7641, 1.0497
	4	2106.3	0.0408	2.9462, 0.6300, 0.4893, 0.3854, 0.1665
	5	2293.7	0.0466	2.9684, 1.1816, 0.7994, 0.7365, 0.3214
Y	1	243.9	0.3568	0.3454, 0.3405, 0.3313, 0.2978, 0.3170
	2	665.3	0.0472	0.7374, 0.6182, 0.5351, 0.5615, 0.4512
	3	1067.2	0.0816	1.4453, 1.2701, 1.0969, 0.9493, 0.8111
	4	1831.0	0.0218	3.0189, 2.7316, 3.1013, 1.5870, 1.1196
	5	2118.2	0.0548	2.6620, 1.8229, 1.0399, 0.6850, 0.5061
Modal parameters of single-DOF dynamic model
Direction		Mass (kg)	Damping (N/(m/s))	Stiffness (N/m)
X		0.081	50.9	1.065 × 107
Y		0.107	53.2	1.416 × 107

DOF: degree of freedom.

Table 2.

System parameters during milling process.

Cutter runout parameters		ρ (μm)	λ (°)		γ (°)	η (°)
		37.41	228.80		4.01 × 10−2	63.59
Specific cutting force coefficients		kts (N/mm2)	krs (N/mm2)		ktp (N/mm)	krp (N/mm)
		823.3	391.7		11.6	18.0
Milling parameters	Test no.	Type	n (r/min)	ap (mm)	ae (mm)	ft (mm/tooth)
	1	Down	9000	40	1	0.1
	2	Down	6000	35	0.3	0.1
	3	Down	9000	10	1	0.15
	4	Down	9000	30	1	0.1
	5	Down	8000	8	0.3	0.1
	6	Down	8000	5	2	0.1
	7	Down	11,000	10	1	0.15
	8	Down	11,000	6	2	0.15
	9	Up	6000	20	1	0.1
	10	Up	8000	3	2	0.1
	11	Up	9000	10	1.5	0.1
	12	Up	11,000	5	2	0.1

Simulation and experiment results

Figure 5(a)–(d) is the simulated results of machined surface topography and form error with the milling parameters in test 1 of Table 2 using the four kinds of simulation methods, respectively. It can be seen that the form error has an evident fluctuation in the cutter axis direction, but the variation in the feed direction is small. The form error distribution of different methods is different. During down milling process, the positive value of form error means under cut, and negative value of form error means over cut. During up milling, the conclusion is just the opposite.

Figure 5.

Simulated and measured form error of Test 1: (a) method 1, (b) method 2, (c) method 3, (d) method 4, (e) simulated and measured data of form error along cutter axis direction and (f) simulated and measured data of scallop height along feed direction

For further analyzing, Figure 5(e) shows the simulated and measured data of form error distribution along the cutter axis direction. Here, we can see that the result of method 1 has a good agreement with the measured data, the results of method 2 and method 3 have a slight distinction compared with the measured data, and there exists a large difference between method 4 and the measured data. The reason for the conclusion is that method 1 considers the tool/workpiece dynamic engagement process, which is more accordant with the actual situation. Method 2 and method 3 ignore the force–displacement coupling effect during cutter/workpiece engagement process, especially method 3 applies the cutter static deflection as the cutter response. Hence, both of them will not reflect the actual cutter response accurately. Only considering the cutter/workpiece geometric engagement, the simulated result of tooth trajectory in method 4 is mainly influenced by the cutter runout. As the cutter runout increases the cutting radius, the machined surface is easily to be over cut. Compared with other three kinds of methods, only method 4 shows the form error as a negative value. Figure 5(f) shows the simulated and measured data of scallop height distribution along the feed direction, and the height position of the contour is located at 1 mm from cutter bottom. The result indicates that the simulated data agrees with the measured one both in variation trend and amplitude. Besides, data variation period is about 0.3 mm, which is equal to the feed per revolution but not feed per tooth. This is because some teeth lost cut during the exit stage resulting from the existence of cutter runout.

Similarly, Figure 6 is the simulated and measured data of form error with its milling parameters given in test 2 of Table 2. From Figure 6(a)–(d), it can be concluded that the form error shows a large variation along cutter axis direction under the effect of cutter runout and dynamic excitation of milling force. Compared with other three kinds of method, the result of method 1 is much closer to the measured data. Besides, as only considering the geometric engagement process between cutter and workpiece, the simulated result of method 4 in Figure 6(e) is the same with that of Figure 5(e) at the same height position, but there is a large difference for the other three methods when the parameters are different, which means that the milling parameters have an important influence on form error. What’s more, the simulated and measured results of scallop height are also very close.

Figure 6.

Simulated and measured form error of Test 2: (a) method 1, (b) method 2, (c) method 3, (d) method 4, (e) simulated and measured data of form error along cutter axis direction and (f) simulated and measured data of scallop height along feed direction.

Finally, Table 3 gives out the evaluation results of form error of the above two milling tests, which includes the average form error (AFE), the mean absolute value of form error (MAE) and the contour roughness Ra. Similar to the above analysis results, the quantitative results show that the method 1 has higher prediction accuracy in predicting the form error when compared with other three kinds of methods. Table 4 further compares the measured data of form error with the predicted one based on method 1 from test 3 to test 12. The comparison indicates that the results of method 1 have a good agreement with the measured values.

Table 3.

Evaluation results of form error.

Test no.	Method	AFE (μm)	MAE (μm)	Ra (μm)
1	Method 1	15.19	4.66	0.36
	Method 2	13.43	2.02	0.36
	Method 3	13.34	3.60	0.35
	Method 4	−11.14	3.24	0.36
	Measured	16.75	6.15	0.35
2	Method 1	−1.26	2.61	0.36
	Method 2	1.96	2.61	0.37
	Method 3	1.54	3.22	0.35
	Method 4	−10.47	3.53	0.36
	Measured	−3.25	3.86	0.36

AFE: average form error; MAE: mean absolute value of form error.

Table 4.

Form error comparison between measured and predicted data with method 1.

Test no.	AFE (μm)		MAE (μm)		Ra (μm)
	Method 1	Measured	Method 1	Measured	Method 1	Measured
3	4.1	7.5	5.7	3.4	0.36	0.35
4	11.3	15.3	2.4	5.8	0.35	0.37
5	−7.9	−3.5	2.5	0.2	0.38	0.41
6	−4.4	−4.9	3.4	8.1	0.42	0.39
7	8.1	5.6	7.5	3.6	0.36	0.36
8	9.8	10.9	14.3	17.6	0.31	0.29
9	−15.4	−10.3	4.8	9.4	0.36	0.35
10	14.5	19.8	3.3	4.7	0.26	0.29
11	−8.4	−14.5	12.9	7.7	0.36	0.36
12	1.6	3.4	7.3	3.2	0.33	0.37

AFE: average form error; MAE: mean absolute value of form error.

The influence of different milling conditions on form error

This section focuses on the influence of different milling conditions on the form error through simulation and experiment studies. In order to analyze the influence trend more conveniently, the first four subsections only apply the single-DOF dynamics model of cutter system with its modal parameters given in Table 1 to study the effect of milling parameters, cutter runout parameters, milling type and dynamics parameters of cutter system, respectively, while the last subsection analyzes the distinction between single-DOF and multi-DOF dynamics model of cutter system.

The influence of spindle speed n and radial immersion ratio aD on form error

Figure 7(a) and (b) is the simulated results of AFE and Ra under different spindle speeds n and radial immersion ratios aD using method 1 without cutter runout effect, respectively. The machining parameters during the down milling are axial depth of cut a_p = 1 mm and feed per tooth f_t = 0.1 mm/tooth. It can be seen that the AFE and Ra show an evident fluctuation within the whole parameter field, which means that the AFE and Ra have a strong dependence on the milling parameters. Besides, both AFE and Ra have multiple local extrema, and the position of these extrema makes the combination of spindle speed and radial immersion ratio change linearly.

Figure 7.

(a) AFE under different n and aD, (b) Ra under different n and aD, (c) AFE under different n with aD = 0.3 and (d) Ra under different n with aD = 0.3.

Furthermore, Figure 7(c) and (d) give out the detailed results of Figure 7(a) and (b) with radial immersion ratio aD = 0.3, respectively. We can see that the AFE and Ra vary evidently with spindle speed. Besides, what should be pointed out is that the extremum is located around some certain spindle speeds, which can easily cause cutter system resonance as a result of harmonic components of milling force. In the single-DOF cutter dynamic system, the natural frequency in the normal feed direction is about 1824 Hz. When the spindle speed is around 12.16, 9.12, 7.29 and 6.08 krpm, the third-, fourth-, fifth- and sixth-order harmonic of tooth passing frequency of milling force is coincident with the natural frequency of the cutter system, respectively. Here, the resonance spindle speed and natural frequency will satisfy $n=60f_n/(k{N}_t)(k=1,2,3,\dots )$ . At this moment, both AFE and Ra will change rapidly and produce extremum, and the variation of AFE will shows alternate change between under cut and over cut around the resonance spindle speeds. Additionally, combining Figure 7(a) and (b), the number of extrema will increase when the spindle speed decreases. The reason for this phenomenon is that the harmonics of milling force causing resonance under low spindle speed are the high-order components, which have a relatively low period, lead to more number of variation cycles within the same range of radial immersion ratio, and finally, increase the number of extrema.

Beyond that, the results of the other three methods are also plotted in the Figure 7(c) and (d). By contrast, the simulated results of method 1 and method 2 have a slight distinction, especially around the resonance and low spindle speed, which indicates that the force–displacement coupling effect has a non-ignored influence on the form error. As the dynamic excitation is neglected, both results of AFE and Ra using method 3 and method 4 are a constant value, in which the AFE of method 4 is equal to zero. The comparison demonstrates that both method 3 and method 4 cannot obtain an accurate simulation results for form error.

The influence of cutter runout on form error

Figure 8(a) and (b) are the simulated results of AFE and Ra under different spindle speeds n and radial immersion ratios aD using method 1 with consideration of cutter runout effect, respectively. The milling parameters are the same with that in section “The influence of spindle speed n and radial immersion ratio aD on form error”, and the cutter runout parameters are $\rho =20\mu \mathrm{m},\lambda =0\circ ,\gamma =0\circ ,\eta =0\circ$ . When comparing with Figure 7(a) and (b), the distribution of the AFE and Ra within the whole milling parameter field here are more complex and the number of extrema has an evident growth.

Figure 8.

(a) AFE under different n and aD, (b) Ra under different n and aD, (c) AFE under different n with aD = 0.3, (d) Ra under different n with aD = 0.3, (e) AFE and MAE under different offset value and (f) AFE and MAE under different inclination value.

Similarly, Figure 8(c) and (d) further detail the results of Figure 8(a) and (b) with radial immersion ratio aD = 0.3. It seems that the number of local extrema of AFE and Ra increases when comparing with Figure 7(c) and (d). This is because that the cutter runout would cause new harmonic components of milling force, which is equal to the spindle passing frequency and its harmonics, and can also excite the cutter system to produce resonance. Hence, when the resonance spindle speed and natural frequency satisfy $n=60f_n/k(k=1,2,3,\dots )$ , both AFE and Ra will vary rapidly. Besides, in Figure 8(c), as the actual cutting radius increase about 20 μm as a result of cutter runout, which leads to the AFE of method 1 and method 2 varying around −20 μm as a kind of over cut. As shown in Figure 8(d), the Ra has a significant increase compared with Figure 7(d), which is due to the lost cut of some teeth during the exit stage and the increase in the contour scallop height as a result of cutter runout.

At last, Figure 8(e) and (f) analyze the effect of cutter runout offset value and inclination value on the AFE and MAE, respectively. Here, the axial depth of cut a_p is 10 mm, and the radial immersion ratio aD is 0.03. We can see that under the milling condition with a very-low radial immersion ratio, the influence of spindle speed on form error will become weak. Besides, the MAE shows a linear increase with the increase in the cutter runout offset and inclination value, while the AFE shows a linear decrease with the same range of cutter runout parameters, which means the over cut condition will become more and more serious with the increase in the offset and inclination value. In addition, the inclination value is more likely to lead to the increase in the form error than offset value.

The influence of milling type on form error

Figure 9(a) and (b) are the simulated results of AFE and Ra under different milling types using method 1 without cutter runout effect, respectively. The machining parameters during the down milling are axial depth of cut a_p = 1 mm and feed per tooth f_t = 0.1 mm/tooth. The results show that the AFE and Ra with the same radial immersion ratio aD under different milling types are different. From Figure 9(a), when the radial immersion ratio is equal to 0.3, the amplitude fluctuation ranges of up and down milling are close to each other. But under the radial immersion ratio with 0.05, the amplitude fluctuation range of up milling is far less than that of down milling. The simulated results of Ra shown in Figure 9(b) have the same conclusion. In fact, the manufacturer always applies down milling during the finishing operation, but the study shows that some potential milling conditions in up milling can obtain a higher form precision than down milling, which can provide more options for machining process.

Figure 9.

(a) AFE under different milling types and aD and (b) Ra under different milling types and aD.

The influence of damping and stiffness of cutter system on form error

Figure 10 is the simulated results of AFE and Ra under different damping and stiffness of cutter system using method 1 and method 2 without cutter runout effect, respectively. The natural frequency of the system remains unchanged when changing the damping and stiffness. The variable C_p in Figure 10(a) and (b) is the proportion between the simulated damping value and experiment value, and the variable K_p in Figure 10(c) and (d) is the proportion between the simulated stiffness value and experiment value. The experiment value of damping and stiffness are given in Table 1. The machining parameters during the down milling are axial depth of cut a_p = 1 mm, radial immersion ratio aD = 0.3 and feed per tooth f_t = 0.1 mm/tooth. It can be evidently seen that with the increase in the damping and stiffness, the amplitude fluctuation range of AFE and Ra within the whole spindle speed field shows a decrease trend, especially around the resonance spindle speed, which means that a large damping and stiffness can restrain the cutter vibration amplitude and reduce the dependence of AFE and Ra on spindle speed.

Figure 10.

(a) AFE under different damping, (b) Ra under different damping, (c) AFE under different stiffness and (d) Ra under different stiffness.

Besides, from the comparison between method 1 and method 2, we can see that when the damping and stiffness are large, the calculation results of the two methods are close to each other, but when the damping and stiffness are low, the deviation of the two methods become remarkable, especially around the resonance spindle speed. This result means that the force–displacement coupling effect during cutter/workpiece engagement process is strong when the damping and stiffness of the system is low. Hence, for a milling system with low damping and stiffness, we should apply the method 1 to analyze the form error.

The influence of multi-DOF cutter system on form error

Figure 11(a) shows the simulated results of AFE under single-DOF and multi-DOF dynamics model of cutter system using method 1 with cutter runout effect. The machining parameters during the down milling are axial depth of cut a_p = 1 mm, radial immersion ratio aD = 0.03 and feed per tooth f_t = 0.1 mm/tooth, and the cutter runout is the same with that in Table 2. As multi-DOF model of cutter dynamics model considers the influence of multi-order natural frequencies, the number of resonance spindle speeds will increase compared with the single-DOF model. Besides, the simulated results of AFE by the multi-DOF model under different spindle speeds are closer to the measured ones than the single-DOF model, which means that it is necessary to consider the multi-order modal characteristics of cutter system.

Figure 11.

(a) Simulated and measured results of AFE, (b) AFE under different a_p with n = 9000 and (c) AFE under different a_p with n = 6000.

For further comparing, Figure 11(b) and (c) shows the simulated results of AFE under different axial depth of cut using method 1 with multi-DOF model and single-DOF model within stable milling, respectively, where the measured data are also plotted in the figures. Here, other machining parameters simulated in Figure 11(b) and (c) are the same with that in test 1 and test 2 of Table 2, respectively. It can be seen that the simulated results of the two kinds of models are close when the axial depth of cut is small, but gradually deviate with each other with the increase in the axial depth of cut. Besides, the simulated results show a good agreement with the measured ones. The above study means that the multi-DOF model can describe the cutter/workpiece engagement process more accurate than the single-DOF model, especially under the milling condition with a large axial depth of cut.

Conclusion

This article presents an enhanced time-domain simulation method for cutter/workpiece dynamic engagement process with consideration of the multi-DOF dynamics model of cutter system and cutter runout and further builds up a reconstruction algorithm of machined surface and evaluates the form error. Besides, the study systematically compares the performances between the proposed method and other three kinds of methods, and the study results mean that the proposed method has the best performance in terms of simulation accuracy. The main conclusions of this article are given as follows:

1. The form error has an obvious dependence characteristic on milling parameters, especially on the spindle speed and radial depth of cut. Around the resonance spindle speed, the form error will change dramatically and lead to multi-extremum, which has been verified by the proposed method and method 2, but method 3 and method 4 cannot describe this kind of variation.

2. The proposed multi-DOF dynamics model of cutter system can more accurately reflect the actual cutter/workpiece engagement process than the single-DOF model, especially under the milling condition with a large axial depth of cut, where the multi-modal characteristics of cutter system will be prominent.

3. For a cutter system with a low damping or stiffness, the force–displacement coupling effect during tool/workpiece engagement is strong, and the proposed method can obtain higher calculation accuracy than other methods. While for a cutter system with large damping or stiffness, the method 2 can also achieve an accurate result, which is close to the proposed method.

4. The existence of cutter runout would not only increase the actual cutting radius, and lead to over cut, but also cause new harmonic components of milling force, then result in the growth of number of resonance spindle speeds. Besides, it would change the tooth cutting state during the start/exit stage, thus seriously deteriorate the surface roughness.