CN102004818B - 自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法 - Google Patents

Abstract

本发明公开了一种自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法，用于解决现有的自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法计算工作量大的技术问题。技术方案是以理论刀轨上的等效刀具位置点对应的等效进给方向、等效法向量和等效曲率代替实际刀轨上的刀具位置点对应的进给方向、法向量和曲率，避免用非二阶连续的实际刀轨近似理论刀轨时所导致的预测铣削力的突变的现象；在考虑偏心的同时，推导出瞬时未变形切屑厚度的解析模型，无需用数值计算的方法计算瞬时未变形切屑厚度。利用本发明的方法，用MATLAB在个人计算机(Intel Core(TM)2Duo Processor，2.4GHz，2GB)计算10个周期的瞬时铣削力的时间为1.9s，比采用文献2的方法计算时间233.3s，计算效率提高122.79倍。

Description

自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法

技术领域

本发明涉及一种圆周铣削过程中瞬时铣削力建模方法，特别是一种自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法。

背景技术

圆周铣削是通过立铣刀在材料表面切除多余的材料层来获得理想的工件形状、尺寸以及表面光洁度的机械加工方法，它是机械制造中加工各种模具、汽车零件以及航空航天零件的最常用的方式之一。在圆周铣削过程中，经常出现刀具磨损、刀具破损、加工颤振以及加工变形等诸多问题，严重影响零件的加工精度和加工效率。为了克服以上问题，需要建立圆周铣削加工中铣削力预测模型。诸多学者已经在具有代表性的平面铣削方面做了广泛而深入的研究。由于自由曲线外形零件的圆周铣削在实际加工中也经常遇到，一些学者在这方面做了相关的研究。

文献1“Z.C.Wei，M.J.Wang，R.G.Ma，L.Wang，Modeling of process geometry in peripheralmilling of curved surfaces，Journal of Materials Processing Technology 210(2010)799-806.”公开了一种自由曲线外形零件圆周铣削过程中铣削力的模型，并系统地给出了铣削力预测的整个流程。然而该模型忽略了刀具偏心对铣削力的影响，并且无法避免利用非二阶连续的实际刀轨预测铣削力时所导致的预测铣削力的突变。

文献2“K.A.Desai，P.K.Agarwal，P.V.M.Rao，Process geometry modeling with cutter runoutfor milling of curved surfaces，International Journal of Machine Tools and Manufacture 49(2009)1015-1028.”公开了一种自由曲线外形零件圆周铣削过程中铣削力的模型，并且考虑了刀具偏心对铣削力的影响。该模型中计算与铣削力相关的瞬时未变形切屑厚度的基本方法如下：

1.计算当前切削刀齿的实际半径值和前m个刀齿的实际半径值。

2.计算直线与圆弧的交点。其中，该直线由当前的切触点和与当前切触点相对应的刀齿位置角确定，圆弧是以对应于当前切触点的前m个刀齿轨迹上的切触点所对应的刀具位置点为圆心，前m个刀齿的实际半径值为半径。

3.计算以上交点到当前刀具位置点的距离。

4.计算可能的瞬时未变形切屑厚度，即以当前刀齿的实际半径值减去交点到当前刀具位置点的距离。

5.重复步骤1-4，直到m等于刀齿数N_f结束。

6.选出N_f个可能的瞬时未变形切屑厚度的最小值，并将该最小值与0比较，取大者作为瞬时未变形切屑厚度。

从以上步骤可以看出，该方法中计算可能的瞬时未变形切削厚度是通过数值计算的方法完成的。这种方法求一次瞬时未变形切屑厚度需要求解N_f个非线性方程，因此计算工作量大。如，用MATLAB在个人计算机(Intel Core(TM)2Duo Processor，2.4GHz，2GB)计算10个周期的瞬时铣削力的时间为233.3s。

发明内容

为了克服现有的自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法计算工作量大的不足，本发明提供一种自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法，该方法首先通过采样时间间隔计算出各时刻所对应的刀具位置点，同时计算出每个刀具位置点对应的几何参数，然后通过数学推导建立局部坐标系下切入/切出角的数值计算模型和瞬时未变形切屑厚度的解析模型，最后建立局部坐标系下瞬时铣削力与瞬时未变形切屑厚度的关系，并将局部坐标系下的瞬时铣削力投影到整体坐标系下，可以减少自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法的计算工作量。

本发明解决其技术问题所采用的技术方案是：一种自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法，其特点是包括下述步骤：

(1)选定立铣刀和工件几何参数，包括立铣刀的半径RAD、螺旋角β₀、刀齿数N_f、刀具偏心参数ρ和λ，得到工件表面矢量方程p(u)＝[X(u)，Y(u)，Z(u)]^T和工件毛坯边界的矢量方程p_W(v)＝[X_W(v)，Y_W(v)，Z_W(v)]^T；设定切削参数，包括进给率V_f、轴向切削深度R_z、径向切削深度R_r、刀具主轴转速n；输入由CAM系统产生的加工工件的实际刀轨信息和采样时间间隔T_s；

