METHOD FOR SIMULATING CHATTER-FREE MILLED SURFACE TOPOGRAPHY

Abstract

Method for simulating chatter-free milled surface topography includes: performing dynamics modeling on milling system and building delay differential equation with multiple delays to represent dynamic model; constructing state transition matrix between two adjacent cutting tool rotation periods by extending GRK method; predicting stable milling parameter zones based on Floquet theory; calculating vibration displacements on discretized time nodes during one cutting tool rotation period by means of fixed mapping point theorem; constructing motion trajectories of cutting tool edges in normal and feed directions with regard to milled surface of workpiece; predicting surface topography of workpiece by using spline interpolation to densify motion trajectories of cutting tool edges that participate in generating final surface of workpiece; calculating surface location error of milled surface based on motion trajectories of cutting tool edges in normal direction with regard to milled surface of workpiece; calculating surface roughness based on the simulated surface topography of workpiece.

Claims

1. A method for simulating chatter-free milled surface topography, comprising steps of: step 1: performing dynamics modeling on milling system and building a delay differential equation with multiple delays to represent dynamic model; step 2: constructing state transition matrix between two adjacent cutting tool rotation periods by extending Generalized Runge-Kutta method; step 3: predicting stable milling parameter zones based on Floquet theory; step 4: calculating vibration displacements on discretized time nodes during one cutting tool rotation period by means of fixed mapping point theorem; step 5: constructing motion trajectories of cutting tool edges in normal and feed directions with regard to milled surface of workpiece; step 6: predicting surface topography of workpiece by using spline interpolation to densify motion trajectories of cutting tool edges that participate in generating final surface of workpiece; step 7: calculating surface location error of milled surface based on motion trajectories of cutting tool edges in normal direction with regard to milled surface of workpiece; step 8: calculating surface roughness based on the simulated surface topography of workpiece.

2. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 1 comprises sub-steps of: step 1.1: considering flexibility of both cutting tool and workpiece, variance of pitch angles and helix angles of cutting tool, and teeth runout, performing dynamics modeling on milling system and building a delay differential equation with multiple delays to represent the dynamic model in physical space: $M\stackrel{.Math.}{q}\left(t\right)+C\stackrel{.}{q}\left(t\right)+\mathrm{Kq}\left(t\right)=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{K_c\left(t,j,k\right)\left[q\left(t\right)-q\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+F_0\left(t,j,k\right)\right\}$ where M C and K are mass, damping and stiffness matrices, respectively; {umlaut over (q)}(t), {dot over (q)}(t) and q(t) are acceleration, velocity and displacement vectors at time t, respectively; K.sub.c(t,j,k) is directional coefficient matrix for Tooth k at j-th axial slice at time t; F.sub.0(t,j,k) is cutting force vector for Tooth k at j-th axial slice at time t, which is independent of chip regeneration; Σ.sub.p=k.sup.kτ(j,p) is time delay with regard to cutting element (j,k) is mathematical summation operator, p is a variable for summation, and is tooth serial number to start delay calculation; N is number of teeth; and N.sub.a is number of axial discretization of tool; step 1.2: performing modal coordinate transformation on the dynamic model in step 1.1, which can be transformed from physical space into modal space; the modal coordinate transformation equation is:
q(t)=PΓ(t) where P is mode shape matrix, and Γ(t) is modal displacement vector at time t; dynamics model in modal space which is represented with delay differential equation with multiple delays turns out to be: $M_{\Gamma }\stackrel{.Math.}{\Gamma }\left(t\right)+C_{\Gamma }\stackrel{.}{\Gamma }\left(t\right)+K_{\Gamma }\Gamma \left(t\right)=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{P^T{K}_c\left(t,j,k\right)P\left[\Gamma \left(t\right)-\Gamma \left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+P^T{F}_0\left(t,j,k\right)\right\}$ where M.sub.Γ, C.sub.Γ and K.sub.Γ are modal mass, modal damping and modal stiffness matrices, respectively; {umlaut over (Γ)}(t), {dot over (Γ)}(t) and Γ(t) are modal acceleration, modal velocity and modal displacement vectors, respectively; for milling system with flexible cutting tool and flexible workpiece, the dynamic model becomes: $\begin{array}{c}\left[\begin{array}{cc}{M}_{\Gamma ;T}& \\ & M_{\Gamma ;W}\end{array}\right]\left\{\begin{array}{c}{\stackrel{.Math.}{\Gamma }}_T\left(t\right)\\ {\stackrel{.Math.}{\Gamma }}_W\left(t\right)\end{array}\right\}+\left[\begin{array}{cc}{C}_{\Gamma ;T}& \\ & C_{\Gamma ;W}\end{array}\right]\left\{\begin{array}{c}{\stackrel{.}{\Gamma }}_T\left(t\right)\\ {\stackrel{.}{\Gamma }}_W\left(t\right)\end{array}\right\}+\left[\begin{array}{cc}{K}_{\Gamma ;T}& \\ & K_{\Gamma ;W}\end{array}\right]\left\{\begin{array}{c}{\Gamma }_T\left(t\right)\\ {\Gamma }_W\left(t\right)\end{array}\right\}=\\ \underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{\left[\begin{array}{c}{P}_T^T\\ -P_W^T\end{array}\right]K_c\left(t,j,k\right)\begin{array}{cc}[P_T& -P_W]\left\{\begin{array}{c}{\Gamma }_T\left(t\right)-{\Gamma }_T\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\\ {\Gamma }_W\left(t\right)-{\Gamma }_W\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\end{array}\right\}\end{array}+\left[\begin{array}{c}{P}_T^T\\ -P_W^T\end{array}\right]F_0\left(t,j,k\right)\right\}\end{array}$ where M.sub.Γ;T, C.sub.Γ;T and K.sub.Γ;T are modal mass, modal damping and modal stiffness matrices of cutting tool; M.sub.Γ;W, C.sub.Γ;W and K.sub.Γ;W are modal mass; modal damping and modal stiffness matrices of workpiece; {umlaut over (Γ)}.sub.T(t), {dot over (Γ)}.sub.T (t) and Γ.sub.T(t) are modal acceleration, modal velocity and modal displacement vectors of cutting tool; {umlaut over (Γ)}.sub.W(t), Γ.sub.W(t) and Γ.sub.W(t) are modal acceleration, modal velocity and modal displacement vectors of workpiece; P.sub.T, and P.sub.W are mode shape matrices of cutting tool and workpiece, respectively; P.sub.T.sup.T and P.sub.W.sup.T are transposed matrices P.sub.T and P.sub.W, respectively; for milling system with flexible cutting tool and rigid workpiece, the dynamic model becomes: $\begin{array}{c}{M}_{\Gamma ;T}{\stackrel{.Math.}{\Gamma }}_T\left(t\right)+C_{\Gamma ;T}{\stackrel{.}{\Gamma }}_T\left(t\right)+K_{\Gamma ;T}{\Gamma }_T\left(t\right)=\\ \underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{P_T^T{K}_c\left(t,j,k\right)\left\{P_T\left[{\Gamma }_T\left(t\right)-{\Gamma }_T\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]\right\}+P_T^T{F}_0\left(t,j,k\right)\right\}\end{array}$ for milling system with rigid cutting tool and flexible workpiece, the dynamic model becomes: $\begin{array}{c}{M}_{\Gamma ;W}{\stackrel{.Math.}{\Gamma }}_W\left(t\right)+C_{\Gamma ;W}{\stackrel{.}{\Gamma }}_W\left(t\right)+K_{\Gamma ;W}{\Gamma }_W\left(t\right)=\\ \underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{P_W^T{K}_c\left(t,j,k\right)\left\{P_W\left[{\Gamma }_W\left(t\right)-{\Gamma }_W\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]\right\}-P_W^T{F}_0\left(t,j,k\right)\right\}.\end{array}$

3. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 2 comprises sub-steps of: step 2.1: performing state space transformation on the dynamic model in modal space by introducing state variable vector x(t)=[(Γ(t)).sup.T,({dot over (Γ)}(t)).sup.T].sup.T, and delay differential equation with multiple delays in state space becomes: $\stackrel{.Math.}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{B\left(t,j,k\right)\left[x\left(t\right)-x\left(t-{.Math.}_{p={\stackrel{\_}{k}}_v}^k\tau \left(j,p\right)\right)\right]+D\left(t,j,k\right)\right\}$ where x(t) denotes state variable vector, {dot over (x)}(t) denotes first derivative of x(t) with regard to time t, $A=\left[\begin{array}{cc}0& I\\ -M_{\Gamma }^{-1}{K}_{\Gamma }& -M_{\Gamma }^{-1}{C}_{\Gamma }\end{array}\right],$ $B\left(t,j,k\right)=\left[\begin{array}{cc}0& 0\\ M_{\Gamma }^{-1}{P}^T{K}_c\left(t,j,k\right)P& 0\end{array}\right],$ $D\left(t,j,k\right)=\left[\begin{array}{c}0\\ M_{\Gamma }^{-1}{P}^T{F}_0\left(t,j,k\right)\end{array}\right],$ and I is unit matrix; step 2.2: dividing one rotation period of cutting tool T into 2i.sub.m+2 equal intervals with a series of discrete nodes {t.sub.0,t.sub.1, . . . t.sub.2i.sub.m.sub.+2}; analytical expression of dynamic responses at interval [t.sub.i,t] with regard to the delay differential equation with multiple delays in step 2.1 is: $x\left(t\right)=e^{A\left(\tau -t_i\right)}x\left(t_i\right)+\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{\int }_{t_i}^t\left\{e^{A\left(t-\xi \right)}B\left(\xi ,j,k\right)\left[x\left(\xi \right)-x\left(\xi -{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+e^{A\left(t-\xi \right)}D\left(\xi ,j,k\right)\right\}d\xi$ where ξ is integration variable; for brevity, define x.sub.i:=x(t.sub.i), E(t,t.sub.i):=e.sup.A(t−t.sup.i.sup.)x(t.sub.i), H.sub.j,k(t,ξ,x(ξ)):=e.sup.A(t−ξ)B(ξ,j,k)[x(ξ)−x(ξ−τ(j,k,p))] and J.sub.j,k(t,ξ):=e.sup.A(t+ξ)D(ξ,j,k); on sub-interval [t.sub.2i,t.sub.2i+2], x(t) is approximated using Simpson's rule: $x_{2i+2}=E\left(t_{2i+2},t_{2i}\right)+\frac{h}{3}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{H_{j,k}\left(t_{2i+2},t_{2i},x_{2i}\right)+4H_{j,k}\left(t_{2i+2},t_{2i+1},x_{2i+1}\right)+H_{j,k}\left(t_{2i+2},t_{2i+2},x_{2i+2}\right)+J_{j,k}\left(t_{2i+2},t_{2i}\right)+4J_{j,k}\left(t_{2i+2},t_{2i+1}\right)+J_{j,k}\left(t_{2i+2},t_{2i+2}\right)\right\}$ on sub-interval [t.sub.2i,t.sub.2i+1], x(i) is approximated using fourth-order Runge-Kutta algorithm: $x_{2i+1}=E\left(t_{2i+1},t_{2i}\right)+\frac{h}{6}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{H_{j,k}\left(t_{2i+1},t_{2i},x_{2i}\right)+2H_{j,k}\left(t_{2i+1},t_{2i+1/2},x_{2i+1/2}\right)+2H_{j,k}\left(t_{2i+1},t_{2i+1/2},x_{2i+1/2}\right)+H_{j,k}\left(t_{2i+1},t_{2i+1},x_{2i+1}\right)+J_{j,k}\left(t_{2i+1},t_{2i}\right)+2J_{j,k}\left(t_{2i+1},t_{2i+1/2}\right)+2J_{j,k}\left(t_{2i+1},t_{2i+1/2}\right)+J_{j,k}\left(t_{2i+1},t_{2i+1}\right)\right\}$ where state variable vector at middle points t.sub.2i+1/2 is approximated with Lagrange formula:
x.sub.2i+1/2=⅜x.sub.2i+¾x.sub.2i−1−⅛x.sub.2i+2 re-organizing the analytical expression of dynamic responses on sub-interval [t.sub.2i,t.sub.2i+2] into: $F_{2i}{x}_{2i}+F_{2i+1}{x}_{2i+1}+F_{2i+2}{x}_{2i+2}=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left(F_{2i-\tau \left(j,k\right)}{x}_{2i-\tau \left(j,k\right)}+F_{2i+1-\tau \left(j,k\right)}{x}_{2i+1-\tau \left(j,k\right)}+F_{2i+2-\tau \left(j,k\right)}{x}_{2i+2-\tau \left(j,k\right)}\right)+\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left(S_{2i}+S_{2i+1/2}+S_{2i+1}\right)$ $\mathrm{where}$ $F_{2i}=-e^{\mathrm{Ah}}-\frac{h}{6}{e}^{\mathrm{Ah}}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i,j,k}-\frac{3}{8}.Math.\frac{4h}{6}{e}^{\mathrm{Ah}/2}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i+1/2,j,k},$ $F_{2i+1}=I-\frac{3}{4}.Math.\frac{4h}{6}{e}^{\mathrm{Ah}/2}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i+1/2,j,k}-\frac{h}{6}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i+1,j,k},$ $F_{2i+2}=\frac{1}{8}.Math.\frac{4h}{6}{e}^{\mathrm{Ah}/2}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i+1/2,j,k},$ $F_{2i-\tau \left(j,k,p\right)}=-\frac{h}{6}{e}^{\mathrm{Ah}}{B}_{2i,j,k}-\frac{3}{8}.Math.\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{2i+1/2,j,k},$ $F_{2i+1-\tau \left(j,k,p\right)}=-\frac{3}{4}.Math.\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{2i+1/2,j,k}-\frac{h}{6}{B}_{2i+1,j,k},$ $F_{2i+2-\tau \left(j,k,p\right)}=\frac{1}{8}.Math.\frac{4h}{6}{e}^{\mathrm{Ah}/2}{B}_{2i+1/2,j,k},$ $S_{2i}=\frac{h}{6}{e}^{A.Math.h}{D}_{2i,j,k},$ $S_{2i+1/2}=\frac{h}{6}.Math.4e^{A.Math.h/2}{D}_{2i+1/2,j,k},$ $S_{2i+1}=\frac{h}{6}{D}_{2i+1,j,k},$ and h is the discretization step; re-organizing the analytical expression of dynamic responses on sub-interval [t.sub.2i,t.sub.2i+2] into: $G_{2i}{x}_{2i}+G_{2i+1}{x}_{2i+1}+G_{2i+2}{x}_{2i+2}=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left(G_{2i-\tau \left(j,k\right)}{x}_{2i-\tau \left(j,k\right)}+G_{2i+1-\tau \left(j,k\right)}{x}_{2i+1-\tau \left(j,k\right)}+G_{2i+2-\tau \left(j,k\right)}{x}_{2i+2-\tau \left(j,k\right)}\right)+\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left(T_{2i}+T_{2i+1}+T_{2i+2}\right)$ $\mathrm{where}$ $G_{2i}=-e^{2\mathrm{Ah}}-\frac{h}{3}{e}^{2\mathrm{Ah}}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i,j,k},$ $G_{2i+1}=-\frac{4h}{3}{e}^{Ah}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i+1,j,k},$ $G_{2i+2}=I-\frac{h}{3}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{B}_{2i+2,j,k},$ $G_{2i-\tau \left(j,k,p\right)}=-\frac{h}{3}{e}^{2Ah}{B}_{2i,j,k},$ $G_{2i+1-\tau \left(j,k,p\right)}=-\frac{4h}{3}{e}^{Ah}{B}_{2i+1,j,k},$ $G_{2i+2-\tau \left(j,k,p\right)}=-\frac{h}{3}{B}_{2i+2,j,k},$ $T_{2i}=\frac{h}{3}{e}^{A.Math.2h}{D}_{2i,j,k},$ $T_{2i+1}=\frac{h}{3}.Math.4e^{A.Math.h}{D}_{2i+1,j,k}$ $\mathrm{and}$ $T_{2i+2}=\frac{h}{3}{D}_{2i+2,j,k};$ step 2.3: based on the formulae in step 2.2, establishing discrete mapping relationship between two adjacent rotation periods of cutting tool [−T,0] and [0,T]: $P_1{y}_{\left[0,T\right]}={\mathrm{Qy}}_{\left[-T,0\right]}+P_2{y}_{\left[0,T\right]}+z_{\left[0,T\right]}$ $\mathrm{where}$ $y_{\left[0,T\right]}=\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ x_4\\ .Math.\\ x_{2i_m+1}\\ x_{2i_m+2}\end{array}\right],y_{\left[-T,0\right]}=\left[\begin{array}{c}{x}_{0-T}\\ x_{1-T}\\ x_{2-T}\\ x_{3-T}\\ x_{4-T}\\ .Math.\\ x_{2i_m+1-T}\\ x_{2i_m+2-T}\end{array}\right],z_{\left[0,T\right]}=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left[\begin{array}{c}0\\ S_0+S_{1/2}+S_1\\ T_0+T_1+T_2\\ S_2+S_{2+1/2}+S_3\\ T_2+T_3+T_4\\ .Math.\\ S_{2i_m}+S_{2i_m+1/2}+S_{2i_m+1}\\ T_{2i_m}+T_{2i_m+1}+T_{2i_m+2}\end{array}\right],$ $P_1=\left[\begin{array}{ccccccccc}I& 0& 0& 0& 0& .Math.& 0& 0& 0\\ F_0& F_1& F_2& 0& 0& .Math.& 0& 0& 0\\ G_0& G_1& G_2& 0& 0& .Math.& 0& 0& 0\\ 0& 0& F_2& F_3& F_4& .Math.& 0& 0& 0\\ 0& 0& G_2& G_3& G_4& .Math.& 0& 0& 0\\ .Math.& .Math.& .Math.& .Math.& .Math.& \ddots & .Math.& .Math.& .Math.\\ 0& 0& 0& 0& 0& .Math.& F_{2i_m}& F_{2i_m+1}& F_{2i_m+2}\\ 0& 0& 0& 0& 0& .Math.& G_{2i_m}& G_{2i_m+1}& G_{2i_m+2}\end{array}\right],$ $Q=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left[\begin{array}{ccccccccc}0& .Math.& 0& 0& 0& .Math.& 0& 0& I/{\mathrm{NN}}_a\\ 0& .Math.& F_{-\tau \left(j,k,p\right)}& F_{1-\tau \left(j,k,p\right)}& F_{2-\tau \left(j,k,p\right)}& .Math.& 0& 0& 0\\ 0& .Math.& G_{-\tau \left(j,k,p\right)}& G_{1-\tau \left(j,k,p\right)}& G_{2-\tau \left(j,k,p\right)}& \ddots & .Math.& 0& 0\\ .Math.& \ddots & .Math.& .Math.& .Math.& \ddots & .Math.& .Math.& .Math.\\ 0& .Math.& 0& 0& 0& .Math.& F_{2i-\tau \left(j,k,p\right)}& F_{2i+1-\tau \left(j,k,p\right)}& F_{2i+2-\tau \left(j,k,p\right)}\\ 0& .Math.& 0& 0& 0& .Math.& G_{2i-\tau \left(j,k,p\right)}& G_{2i+2-\tau \left(j,k,p\right)}& G_{2i+2-\tau \left(j,k,p\right)}\\ .Math.& \ddots & .Math.& .Math.& .Math.& \ddots & .Math.& .Math.& .Math.\\ 0& .Math.& 0& 0& 0& .Math.& 0& 0& 0\end{array}\right]$ $\mathrm{and}$ $P_2=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left[\begin{array}{ccccccccc}0& 0& 0& .Math.& 0& 0& 0& .Math.& 0\\ .Math.& .Math.& .Math.& \ddots & .Math.& .Math.& .Math.& .Math.& 0\\ F_{2i+2+\tau \left(j,k,p\right)}& F_{2i+3+\tau \left(j,k,p\right)}& F_{2i+4-\tau \left(j,k,p\right)}& .Math.& 0& 0& 0& .Math.& 0\\ G_{2i+2-\tau \left(j,k,p\right)}& G_{2i+3-\tau \left(j,k,p\right)}& G_{2i+4-\tau \left(j,k,p\right)}& \ddots & 0& 0& 0& .Math.& 0\\ 0& 0& .Math.& \ddots & 0& 0& 0& .Math.& 0\\ .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.& .Math.\\ 0& 0& 0& .Math.& F_{2i_m-\tau \left(j,k,p\right)}& F_{2i_m+1-\tau \left(j,k,p\right)}& F_{2i_m+2-\tau \left(j,k,p\right)}& .Math.& 0\\ 0& 0& 0& .Math.& G_{2i_m-\tau \left(j,k,p\right)}& G_{2i_m+1-\tau \left(j,k,p\right)}& G_{2i_m+2-\tau \left(j,k,p\right)}& .Math.& 0\end{array}\right].$

4. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 3 is: constructing transition matrix Φ of the milling system: $\Phi =\{\begin{array}{cc}{\left(P_1-P_2\right)}^{-1}Q,& \mathrm{if}.Math.P_1-P_2.Math.\ne 0\\ {\left(P_1-P_2\right)}^{†}Q,& \mathrm{if}.Math.P_1-P_2.Math.=0\end{array}$ where |P.sub.1−P.sub.2| denotes determinant of matrix (P.sub.1−P.sub.2); † denotes Moore-Penrose pseudoinverse of matrix (P.sub.1−P.sub.2); based on Floquet theory, stability of milling system depends on eigenvalue of Φ:it is stable if all modules of eigenvalues of Φ are less than 1; it is unstable otherwise.

5. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 4 is: based on fixed mapping point theorem, the state variable vector x(t) satisfies x(t.sub.i)=x(t.sub.i−T) for chatter-free cutting conditions; therefore, it has y.sub.[0,T]=y.sub.−T,0; calculating state variables on all discrete points by: $y^{*}=y_{\left[0,T\right]}=y_{\left[-T,0\right]}=\{\begin{array}{cc}{\left(P_1-P_2-Q\right)}^{-1}{z}_{\left[0,T\right]},& \mathrm{if}.Math.P_1-P_2-Q.Math.\ne 0\\ {\left(P_1-P_2-Q\right)}^{†}{z}_{\left[0,T\right]},& \mathrm{if}.Math.P_1-P_2-Q.Math.=0\end{array}$ where y* consists of modal displacement vector Γ(t.sub.i) and modal velocity vector {dot over (Γ)}(t.sub.i) at all discrete time nodes i=0,1, . . . ,2i.sub.m+2 during one rotation period of cutting tool.

6. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 5 comprises sub-steps of: step 5.1: transforming modal displacements into physical space, and calculating relative vibration displacements q(t.sub.i) between cutting tool and workpiece:
q(t.sub.i)=q.sub.T(t.sub.i)−q.sub.W(t.sub.i)=P.sub.TΓ.sub.T(t.sub.i)−P.sub.WΓ.sub.W(t.sub.i) where q.sub.T (t.sub.i) and q.sub.W(t.sub.i) are vibration displacements of cutting tool and workpiece, respectively; step 5.2: extracting vibration displacements in direction normal to milled surface of workpiece (defined as {right arrow over (Oy)} direction) from q(t.sub.i), which constructs a new vector Q*.sub.y:
Q.sub.y*={[q.sub.y(t.sub.1)]T,[q.sub.y(t.sub.2)].sup.T, . . . ,[q.sub.x(t.sub.2i.sub.m.sub.+2)]T} extracting vibration displacements in feed direction (defined as {right arrow over (Ox)} direction), which constructs a new vector Q.sub.y*:
Q.sub.x*={[q.sub.x(t.sub.1)].sup.T,[q.sub.x(t.sub.2)].sup.T, . . . ,[q.sub.x(t.sub.2i.sub.m.sub.+2)]T} step 53: based on kinematic synthesis of relative feed movement between cutting tool and workpiece, rotation of cutting tool and relative vibrations between cutting tool and workpiece, calculating motion trajectories in {right arrow over (Oy)} direction and {right arrow over (Ox)} direction for cutting element (j,k), respectively: ${\Theta }_y^{*}\left(i,j,k\right)=\underset{\mathrm{feed}}{\underbrace{\frac{f.Math.t_i}{1000\times 60}\sin \left({\phi }_{\mathrm{fy}}\right)}}-\underset{\mathrm{rotation}}{\underbrace{R\left(j,k\right)\cos \left({\phi }_{\mathrm{ref}}-\varphi \left(j,k\right)-\frac{z_j\tan \left({\beta }_k\right)}{R}+2\pi \frac{\Omega t_i}{60}\right)}}+\underset{\mathrm{vibration}}{\underbrace{Q_y^{*}\left(i,j\right)}}$ ${\Theta }_y^{*}\left(j,k,i\right)=\underset{\mathrm{feed}}{\underbrace{\frac{f.Math.t_i}{1000\times 60}\cos \left({\phi }_{\mathrm{fy}}\right)}}+\underset{\mathrm{rotation}}{\underbrace{R\left(j,k\right)\sin \left({\phi }_{\mathrm{ref}}+\varphi \left(j,k\right)-\frac{z_j\tan \left({\beta }_k\right)}{R}+\frac{60t_i}{\Omega }\right)}}+\underset{\mathrm{vibrations}}{\underbrace{Q_x^{*}\left(j,i\right)}}$ where f is relative feed rate between cutting tool and workpiece in mm/min; φ.sub.fy is angle between feed direction of cutting tool and normal direction of workpiece; R(j,k) is actual cutting radius of cutting element (j,k); φ.sub.ref is reference angle used in the Generalized Runge-Kutta method; ϕ(j,k) is pitch angle between Tooth k and Tooth 1 for j-th discrete axial slice; β.sub.k is helix angle with regard to Tooth k; z.sub.f is axial height of j-th discrete axial slice; R is geometric radius of cutting tool; Ω is spindle rotation speed; Q.sub.y*(j,i) and Q.sub.x*(j,i) are vibration displacements of cutting element (j,k) in {right arrow over (Oy)} and {right arrow over (Ox)} directions, respectively.

7. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 6 comprises sub-steps of: step 6.1: selecting part of motion trajectories of cutting tool that are near milled surface of workpiece based on following equation, which construct motion trajectory points to be densified with interpolation, i.e., (Θ.sub.x_trim*(j,k,i), Θ.sub.y_trim*(j,k,i)):
Θ.sub.y*(j,k,i)≥max{Θ.sub.y*(j,k,i)}−αR, down-milling
Θ.sub.y*(j,k,i)≤min{(Θ.sub.y*(j,k,i)}+αR, up-milling where α is a coefficient to adjust selecting range; step 6.2: determining lower and upper limits of the selected trajectory points in {right arrow over (Ox)} direction, i.e., x.sub.min=min{Θ.sub.x_trim*(j,k,i)} and x.sub.max=max{Θ.sub.x_trim*(j,k,i)}; discretize interval [x.sub.min,x.sub.max] equally with a step δx to obtain x coordinate set of interpolation points {x.sub.min, x.sub.min+δx,x.sub.min+2.Math.δx, . . . ,x.sub.max}, where number of points is denoted by N.sub.s; step 6.3: taking the selected points (Θ.sub.x_trim*(j,k,i), Θ.sub.y_trim*(j,k,i)) in step 6.1 as known points, using spline interpolation algorithm to determine y coordinates of interpolation points with regard to x coordinates {x.sub.min, x.sub.min+δx,x.sub.min+2.Math.δx, . . . ,x.sub.max}, which results in densified motion trajectories of cutting tool edges during one rotation period T of cutting tool, i.e., (x.sub.s(l), y.sub.s(l)),l=1, 2 . . . N.sub.s; step 6.4: based on periodicity of chatter-free milling dynamic responses, i.e., x.sub.s(l) and x.sub.s(l)+n.sub.revT share same y coordinate y.sub.s(l), obtaining densified motion trajectory points of cutting tool during n.sub.rev rotation periods of cutting tool, i.e., (x.sub.s_n(l), y.sub.s_n(l)), l=1, 2, . . . N.sub.s_n, where N.sub.s_n, is number of interpolation points on n.sub.rev rotation periods of cutting tool, and n.sub.rev≥4; step 6.5: selecting densified motion trajectory points of cutting tool edges on middle periods that satisfy x.sub.min+T≤x.sub.s_n(l)≤n.sub.rev.Math.x.sub.max−T to construct a new densified point set (x.sub.s_n_trim(l), y.sub.s_n_trim(l), l=1, 2, . . . N.sub.s_n_trim, where N.sub.s_n_trim is number of trimmed interpolation points; step 6.6: sorting the selected set (x.sub.s_n_trim(l), y.sub.s_n_trim(l)), in order of axial height to construct N.sub.a sub-sets of densified points, i.e., (x.sub.s_n_trim(j,l), y.sub.s_n_trim(j,l)), j=1,2, . . . ,N.sub.a; step 6.7: at each axial height, comparing values of v coordinates for all teeth that have same x coordinate; based on following equation, selecting densified points that are nearest to milled surface of workpiece to construct final surface topography points of workpiece (x.sub.surf(j,l), y.sub.surf (j,l)), l=1, 2, . . . N.sub.s_n_trim: $\begin{array}{cc}{y}_{\mathrm{surf}}\left(j,l\right)=\underset{k}{\max }\left\{y_{\mathrm{s\_n}\mathrm{\_trim}}\left(j,l\right)\right\},& \mathrm{down}-\mathrm{milling}\\ y_{\mathrm{surf}}\left(j,l\right)=\underset{k}{\min }\left\{y_{\mathrm{s\_n}\mathrm{\_trim}}\left(j,l\right)\right\},& \mathrm{up}-\mathrm{milling}\end{array}$

8. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 7 is: calculating surface location error (SLE) based on motion trajectories Θ.sub.y*(i,j,k) in normal direction of workpiece; surface location error SLE(j,k) with regard to cutting element (j,k) is: $\mathrm{SLE}\left(j,k\right)=\{\begin{array}{cc}\underset{i}{\max }\left\{{\Theta }_y^{*}\left(j,k,i\right).Math.1\le i\le 2i_m+2\right\}-D/2,& \mathrm{down}-\mathrm{milling}\\ \underset{i}{-\min }\left\{{\Theta }_y^{*}\left(j,k,i\right).Math.1\le i\le 2i_m+2\right\}-D/2,& \mathrm{up}-\mathrm{milling}\end{array}$ surface location error with regard to axial height z.sub.j is: $\mathrm{SLE}\left(j\right)=\{\begin{array}{cc}\underset{i,k}{\max }\left\{{\Theta }_y^{*}\left(j,k,i\right).Math.1\le i\le 2i_m+2,1\le k\le N\right\}-D/2,& \mathrm{down}-\mathrm{milling}\\ \underset{i,k}{-\min }\left\{{\Theta }_y^{*}\left(j,k,i\right).Math.1\le i\le 2i_m+2,1\le k\le N\right\}-D/2,& \mathrm{up}-\mathrm{milling}\end{array}$ where positive value mean overcut, and negative value mean undercut.

9. The method for simulating chatter-free milled surface topography of claim 1, wherein the step 8 is: based on points (x.sub.surf(j,l), y.sub.surf(j,l), y.sub.surf(j,l)) that participate in generating surface topography of workpiece, calculating surface roughness: $\mathrm{Ra}\left(j\right)=\frac{1}{N_{\mathrm{s\_n}\mathrm{\_trim}}}\underset{i^{\prime }=1}{\overset{N_{\mathrm{s\_n}\mathrm{\_trim}}}{.Math.}}.Math.y_{\mathrm{surf}}\left(j,l\right)-{\stackrel{\_}{y}}_{\mathrm{surf}}\left(j,l\right).Math..$

Description

BRIEF DESCRIPTION OF DRAWINGS

[0058] Other features, objects and advantages of the present invention will become apparent upon reading following detailed description of unrestricted examples, while referring to attached drawings in which:

[0059] FIGS. 1(a) to 1(d) are a flow diagram representing simulation procedure of milled surface topography of workpiece, where FIG. 1(a) represents motion trajectories of cutting tool edges, FIG. 1(b) represents densified motion trajectories of cutting tool using spline interpolation for each tooth during one rotation period of cutting tool, FIG. 1(c) represents densified motion trajectories of cutting tool during several rotation periods of cutting tool, and FIG. 1(d) represents the final surface topography of workpiece, which is determined by comparing motion trajectories of cutting tool edges for different teeth.

[0060] FIGS. 2(a) to 2(b) are comparison between experimental surface topography and simulated results; FIG. 2(a) is microphotograph of milled surface topography, and FIG. 2(b) is simulated milled surface topography.

[0061] FIG. 3 is comparison between simulated and measured surface location errors along tool axial height.

[0062] FIG. 4 is simulated surface roughness value along tool axial height.

EXAMPLES

[0063] The present invention will be described in detail below in conjunction with specific embodiments. The following examples will help those skilled in the art to further understand the present invention, but do not limit the present invention in any way. It should be pointed out that for those of ordinary skill in the art, several modifications and improvements can be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0064] Refer to FIG. 1(a), FIG. 1(b), FIG. 1(c), FIG. 1(d), FIG. 2(a) and FIG. 2(b).

[0065] This example provides a method for simulating chatter-free milled surface topography, comprising steps of:

[0066] Step 1 comprises sub-steps of:

[0067] step 1.1: considering flexibility of both cutting tool and workpiece, variance of pitch angles and helix angles of cutting tool, and teeth runout, perform dynamics modeling on milling system and build a delay differential equation with multiple delays to represent the dynamic model in physical space:

[00021]$\begin{array}{cc}M\stackrel{.Math.}{q}\left(t\right)+C\dot{q}\left(t\right)+Kq\left(t\right)=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{K_c\left(t,j,k\right)\left[q\left(t\right)-q\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+F_0\left(t,j,k\right)\right\}& \left(1\right)\end{array}$

where M, C and K are mass, damping and stiffness matrices, respectively; {umlaut over (q)}(t) {dot over (q)}(t) and q(t) are acceleration, velocity and displacement vectors at time t, respectively; K.sub.c(t,j,k) is directional coefficient matrix for Tooth k at j-th axial slice at time t; F.sub.0(t,j,k) is cutting force vector for Tooth k at j-th axial slice at time t, which is independent of chip regeneration; Σ.sub.p=k.sub.v.sup.kτ(j,p) is time delay with regard to cutting element (j,k), Σ is mathematical summation operator, p is a variable for summation, and k.sub.v is tooth serial number to start delay calculation; N is number of teeth; and N.sub.a is number of axial discretization of tool;

[0068] step 1.2: perform modal coordinate transformation on the dynamic model in step 1.1, which can be transformed from physical space into modal space; the modal coordinate transformation equation is:

q(t)=PΓ(t)  (2)

where P is the mode shape matrix, and Γ(t) is the modal displacement at time t;

[0069] dynamics model in modal space which is represented with delay differential equation with multiple delays turns out to be:

[00022]$\begin{array}{cc}{M}_{\Gamma }\stackrel{.Math.}{\Gamma }\left(t\right)+C_{\Gamma }\stackrel{.}{\Gamma }\left(t\right)+K_{\Gamma }\Gamma \left(t\right)=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{P^T{K}_c\left(t,j,k\right)P\left[\Gamma \left(t\right)-\Gamma \left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+P^T{F}_0\left(t,j,k\right)\right\}& \left(3\right)\end{array}$

[0070] where M.sub.Γ, C.sub.Γ and K.sub.Γ are modal mass, modal damping and modal stiffness matrices, respectively; {umlaut over (Γ)}(t), {dot over (Γ)}(t) and Γ(t) are modal acceleration, modal velocity and modal displacement vectors, respectively;

[0071] for milling system with flexible cutting tool and flexible workpiece, the dynamic model becomes:

[00023]$\begin{array}{cc}\left[\begin{array}{cc}{M}_{\Gamma ;T}& \\ & M_{\Gamma ;W}\end{array}\right]\left\{\begin{array}{c}{\stackrel{.Math.}{\Gamma }}_T\left(t\right)\\ {\stackrel{.Math.}{\Gamma }}_W\left(t\right)\end{array}\right\}+\left[\begin{array}{cc}{C}_{\Gamma ;T}& \\ & C_{\Gamma ;W}\end{array}\right]\left\{\begin{array}{c}{\stackrel{.}{\Gamma }}_T\left(t\right)\\ {\stackrel{.}{\Gamma }}_W\left(t\right)\end{array}\right\}+\left[\begin{array}{cc}{K}_{\Gamma ;T}& \\ & K_{\Gamma ;W}\end{array}\right]\left\{\begin{array}{c}{\Gamma }_T\left(t\right)\\ {\Gamma }_W\left(t\right)\end{array}\right\}=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{\left[\begin{array}{c}{P}_T^T\\ -P_W^T\end{array}\right]K_c\left(t,j,k\right)\left[\begin{array}{cc}{P}_T& -P_W\end{array}\right]\left\{\begin{array}{c}{\Gamma }_T\left(t\right)-{\Gamma }_T\left({.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\\ {\Gamma }_W\left(t\right)-{\Gamma }_W\left({.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\end{array}\right\}+\left[\begin{array}{c}{P}_T^T\\ -P_W^T\end{array}\right]F_0\left(t,j,k\right)\right\}& \left(4\right)\end{array}$

where M.sub.Γ,T, C.sub.Γ,T and K.sub.Γ,T are modal mass, modal damping and modal stiffness matrices of cutting tool; M.sub.Γ;W, C.sub.Γ;W and K.sub.Γ;W are modal mass, modal damping and modal stiffness matrices of workpiece; {umlaut over (Γ)}.sub.T(t), {dot over (Γ)}.sub.T(t) and Γ.sub.T(t) are modal acceleration, modal velocity and modal displacement vectors of cutting tool; {umlaut over (Γ)}.sub.W(t), {dot over (Γ)}.sub.W(t) and Γ.sub.W(t) are modal acceleration, modal velocity and modal displacement vectors of workpiece; P.sub.T and P.sub.W are mode shape matrices of cutting tool and workpiece, respectively; P.sub.T.sup.T and P.sub.W.sup.T are transposed matrices of P.sub.T and P.sub.W, respectively;

[0072] for milling system with flexible cutting tool and rigid workpiece, the dynamic model becomes:

[00024]$\begin{array}{cc}{M}_{\Gamma ;T}{\stackrel{.Math.}{\Gamma }}_T\left(t\right)+C_{\Gamma ;T}{\dot{\Gamma }}_T\left(t\right)+K_{\Gamma ;T}{\Gamma }_T\left(t\right)=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{P_T^T{K}_c\left(t,j,k\right)\left\{P_T\left[{\Gamma }_T\left(t\right)-{\Gamma }_T\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]\right\}+P_T^T{F}_0\left(t,j,k\right)\right\}& \left(5\right)\end{array}$

[0073] for milling system with rigid cutting tool and flexible workpiece, the dynamic model becomes:

[00025]$\begin{array}{cc}{M}_{\Gamma ;W}{\stackrel{.Math.}{\Gamma }}_W\left(t\right)+C_{\Gamma ;W}{\dot{\Gamma }}_W\left(t\right)+K_{\Gamma ;W}{\Gamma }_W\left(t\right)=\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{P_W^T{K}_c\left(t,j,k\right)\left\{P_W\left[{\Gamma }_W\left(t\right)-{\Gamma }_W\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]\right\}-P_W^T{F}_0\left(t,j,k\right)\right\}& \left(6\right)\end{array}$

[0074] Step 2 comprises sub-steps of:

[0075] step 2.1: perform state space transformation on the dynamic model in modal space by introducing state variable vector x(t)=[(Γ(t)).sup.T, ({dot over (Γ)}(t)).sup.T].sup.T, and delay differential equation with multiple delays in state space becomes:

[00026]$\begin{array}{cc}\dot{x}\left(t\right)=Ax\left(t\right)+\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{B\left(t,j,k\right)\left[x\left(t\right)-x\left(t-{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+D\left(t,j,k\right)\right\}& \left(7\right)\end{array}$

where x(t) denotes state variable vector, {dot over (x)}(t) denotes first derivative of x(t) with regard to time t,

[00027]$A=\left[\begin{array}{cc}0& I\\ -M_{\Gamma }^{-1}{K}_{\Gamma }& -M_{\Gamma }^{-1}{C}_{\Gamma }\end{array}\right],$$B\left(t,j,k\right)=\left[\begin{array}{cc}0& 0\\ M_{\Gamma }^{-1}{P}^T{K}_c\left(t,j,k\right)P& 0\end{array}\right],$$D\left(t,j,k\right)=\left[\begin{array}{c}0\\ M_{\Gamma }^{-1}{P}^T{F}_0\left(t,j,k\right)\end{array}\right],$

and I is unit matrix;

[0076] step 2.2: divide one rotation period of cutting tool T into 2i.sub.m+2 equal intervals with a series of discrete nodes {t.sub.0,t.sub.1, . . . ,t.sub.2i.sub.m.sub.+2}; analytical expression of dynamic responses at interval [t.sub.i,t] with regard to the delay differential equation with multiple delays in step 2.1 is:

[00028]$\begin{array}{cc}x\left(t\right)=e^{A\left(t-t_i\right)}x\left(t_i\right)+\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}{\int }_{t_i}^t\left\{e^{A\left(t-\xi \right)}B\left(\xi ,j,k\right)\left[x\left(\xi \right)-x\left(\xi -{.Math.}_{p=k_v}^k\tau \left(j,p\right)\right)\right]+e^{A\left(t-\xi \right)}D\left(\xi ,j,k\right)\right\}d\xi .& \left(8\right)\end{array}$

[0077] for brevity, define x.sub.i:=x(t.sub.i), E(t,t.sub.i):=e.sup.A(t−t.sup.i.sup.)x(t.sub.i), H.sub.j,k(t,ξ,x(ξ)):=e.sup.A(t−ξ)B(ξ, j,k)[x(ξ)−x(ξ−τ(j,k, p))] and J.sub.j,k(t,ξ):=e.sup.A(t−ξ)D(ξ,j, j,k);

[0078] on sub-interval [t.sub.2i,t.sub.2i+2], x(t) is approximated using Simpson's rule:

[00029]$\begin{array}{cc}{x}_{2i+2}=E\left(t_{2i+2},t_{2i}\right)+\frac{h}{3}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{H_{j,k}\left(t_{2i+2},t_{2i},x_{2i}\right)+4H_{j,k}\left(t_{2i+2},t_{2i+1},x_{2i+1}\right)+H_{j,k}\left(t_{2i+2},t_{2i+2},x_{2i+2}\right)+J_{j,k}\left(t_{2i+2},t_{2i}\right)+4J_{j,k}\left(t_{2i+2},t_{2i+1}\right)+J_{j,k}\left(t_{2i+2},t_{2i+2}\right)\right\}& \left(9\right)\end{array}$

[0079] on sub-interval [t.sub.2i,t.sub.2i+1], x(t) is approximated using fourth-order Runge-Kutta algorithm:

[00030]$\begin{array}{cc}{x}_{2i+1}=E\left(t_{2i+1},t_{2i}\right)+\frac{h}{6}\underset{k=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{N_a}{.Math.}}\left\{H_{j,k}\left(t_{2i+1},t_{2i},x_{2i}\right)+2H_{j,k}\left(t_{2i+1},t_{2i+1/2},x_{2i+1/2}\right)+2H_{j,k}\left(t_{2i+1},t_{2i+1/2},x_{2i+1/2}\right)+H_{j,k}\left(t_{2i+1},t_{2i+1},x_{2i+1}\right)+J_{j,k}\left(t_{2i+1},t_{2i}\right)+2J_{j,k}\left(t_{2i+1},t_{2i+1/2}\right)+2J_{j,k}\left(t_{2i+1},t_{2i+1/2}\right)+J_{j,k}\left(t_{2i+1},t_{2i+1}\right)\right\}& \left(10\right)\end{array}$