(2)根据步骤(1)设定的采样时间间隔T_s和实际刀轨信息，按照下式

在实际刀轨上计算出各采样时刻对应的刀具位置点；

式中，t表示采样时刻；

p_en表示直线刀轨段的起点；p_st表示直线刀轨段的终点；o_CTP表示圆弧刀轨段的圆心；R_CTP表示圆弧刀轨段的半径；

(3)根据步骤(1)给定的工件表面矢量方程，通过偏置，计算出理论刀轨的矢量方程p_e(u)＝[X_t(u)，Y_t(u)，Z_t(u)]^T，

式中，

$X_t\left(u\right)=X\left(u\right)+\frac{{\mathrm{dY}}^{\′}\left(u\right)}{\sqrt{{\left[X^{\′}\left(u\right)\right]}^2+{\left[Y^{\′}\left(u\right)\right]}^2}};$

$Y_t\left(u\right)=Y\left(u\right)+\frac{{\mathrm{dX}}^{\′}\left(u\right)}{\sqrt{{\left[X^{\′}\left(u\right)\right]}^2+{\left[Y^{\′}\left(u\right)\right]}^2}}$

Z_t(u)等于实际刀轨信息中刀心点的Z坐标值；

(4)按照下式，将步骤(2)中得到的刀具位置点p_a(t)投影到理论刀轨上，得到等效刀具位置点p_e(u(t))，

(p_e(u(t))-p_a(t))×n_a(t)＝0

解非线性方程组，得到参数u(t)，将参数u(t)带入理论刀轨矢量方程p_e(u)中，即得到等效刀具位置点p_e(u(t))；

(5)将步骤(4)中解得的参数u(t)代入下式，即得到等效进给方向、等效法向量以及等效曲率，以等效进给方向、等效法向量以及等效曲率作为刀具位置点处的实际进给方向、实际法向量以及实际曲率；

f_e(u(t))＝[X′_t(t)  Y′_t(t)  0]^T

n_e(u(t))＝[0  0  1]^T×f_e(u(t))

$K_e\left(u\left(t\right)\right)=\frac{X_t^{\′}\left(u\left(t\right)\right)Y_t^{\′\′}\left(u\left(t\right)\right)-X_t^{\′\′}\left(u\left(t\right)\right)Y_t^{\′}\left(u\left(t\right)\right)}{{({\left(X_t^{\′}\left(u\left(t\right)\right)\right)}^2+{\left(Y_t^{\′}\left(u\left(t\right)\right)\right)}^2)}^{\frac{3}{2}}}$

同时计算出等效进给方向在整体坐标系中的角位置；

$\mathrm{\θ}\left(t\right)=\arccos \left(\frac{f_e\left(u\left(t\right)\right)\·i}{\left|f_e\left(u\left(t\right)\right)\right|}\right)$

式中，i＝[1  0  0]^T；

(6)将刀具参与切削的区域沿轴向划分为N个等高梁段，根据步骤(1)给定的工件毛坯边界矢量方程、立铣刀几何参数和步骤(2)中得到的刀具位置点，通过下式计算刀刃片{i，j}的切入切出角和切出角；

(a)切入角：

式中，p_W(v_i，j，en(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点，

${\mathrm{RAD}}_{i,j}=\mathrm{RAD}+\mathrm{\ρ}\cos [\mathrm{\λ}-\frac{\tan \left({\mathrm{\β}}_0\right)}{\mathrm{RAD}}z-\frac{2(i-1)\mathrm{\π}}{N_f}],i=\mathrm{1,2},\·\·\·N_f,$ N_f是刀刃数，j＝1，2，…，N，z是刀刃片{i，j}中点的Z向高度；

(b)切出角：

刀具切入工件阶段；

刀具在其他切削阶段；

式中，

p_D(t，m)为满足方程组

的点，

m表示当前刀齿之前的第m个刀齿，m＝1，…，N_f

p_W(v_{i，j，B，ex}(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点；

(7)通过下式计算作用在刀刃片{i，j}的铣削力：

式中，K_i，j，T(t)，K_i，j，R(t)分别是与h_i，j(t)相关的切向合径向瞬时铣削力系数，

Δa是等高梁段的高度，

h_i，j(t，m)是利用步骤(5)中得到的等效曲率K_e(u(t))，通过下式计算

凸型曲面；

凹型曲面；

式中， $R_{\mathrm{TP}}\left(t\right)=\left|\frac{1}{K_e\left(u\left(t\right)\right)}\right|,$ $f_t=\frac{V_f}{n{N}_f};$

(8)将各个刀刃上的力转化到X_S，Y_S和Z_S方向：

式中，

是刀具在t时刻处与刀刃片{i，j}对应的切削角度，被定义为从Y_S向顺时针到刀刃片{i，j}的中点所转过的角度；

(9)对于每个侧刃，将作用在所有刀刃片上的微元力求和，求得t时刻局部坐标系下作用于各个侧刃的铣削合力：

$F_{X_S}\left(t\right)=\underset{i,j}{\mathrm{\Σ}}{F}_{i,j,X_S}\left(t\right)$

$F_{Y_S}\left(t\right)=\underset{i,j}{\mathrm{\Σ}}{F}_{i,j,Y_S}\left(t\right)$

(10)将局部坐标系下的铣削合力转化到X，Y和Z方向：

$F_X\left(t\right)=F_{X_S}\left(t\right)\cos \mathrm{\θ}\left(t\right)-F_{Y_S}\left(t\right)\sin \mathrm{\θ}\left(t\right)$

$F_Y\left(t\right)=F_{X_S}\left(t\right)\sin \mathrm{\θ}\left(t\right)+F_{Y_S}\left(t\right)\cos \mathrm{\θ}\left(t\right)$

本发明的有益效果是：以理论刀轨上的等效刀具位置点对应的等效进给方向、等效法向量和等效曲率代替实际刀轨上的刀具位置点对应的进给方向、法向量和曲率，避免用非二阶连续的实际刀轨近似理论刀轨时所导致的预测铣削力的突变的现象；在考虑偏心的同时，推导出瞬时未变形切屑厚度的解析模型，无需用数值计算的方法计算瞬时未变形切屑厚度。利用本发明的方法，用MATLAB在个人计算机(Intel Core(TM)2Duo Processor，2.4GHz，2GB)计算10个周期的瞬时铣削力的时间为1.9s，比采用文献2的方法计算时间233.3s，计算效率提高122.79倍。

下面结合附图和具体实施方式对本发明进一步说明。

附图说明

图1是实施例1的铣削力随时间变化图。

图2是图1中3.7-3.75s时的放大图。

图3是图1中8.6-8.65s时的放大图。

图4是实施例1的刀轨曲率随时间变化图和进给方向的角位置θ(t)随时间变化图。

图5是实施例2的铣削力随时间变化图。

图6是图5中3.2-3.25s时的放大图。

图7是图5中5.7-5.75s时的放大图。

图中，实线代表测量力，点画线代表预测力，点线代表切入阶段，虚线代表稳定切削阶段。

具体实施方式

实施例1：(1)选定半径RAD＝6mm、螺旋角β₀＝30度的三齿硬质合金立铣刀在三坐标立铣床上对铝合金Al7050进行顺铣切削，刀具偏心参数ρ＝0.0026和λ＝31.8°，主轴转速n＝2000RPM，进给率V_f＝300mm/min，轴向切削深度R_z＝10mm，径向切削深度R_r＝3mm，采样时间间隔T_s＝0.0002s，输入工件欲得到表面的矢量方程p(u)＝[X(u)，Y(u)，Z(u)]^T，式中

X(u)＝20+105u(1-u)²+15u²(1-u)+40u³

Y(u)＝5+90u²(1-u)+30u³        u∈[0，1]

0≤Z(u)≤10

毛坯边界的矢量方程p_W(v)＝[X_W(v)，Y_W(v)，Z_W(v)]^T，式中

$X_W\left(v\right)=20+105v{(1-v)}^2+{15v}^2(1-v)+{40v}^3-\frac{36v(1-v)}{\sqrt{49-364v+{1184v}^2-{1640v}^3+{820v}^4}}$

$Y_W\left(v\right)=5+{90v}^2(1-v)+{30v}^3+\frac{21{(1-v)}^2-36v(1-v)+{21v}^2}{\sqrt{49-364v+{1184v}^2-{1640v}^3+{820v}^4}}$ v∈[0，1]

0≤Z_W(v)≤10

X_W(v)＝20

Y_W(v)＝8(1+v)-5v v∈[-1，0]

0≤Z_W(v)≤10

X_W(v)＝60

Y_W(v)＝38(v-1)+35(2-v)v∈[1，2]

0≤Z_W(v)≤10

实际刀轨信息

…

GOT0/13.9002，16.9496，0.0000

CIRCLE/19.9000，16.9996，0.0000，0.0000000，0.0000000，-1.0000000，6.0000，0.0100，0.5000，12.0000，0.0000

GOTO/19.9500，10.9998，0.0000

CIRCLE/20.2133，52.2325，0.0000，0.0000000，0.0000000，-1.0000000，41.2336，0.0100，0.5000，12.0000，0.0000

GOTO/24.7369，11.2478，0.0000

CIRCLE/22.4035，32.7858，0.0000，0.0000000，0.0000000，-1.0000000，21.6640，0.0100，0.5000，12.0000，0.0000

GOTO/28.0904，11.8815，0.0000

CIRCLE/25.3711，22.0905，0.0000，0.0000000，0.0000000，-1.0000000，10.5649，0.0100，0.5000，12.0000，0.0000

GOTO/30.4537，12.8285，0.0000

CIRCLE/28.0568，17.2549，0.0000，0.0000000，0.0000000，-1.0000000，5.0337，0.0100，0.5000，12.0000，0.0000

GOTO/32.2469，14.4654，0.0000

CIRCLE/25.8212，18.5149，0.0000，0.0000000，0.0000000，-1.0000000，7.5952，0.0100，0.5000，12.0000，0.0000

GOTO/33.0908，16.3149，0.0000

GOTO/33.2840，16.9667，0.0000

GOTO/33.4579，17.6565，0.0000

GOTO/33.6215，18.4079，0.0000

GOTO/33.7774，19.2173，0.0000

GOTO/34.4212，22.9586，0.0000

GOTO/34.6224，24.0036，0.0000

GOTO/34.8583，25.0871，0.0000

CIRCLE/58.6472，19.7687，0.0000，0.0000000，0.0000000，1.0000000，24.3761，0.0100，0.5000，12.0000，0.0000

GOTO/36.3976，29.7262，0.0000

CIRCLE/52.0438，22.6689，0.0000，0.0000000，0.0000000，1.0000000，17.1642，0.0100，0.5000，12.0000，0.0000

GOTO/40.1300，35.0248，0.0000

CIRCLE/52.4252，22.0901，0.0000，0.0000000，0.0000000，1.0000000，17.8460，0.0100，0.5000，12.0000，0.0000

GOTO/44.3945，38.0271，0.0000

CIRCLE/55.0895，16.6900，0.0000，0.0000000，0.0000000，1.0000000，23.8675，0.0100，0.5000，12.0000，0.0000

GOTO/49.4660，39.8855，0.0000

CIRCLE/58.0277，4.7902，0.0000，0.0000000，0.0000000，1.0000000，36.1246，0.0100，0.5000，12.0000，0.0000

GOTO/54.8228，40.7723，0.0000

CIRCLE/59.7855，-13.1822，0.0000，0.0000000，0.0000000，1.0000000，54.1822，0.0100，0.5000，12.0000，0.0000

GOTO/59.9500，40.9998，0.0000

GOTO/60.0000，41.0000，0.0000

(2)根据步骤(1)设定的采样时间间隔T_s＝0.0002s，按照下式，在实际刀轨上计算出刀具位置点

式中

t表示采样时刻

$f_a\left(t\right)=\frac{p_{\mathrm{en}}-p_{\mathrm{st}}}{|p_{\mathrm{en}}-p_{\mathrm{st}}|}$

$\left[R\right]=\left[\begin{array}{ccc}{\cos \mathrm{\α}}_s& -\sin {\mathrm{\α}}_s& 0\\ {\sin \mathrm{\α}}_s& \cos {\mathrm{\α}}_s& 0\\ 0& 0& 1\end{array}\right]$

${\mathrm{\α}}_s=\frac{V_f{T}_s}{R_{\mathrm{CTP}}}$

p_en和p_st分别表示直线刀轨段的起点和终点，o_CTP和R_CTP分别表示圆弧刀轨段的圆心和半径。

(3)根据步骤(1)给定的欲得到表面的矢量方程，通过偏置，计算出理论刀轨的矢量方程p_e(u)＝[X_t(u)，Y_t(u)，Z_t(u)]^T，式中

$X_t\left(u\right)=20+105u{(1-u)}^2+{15u}^2(1-u)+{40u}^3-\frac{72u(1-u)}{\sqrt{49-364u+{1184u}^2-{1640u}^3+{820u}^4}}$

$Y_t\left(u\right)=5+{90u}^2(1-u)+{30u}^3+\frac{42{(1-u)}^2-72u(1-u)+{42u}^2}{\sqrt{49-364u+{1184u}^2-{1640u}^3+{820u}^4}}$ u∈[0，1]。

Z_t(u)＝0

(4)按照下式，将步骤(2)中得到的刀具位置点p_a(t)投影到理论刀轨上，得到等效刀具位置点，

(p_e(u(t))-p_a(t))×n_a(t)＝0

解非线性方程组，得到参数u(t)，将参数u(t)带入理论导轨矢量方程p_e(u)中，即得到等效刀具位置点p_e(u(t))。

(5)将步骤(4)中解得的参数u(t)代入下式，即得到等效进给方向、等效法向量以及等效曲率，以等效进给方向、等效法向量以及等效曲率作为刀具位置点处的实际进给方向、实际法向量以及实际曲率。

f_e(u(t))＝[X′_t(t)  Y′_t(t)  0]^T

n_e(u(t))＝[0  0  1]^T×f_e(u(t))

$K_e\left(u\left(t\right)\right)=\frac{X_t^{\′}\left(u\left(t\right)\right)Y_t^{\′\′}\left(u\left(t\right)\right)-X_t^{\′\′}\left(u\left(t\right)\right)Y_t^{\′}\left(u\left(t\right)\right)}{{({\left(X_t^{\′}\left(u\left(t\right)\right)\right)}^2+{\left(Y_t^{\′}\left(u\left(t\right)\right)\right)}^2)}^{\frac{3}{2}}}$

同时计算出等效进给方向在整体坐标系中的角位置

$\mathrm{\θ}\left(t\right)=\arccos \left(\frac{f_e\left(u\left(t\right)\right)\·i}{\left|f_e\left(u\left(t\right)\right)\right|}\right)$

式中，i＝[1  0  0]^T。

(6)将刀具参与切削的区域沿轴向划分为N个等高梁段，根据步骤(1)给定的毛坯的边界的矢量方程、立铣刀几何参数和步骤(2)中得到的刀具位置点，通过下式计算刀刃片{i，j}的切入切出角。

(a)切入角的计算

式中p_W(v_i，j，en(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点，

${\mathrm{RAD}}_{i,j}=\mathrm{RAD}+\mathrm{\ρ}\cos [\mathrm{\λ}-\frac{\tan \left({\mathrm{\β}}_0\right)}{\mathrm{RAD}}z-\frac{2(i-1)\mathrm{\π}}{N_f}],i=\mathrm{1,2},\·\·\·N_f,$ N_f是刀刃数，j＝1，2，…，N，z是刀刃片{i，j}中点的Z向高度。

(b)切出角的计算

刀具切入工件阶段

刀具在其他切削阶段

式中，

p_D(t，m)为满足方程组

的点，

m表示当前刀齿之前的第m个刀齿，m＝1，…，N_f

p_W(v_{i，j，B，ex}(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点。

(7)通过下式计算作用在刀刃片{i，j}的铣削力：

式中

K_i，j，T(t)，K_i，j，R(t)分别是与h_i，j(t)相关的切向合径向瞬时铣削力系数，

Δa是等高梁段的高度，

h_i，j(t，m)利用步骤5)中得到的等效曲率K_e(u(t))，通过下式计算

凸型曲面

凹型曲面

式中， $R_{\mathrm{TP}}\left(t\right)=\left|\frac{1}{K_e\left(u\left(t\right)\right)}\right|,$ $f_t=\frac{V_f}{n{N}_f}.$

(8)将各个刀刃上的力转化到X_S，Y_S和Z_S方向：

式中，

是刀具在t时刻处与刀刃片{i，j}对应的切削角度，被定义为从Y_S向顺时针到刀刃片{i，j}的中点所转过的角度。

(9)对于每个侧刃，将作用在所有刀刃片上的微元力求和，求得t时刻局部坐标系下作用于各个侧刃的铣削合力：

$F_{X_S}\left(t\right)=\underset{i,j}{\mathrm{\Σ}}{F}_{i,j,X_S}\left(t\right)$

$F_{Y_S}\left(t\right)=\underset{i,j}{\mathrm{\Σ}}{F}_{i,j,Y_S}\left(t\right)$

(10)将局部坐标系下的铣削合力转化到X，Y和Z方向：

$F_X\left(t\right)=F_{X_S}\left(t\right)\cos \mathrm{\θ}\left(t\right)-F_{Y_S}\left(t\right)\sin \mathrm{\θ}\left(t\right)$

$F_Y\left(t\right)=F_{X_S}\left(t\right)\sin \mathrm{\θ}\left(t\right)+F_{Y_S}\left(t\right)\cos \mathrm{\θ}\left(t\right)$

通过以上方法，即得到自由曲线外形零件圆周铣削过程中瞬时铣削力的模型。

从图1、2、3可以看出，本发明的方法有效地考虑刀具偏心对铣削力的影响，其结果与实际测量结果吻合，验证了本发明的有效性。

从图1、2、3、4可以看出，在稳定切削阶段，切削力、刀轨曲率和进给方向的角位置θ(t)随时间变化连续，没有文献1中用非二阶连续的实际刀轨近似理论刀轨时所导致的切削过程中的几何量和预测铣削力的突变的现象。

实施例2：(1)选定半径RAD＝6mm、螺旋角β₀＝30度的三齿硬质合金立铣刀在三坐标立铣床上对铝合金Al7050进行顺铣切削，刀具偏心参数ρ＝0.0026和λ＝31.8°，主轴转速n＝2000RPM，进给率V_f＝450mm/min，轴向切削深度R_z＝10mm，径向切削深度R_r＝3mm，采样时间间隔T_s＝0.0002s，输入工件欲得到表面的矢量方程p(u)＝[X(u)，Y(u)，Z(u)]^T，式中

X(u)＝20+105u(1-u)²+15u²(1-u)+40u³

Y(u)＝5+90u²(1-u)+30u³        u∈[0，1]

0≤Z(u)≤10

毛坯边界的矢量方程p_W(v)＝[X_W(v)，Y_W(v)，Z_W(v)]^T，式中

$X_W\left(v\right)=20+105v{(1-v)}^2+{15v}^2(1-v)+{40v}^3-\frac{36v(1-v)}{\sqrt{49-364v+{1184v}^2-{1640v}^3+{820v}^4}}$

$Y_W\left(v\right)=5+{90v}^2(1-v)+{30v}^3+\frac{21{(1-v)}^2-36v(1-v)+{21v}^2}{\sqrt{49-364v+{1184v}^2-{1640v}^3+{820v}^4}}$ v∈[0，1]

0≤Z_W(v)≤10

X_W(v)＝20

Y_W(v)＝8(1+v)-5v v∈[-1，0]

0≤Z_W(v)≤10

X_W(v)＝60

Y_W(v)＝38(v-1)+35(2-v)v∈[1，2]

0≤Z_W(v)≤10

实际刀轨信息

…

GOTO/13.9002，16.9496，0.0000

CIRCLE/19.9000，16.9996，0.0000，0.0000000，0.0000000，-1.0000000，6.0000，0.0100，0.5000，12.0000，0.0000

GOTO/19.9500，10.9998，0.0000

CIRCLE/20.2133，52.2325，0.0000，0.0000000，0.0000000，-1.0000000，41.2336，0.0100，0.5000，12.0000，0.0000

GOTO/24.7369，11.2478，0.0000

CIRCLE/22.4035，32.7858，0.0000，0.0000000，0.0000000，-1.0000000，21.6640，0.0100，0.5000，12.0000，0.0000

GOTO/28.0904，11.8815，0.0000

CIRCLE/25.3711，22.0905，0.0000，0.0000000，0.0000000，-1.0000000，10.5649，0.0100，0.5000，12.0000，0.0000

GOTO/30.4537，12.8285，0.0000

CIRCLE/28.0568，17.2549，0.0000，0.0000000，0.0000000，-1.0000000，5.0337，0.0100，0.5000，12.0000，0.0000

GOT0/32.2469，14.4654，0.0000

CIRCLE/25.8212，18.5149，0.0000，0.0000000，0.0000000，-1.0000000，7.5952，0.0100，0.5000，12.0000，0.0000

GOTO/33.0908，16.3149，0.0000

GOTO/33.2840，16.9667，0.0000

GOTO/33.4579，17.6565，0.0000

GOTO/33.6215，18.4079，0.0000

GOTO/33.7774，19.2173，0.0000

GOTO/34.4212，22.9586，0.0000

GOTO/34.6224，24.0036，0.0000

GOTO/34.8583，25.0871，0.0000

CIRCLE/58.6472，19.7687，0.0000，0.0000000，0.0000000，1.0000000，24.3761，0.0100，0.5000，12.0000，0.0000

GOTO/36.3976，29.7262，0.0000

CIRCLE/52.0438，22.6689，0.0000，0.0000000，0.0000000，1.0000000，17.1642，0.0100，0.5000，12.0000，0.0000

GOTO/40.1300，35.0248，0.0000

CIRCLE/52.4252，22.0901，0.0000，0.0000000，0.0000000，1.0000000，17.8460，0.0100，0.5000，12.0000，0.0000

GOTO/44.3945，38.0271，0.0000

CIRCLE/55.0895，16.6900，0.0000，0.0000000，0.0000000，1.0000000，23.8675，0.0100，0.5000，12.0000，0.0000

GOTO/49.4660，39.8855，0.0000

CIRCLE/58.0277，4.7902，0.0000，0.0000000，0.0000000，1.0000000，36.1246，0.0100，0.5000，12.0000，0.0000

GOTO/54.8228，40.7723，0.0000

CIRCLE/59.7855，-13.1822，0.0000，0.0000000，0.0000000，1.0000000，54.1822，0.0100，0.5000，12.0000，0.0000

GOTO/59.9500，40.9998，0.0000

GOTO/60.0000，41.0000，0.0000

(2)根据步骤(1)设定的采样时间间隔T_s＝0.0002s，按照下式，在实际刀轨上计算出刀具位置点

式中

t表示采样时刻

$f_a\left(t\right)=\frac{p_{\mathrm{en}}-p_{\mathrm{st}}}{|p_{\mathrm{en}}-p_{\mathrm{st}}|}$

$\left[R\right]=\left[\begin{array}{ccc}{\cos \mathrm{\α}}_s& -\sin {\mathrm{\α}}_s& 0\\ {\sin \mathrm{\α}}_s& \cos {\mathrm{\α}}_s& 0\\ 0& 0& 1\end{array}\right]$

${\mathrm{\α}}_s=\frac{V_f{T}_s}{R_{\mathrm{CTP}}}$

p_en和p_st分别表示直线刀轨段的起点和终点，o_CTP和R_CTP分别表示圆弧刀轨段的圆心和半径。

(3)根据步骤(1)给定的欲得到表面的矢量方程，通过偏置，计算出理论刀轨的矢量方程p_e(u)＝[X_t(u)，Y_t(u)，Z_t(u)]^T，式中

$X_t\left(u\right)=20+105u{(1-u)}^2+{15u}^2(1-u)+{40u}^3-\frac{72u(1-u)}{\sqrt{49-364u+{1184u}^2-{1640u}^3+{820u}^4}}$

$Y_t\left(u\right)=5+{90u}^2(1-u)+{30u}^3+\frac{42{(1-u)}^2-72u(1-u)+{42u}^2}{\sqrt{49-364u+{1184u}^2-{1640u}^3+{820u}^4}}$ u∈[0，1]。

Z_t(u)＝0

(4)按照下式，将步骤(2)中得到的刀具位置点p_a(t)投影到理论刀轨上，得到等效刀具位置点，

(p_e(u(t))-p_a(t))×n_a(t)＝0

解非线性方程组，得到参数u(t)，将参数u(t)带入理论导轨矢量方程p_e(u)中，即得到等效刀具位置点p_e(u(t))。

(5)将步骤(4)中解得的参数u(t)代入下式，即得到等效进给方向、等效法向量以及等效曲率，以等效进给方向、等效法向量以及等效曲率作为刀具位置点处的实际进给方向、实际法向量以及实际曲率。

f_e(u(t))＝[X′_t(t)  Y′_t(t)  0]^T

n_e(u(t))＝[0  0  1]^T×f_e(u(t))

$K_e\left(u\left(t\right)\right)=\frac{X_t^{\′}\left(u\left(t\right)\right)Y_t^{\′\′}\left(u\left(t\right)\right)-X_t^{\′\′}\left(u\left(t\right)\right)Y_t^{\′}\left(u\left(t\right)\right)}{{({\left(X_t^{\′}\left(u\left(t\right)\right)\right)}^2+{\left(Y_t^{\′}\left(u\left(t\right)\right)\right)}^2)}^{\frac{3}{2}}}$

同时计算出等效进给方向在整体坐标系中的角位置

$\mathrm{\θ}\left(t\right)=\arccos \left(\frac{f_e\left(u\left(t\right)\right)\·i}{\left|f_e\left(u\left(t\right)\right)\right|}\right)$

式中，i＝[1  0  0]^T。

(6)将刀具参与切削的区域沿轴向划分为N个等高梁段，根据步骤(1)给定的毛坯的边界的矢量方程、立铣刀几何参数和步骤2)中得到的刀具位置点，通过下式计算刀刃片{i，j}的切入切出角。

(a)切入角的计算：

式中p_W(v_i，j，en(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点，

${\mathrm{RAD}}_{i,j}=\mathrm{RAD}+\mathrm{\ρ}\cos [\mathrm{\λ}-\frac{\tan \left({\mathrm{\β}}_0\right)}{\mathrm{RAD}}z-\frac{2(i-1)\mathrm{\π}}{N_f}],i=\mathrm{1,2},\·\·\·N_f,$ N_f是刀刃数，j＝1，2，…，N，z是刀刃片{i，j}中点的Z向高度。

(b)切出角的计算：

刀具切入工件阶段

刀具在其他切削阶段

式中，

p_D(t，m)为满足方程组

的点，

m表示当前刀齿之前的第m个刀齿，m＝1，…，N_f

p_W(v_{i，j，B，ex}(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点。

(7)通过下式计算作用在刀刃片{i，j}的铣削力：

式中，K_i，j，T(t)，K_i，j，R(t)分别是与h_i，j(t)相关的切向合径向瞬时铣削力系数，

Δa是等高梁段的高度，

h_i，j(t，m)利用步骤5)中得到的等效曲率K_e(u(t))，通过下式计算

凸型曲面

凹型曲面

式中， $R_{\mathrm{TP}}\left(t\right)=\left|\frac{1}{K_e\left(u\left(t\right)\right)}\right|,$ $f_t=\frac{V_f}{n{N}_f}.$

(8)将各个刀刃上的力转化到X_S，Y_S和Z_S方向：

式中，

是刀具在t时刻处与刀刃片{i，j}对应的切削角度，被定义为从Y_S向顺时针到刀刃片{i，j}的中点所转过的角度。

(9)对于每个侧刃，将作用在所有刀刃片上的微元力求和，求得t时刻局部坐标系下作用于各个侧刃的铣削合力：

$F_{X_S}\left(t\right)=\underset{i,j}{\mathrm{\Σ}}{F}_{i,j,X_S}\left(t\right)$

$F_{Y_S}\left(t\right)=\underset{i,j}{\mathrm{\Σ}}{F}_{i,j,Y_S}\left(t\right)$

(10)将局部坐标系下的铣削合力转化到X，Y和Z方向：

$F_X\left(t\right)=F_{X_S}\left(t\right)\cos \mathrm{\θ}\left(t\right)-F_{Y_S}\left(t\right)\sin \mathrm{\θ}\left(t\right)$

$F_Y\left(t\right)=F_{X_S}\left(t\right)\sin \mathrm{\θ}\left(t\right)+F_{Y_S}\left(t\right)\cos \mathrm{\θ}\left(t\right)$

通过以上方法，即得到自由曲线外形零件圆周铣削过程中瞬时铣削力的模型。

从图5、6、7可以看出，本发明的方法有效地考虑刀具偏心对铣削力的影响，其结果与实际测量结果吻合，验证了本发明的有效性。

利用本发明的方法，用MATLAB在个人计算机(Intel Core(TM)2Duo Processor，2.4GHz，2GB)计算10个周期的瞬时铣削力的时间为1.9s，比采用文献2的方法计算时间233.3s，计算效率提高122.79倍。

Claims (1)

1.一种自由曲线外形零件圆周铣削过程中瞬时铣削力建模方法，其特征在于包括下述步骤：

(1)选定立铣刀和工件几何参数，包括立铣刀的半径RAD、螺旋角β₀、刀齿数N_f、刀具偏心参数ρ和λ，得到工件表面矢量方程p(u)＝[X(u)，Y(u)，Z(u)]^T和工件毛坯边界的矢量方程P_W(v)＝[X_W(v)，Y_W(v)，Z_W(v)]^T；设定切削参数，包括进给率V_f、轴向切削深度R_z、径向切削深度R_r、刀具主轴转速n；输入由CAM系统产生的加工工件的实际刀轨信息和采样时间间隔T_s；

(2)根据步骤(1)设定的采样时间间隔T_s和实际刀轨信息，按照下式

在实际刀轨上计算出各采样时刻对应的刀具位置点；

式中，t表示采样时刻； 

p_en表示直线刀轨段的起点；p_st表示直线刀轨段的终点；o_CTP表示圆弧刀轨段的圆心；R_CTP表示圆弧刀轨段的半径；

(3)根据步骤(1)给定的工件表面矢量方程，通过偏置，计算出理论刀轨的矢量方程p_e(u)＝[X_t(u)，Y_t(u)，Z_t(u)]^T，

式中，

Z_t(u)等于实际刀轨信息中刀心点的Z坐标值；

(4)按照下式，将步骤(2)中得到的刀具位置点p_a(t)投影到理论刀轨上，得到等效刀具位置点p_e(u(t))，

(p_e(u(t))-p_a(t))×n_a(t)＝0

解非线性方程组，得到参数u(t)，将参数u(t)带入理论刀轨矢量方程p_e(u)中，即得到等效刀具 位置点p_e(u(t))；

(5)将步骤(4)中解得的参数u(t)代入下式，即得到等效进给方向、等效法向量以及等效曲率，以等效进给方向、等效法向量以及等效曲率作为刀具位置点处的实际进给方向、实际法向量以及实际曲率；

f_e(ut))＝[X′_t(t)  Y′_t(t)  0]^T

n_e(u(t))＝[0 0 1]^T×f_e(u(t))

同时计算出等效进给方向在整体坐标系中的角位置；

式中，i＝[1 0 0]^T；

(6)将刀具参与切削的区域沿轴向划分为N个等高梁段，根据步骤(1)给定的工件毛坯边界矢量方程、立铣刀几何参数和步骤(2)中得到的刀具位置点，通过下式计算刀刃片{i，j}的切入切出角和切出角；

(a)切入角：

式中，p_W(v_i，j，en(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点，

N_f是刀刃数，j＝1，2，…，N，z是刀刃片{i，j}中点的Z向高度；

(b)切出角：

刀具切入工件阶段；

刀具在其他切削阶段；

式中， 

p_D(t，m)为满足方程组 的点，

m表示当前刀齿之前的第m个刀齿，m＝1，...，N_f

p_W(v_{i，j，B，ex}(t))为满足方程|p_a(t)-p_W(v)|＝RAD_i，j ²的点；

(7)通过下式计算作用在刀刃片{i，j}的铣削力：

式中，K_i，j，T(t)，K_i，j，R(t)分别是与h_i，j(t)相关的切向合径向瞬时铣削力系数，

Δa是等高梁段的高度，

h_i，j(t，m)是利用步骤(5)中得到的等效曲率K_e(u(t))，通过下式计算

凸型曲面；

凹型曲面；

式中，

(8)将各个刀刃上的力转化到X_S，Y_S和Z_S方向：

式中， 是刀具在t时刻处与刀刃片{i，j}对应的切削角度，被定义为从Y_S向顺时针到刀刃片{i，j}的中点所转过的角度；

(9)对于每个侧刃，将作用在所有刀刃片上的微元力求和，求得t时刻局部坐标系下作用于各个侧刃的铣削合力：

(10)将局部坐标系下的铣削合力转化到X，Y和Z方向：

Images