CN105975766A - 一种微铣削热力耦合解析建模方法 - Google Patents

Abstract

本发明一种微铣削热力耦合解析建模方法属于微铣削切削领域，涉及一种考虑刃口圆弧半径的微铣削力和温度预测模型，以及通过解析法对切削力和温度进行耦合计算的方法。建模方法基于微铣削瞬时切削厚度模型，以最小切削厚度为分界点分别建立以剪切效应和以耕犁效应为主导的微切削力解析模型。在耕犁效应为主导的微切削力则是基于耕犁力和过盈体积之间的关系建立；将微铣削切削区域简化为移动的有限长热源，建立切削温度模型；经过热耦合计算，实现对微铣削力和温度的预测。建模方法将微铣削温度的作用加入到微铣削力模型的建立过程中，通过解析法建立微铣削力和温度预测模型，并通过热力耦合计算，实现微铣削力和温度的准确快速预测。

Description

一种微铣削热力耦合解析建模方法

技术领域

本发明属于微铣削切削领域，涉及一种考虑刃口圆弧半径的微铣削力和温度预测模型，以及通过解析法对切削力和温度进行耦合计算的方法。

背景技术

基于微小型机床的微铣削技术是加工微小零件和高精密零件的一种新兴加工技术，具有加工材料范围广、能实现三维曲面铣削、加工精度高、能耗小、设备投资少、效率高等突出优点。由于加工尺度的急剧减小，微铣削加工过程中存在最小切削厚度现象，当实际切削厚度大于最小切削厚度时工件材料发生剪切滑移变形，形成连续切屑；当实际切削厚度小于最小切削厚度时，没有切屑产生，仅对工件材料进行耕犁和摩擦，导致切削厚度积累，当积累的瞬时切削厚度大于最小切削厚度时，产生切屑，由此导致切削力动态变化。微铣削加工过程中，刀具和工件的相对运动引起工件微观组织各部分的应力分布变化和接触界面摩擦热变化，局部高应力引发塑性应变而产生热并与摩擦产生的热共同构成并影响切削区温度场，温度场又反过来通过材料本构关系影响材料的力学性能，材料力学性能的改变又导致切削区应力变化。因此，在微铣削加工过程中，切削力和切削热相互影响最终达到动态平衡，是一个典型的热力耦合作用过程。

目前已经存在不少关于微铣削力的解析预测模型，但这些微铣削力模型几乎没有或很少将微铣削温度的作用加入到微铣削力模型的建立过程中，而微铣削加工是一个切削热与切削力耦合作用的过程。因此力模型中忽略温度的影响，势必影响力模型的预测精度。而关于微铣削热力耦合作用的研究大多基于有限元分析，有限元方法在用到微铣削领域时，由于需要划分较密的网格，因此具有耗时长的特点。还有些研究通过实验来观察切削力和温度的变化规律，难以反映造成这种变化规律的内在机制。因此有必要通过解析法建立微铣削力和温度预测模型，并通过热力耦合计算，实现微铣削力和温度的准确快速预测，进而揭示微铣削的成形机理，为后续表面残余应力和加工硬化的研究奠定理论基础。

来新民等2008年在International Journal of Machine Tools&Manufacture期刊第48卷第1期第1～14页发表的论文“Modeling and analysis of micro scale millingconsidering size effect,micro cutter edge radius and minimum chip thickness”中对微铣削力进行了理论建模，在微细切削有限元仿真分析的基础上，考虑切削刃刃口圆弧半径及积削瘤对切削力的影响，应用滑移线场分析法，以最小切削厚度为分界点，分别建立以剪切效应为主和以犁切效应为主的微细切削力预测模型，此模型同时考虑了切削厚度累计对切削力的影响。但该模型中没有考虑切削温度的影响，也无法对切削温度进行预测，而且该模型属于半经验模型，无法从理论上揭示微铣削的成形原理。而本发明中所阐述的方法则有效的弥补了这些缺点，更好的揭示微铣削成形原理，并准确预测微铣削力和温度。

发明内容

本发明所要解决的技术问题是克服现有技术的不足，考虑微铣削中存在的尺度效应，基于微铣削瞬时切削厚度模型，以最小切削厚度为分界点分别建立以剪切效应和以耕犁效应为主导的微切削力解析模型；并考虑切削温度的影响，建立了切削温度模型；将所建立的微铣削力模型和切削温度模型进行耦合计算，得到了考虑切削温度的微铣削力模型，更加符合微切削的实际情况，实现了微铣削力和温度的准确快速预测。

本发明采用的技术方案是一种微铣削热力耦合解析建模方法，其特征是，建模方法基于微铣削瞬时切削厚度模型，以最小切削厚度为分界点分别建立以剪切效应和以耕犁效应为主导的微切削力解析模型，当切削厚度大于最小切削厚度时，利用以剪切效应为主导的切削力模型进行预测；当切削厚度小于最小切削厚度时，利用以耕犁效应为主导的微铣削力模型进行预测；在耕犁效应为主导的微切削力则是基于耕犁力和过盈体积之间的关系建立；将微铣削切削区域简化为移动的有限长热源，建立切削温度模型；经过热耦合计算，实现对微铣削力和温度的预测；方法的具体步骤如下：

步骤1、计算微铣削瞬时切削厚度

首先，根据公式(1)判断是否发生单齿切削效应；

其中，f为进给率，单位为mm/s；n_s为主轴转速，单位为r/min；K_t为铣刀总齿数；R_t为刀齿齿尖径向跳动，单位为mm；为跳动初始角；f_t为每齿进给量，单位为mm/z；

当切削参数满足ξ(f,n_s,K_t,R_t,)＜0时，即出现单齿切削现象；当切削参数满足ξ(f,n_s,K_t,R_t,)≥0时，即为多齿交替切削；

然后计算与在t时刻刀具切削位置对应的上一次刀具切削到该位置时所对应的时刻t′；当出现单齿切削现象时，t′采用公式(2)计算；而当多齿交替切削时，t′采用公式(3)计算；

其中：R为刀具半径，单位为mm；t表示时间，单位为s；k为刀齿编号，k＝0,1,2,...,K_t-1；K_t为铣刀总齿数；ω为刀具转动的角速度，单位为rad/s；

采用Newton-Raphson迭代算法求解公式(2)和公式(3)，给定初始迭代值：

t′₁＝t-2π/(ωK)t′_i+1＝t′_i-F(t′_i)/F′(t′_i) (4)

其中F′(t′)为公式(2)和公式(3)的导函数；

求解出t′之后，即可计算实际切削厚度；t时刻切削厚度t_c可表示为公式(5)；

$t_c=\left\{\begin{array}{cc}{t}_c(t-2\mathrm{\&pi;}/({\mathrm{\&omega;K}}_t),k-1)+t_c(t,k)& t_c(t-2\mathrm{\&pi;}/({\mathrm{\&omega;K}}_t),k-1)<t_{\min }\\ t_c(t,k)& t_c(t-2\mathrm{\&pi;}/({\mathrm{\&omega;K}}_t),k-1)\&GreaterEqual;t_{\min }\end{array}\right.---\left(5\right)$

其中，t_min为最小切削厚度，单位为mm；t_c(t,k)表示第k齿在t时刻的切削厚度，公式(1)中的t_c(t,k)可根据公式(6)计算：

$t_c(t,k)=R+f_c\sin (\mathrm{\&omega;}t-2k\mathrm{\&pi;}/K_t+{\mathrm{\&omega;}}_0)-\sqrt{R^2-f_c^2{\cos }^2(\mathrm{\&omega;}t-2k\mathrm{\&pi;}/K_t+{\mathrm{\&omega;}}_0)}---\left(6\right)$

其中：R为微径铣刀半径，单位为mm；

步骤2、比较实际切削厚度与最小切削厚度的大小；

当实际切削厚度大于最小切削厚度时，建立以剪切效应为主导的切削过程中，使材料变形形成切屑的剪切力shear_F_s的计算模型；

在形成切屑的剪切力计算过程中，将刀具看作绝对锋利的刀具，单独考虑剪切作用；假设剪切区主剪切面上的剪应力是均匀分布的；以剪切效应为主导的切削过程中，剪切力的切向分力shear_F_sc、径向分力shear_F_sr和轴向分力shear_F_sa可分别表示为公式(8)；

$\begin{array}{c}shear\_F_{sc}=\cos B_h\&CenterDot;\frac{t_c{\mathrm{w\&tau;}}_{sm}\cos ({\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e)}{\sin \mathrm{\&phi;}\cos (\mathrm{\&phi;}+{\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e)}\\ shear\_F_{sr}=\frac{t_c{\mathrm{w\&tau;}}_{sm}\sin ({\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e)}{\sin \mathrm{\&phi;}\cos (\mathrm{\&phi;}+{\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e)}\\ shear\_F_{sa}=\sin B_h\&CenterDot;\frac{t_c{\mathrm{w\&tau;}}_{sm}\cos ({\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e)}{\sin \mathrm{\&phi;}\cos (\mathrm{\&phi;}+{\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e)}\end{array}---\left(8\right)$

其中，B_h为刀具螺旋角；t_c为切削厚度，单位为mm，在步骤1中已计算得到，w为切削宽度，单位为mm，τ_sm是主剪切面上的剪应力，单位为MPa，φ是剪切角，β_f是刀具与切屑的摩擦角，α_e是有效刀具前角；

计算切削宽度w，切削宽度为w＝a_p/cosB_h，a_p为切削深度，单位为mm；

计算剪切角φ，利用Merchant公式计算剪切角，如公式(9)；

$\mathrm{\&phi;}=\frac{\mathrm{\&pi;}}{4}-\frac{{\mathrm{\&beta;}}_f}{2}+\frac{{\mathrm{\&alpha;}}_e}{2}---\left(9\right)$

其中β_f为摩擦角，根据刀具与工件材料之间的摩擦特性获得；α_e为有效刀具前角，有效刀具前角由公式(10)获得；

${\mathrm{\&alpha;}}_e=\left\{\begin{array}{cc}\arcsin \left(\frac{t_c-r_e}{r_e}\right)& t_c\&le;r_e\\ {\mathrm{\&alpha;}}_0& t_c>r_e\end{array}\right.---\left(11\right)$

其中，t_c为切削厚度，单位为mm；r_e为刃口圆弧半径；α₀为刀具名义前角；

主剪切区的剪应力τ_s和主剪切面上的剪应力τ_sm的计算，主剪切区的剪应力τ_s根据Johnson-Cook本构模型进行计算，按公式(11)计算：

${\mathrm{\&tau;}}_s=\frac{1}{\sqrt{3}}\&lsqb;A+B\&CenterDot;{\mathrm{\&epsiv;}}^n\&rsqb;\&CenterDot;\&lsqb;1+C\ln \frac{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}{{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_0}\&rsqb;\&lsqb;1-{\left(\frac{T_{start}-T_r}{T_{melt}-T_r}\right)}^m\&rsqb;---\left(11\right)$

A，B,C,n，m为工件材料的Johnson-Cook本构模型参数；T_start为切削温度，T_r为参考温度，T_melt为工件材料的熔点；主剪切区剪应变率和剪应变ε由公式(12)和公式(13)获得；

$\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}=\left\{\begin{array}{cc}\frac{y{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m}{lh}& (y\&Element;\&lsqb;0,lh\&rsqb;)\\ \frac{h-y}{(1-l)h}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m& (y\&Element;(lh,h\&rsqb;)\end{array}\right.---\left(12\right)$

$\mathrm{\&epsiv;}\left\{\begin{array}{cc}\frac{{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m{y}^2}{2lhv\sin \mathrm{\&phi;}}& (y\&Element;\&lsqb;0,lh\&rsqb;)\\ \frac{{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m(-y^2+2hy-{\mathrm{lh}}^2)}{2(1-l)hv\sin \mathrm{\&phi;}}& (y\&Element;(lh,h\&rsqb;)\end{array}\right.---\left(13\right)$

其中：v为切削速度，单位为mm/s；

$\begin{array}{ccc}l=\frac{cos\mathrm{\&phi;}cos(\mathrm{\&phi;}-{\mathrm{\&alpha;}}_e)}{{\mathrm{cos\&alpha;}}_e}& {\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m=\frac{2v{\mathrm{cos\&alpha;}}_e}{hcos(\mathrm{\&phi;}-{\mathrm{\&alpha;}}_e)}& h=\frac{t_c}{10\sin \mathrm{\&phi;}}\end{array}---\left(14\right)$

当y＝lh时，带入公式(13)即为主剪切面上应变ε_m；将和ε_m及其余参数带入公式(11)计算得到τ_sm；

将计算得到的切削厚度t_c、切削宽度w、剪切角φ、摩擦角β_f、刀具有效前角α_e以及主剪切面上的剪应力τ_sm带入公式(1)计算得到以剪切效应为主导的微铣削力中形成切屑的剪切力shear_F_s；

步骤3、建立以剪切效应为主导的切削过程中，刃口圆弧对工件材料的耕犁力shear_F_p的计算模型；

应用Waldorf滑移线场理论计算当切削厚度大于最小切削厚度时，刃口圆弧对工件材料的耕犁力的切向分力shear_F_pc、径向分力shear_F_pr和轴向分力shear_F_pa分别表示为公式(15)；

shear_F_pc＝cos B_h·μ·τ_sm·w·[(1+2θ_w+2γ_w+sin(2η))sin(φ-γ_w+η)+cos(2η)cos(φ-γ_w+η)]·AC

shear_F_pr＝μ·τ_sm·w·[(1+2θ_w+2γ_w+sin(2η))cos(φ-γ_w+η)-cos(2η)sin(φ-γ_w+η)]·AC (15)

shear_F_pa＝sinB_h·μ·τ_sm·w·[(1+2θ+2γ_w+sin(2η))sin(φ-γ_w+η)+cos(2η)cos(φ-γ_w+η)]·AC

其中:滑移线和积屑区底部AC之间的角度：η＝0.5arccosμ

以A为顶点的扇形圆心角：

以B为顶点的扇形圆心角：

以A为顶点的扇形和以B为顶点的扇形的半径：

$R_w=\sin \mathrm{\&eta;}\sqrt{{\&lsqb;r_e\tan (\frac{\mathrm{\&pi;}}{4}+\frac{{\mathrm{\&alpha;}}_e}{2})+\frac{\sqrt{2}{R}_w{\mathrm{sin\&rho;}}_0}{\tan (\frac{\mathrm{\&pi;}}{4}+\frac{{\mathrm{\&alpha;}}_e}{2})}\&rsqb;}^2+2{\left(R_w{\mathrm{sin\&rho;}}_0\right)}^2}$

积屑区的底线长度AC则为：

上式中，μ为刀具与切屑之间的摩擦因子；ρ₀为切屑与为加工表面之间过渡斜面倾斜角；φ为剪切角；

步骤4、叠加步骤2和步骤3中计算得到的剪切力shear_F_s和耕犁力shear_F_p，获得以剪切效应为主导的微铣削力shear_F，其切向分力shear_F_c、径向分力shear_F_r和轴向分力shear_F_a按公式(16)计算：

shear_F_c＝shear_F_sc+shear_F_pc

shear_F_r＝shear_F_sr+shear_F_pr (16)

shear_F_a＝shear_F_sa+shear_F_pa

步骤5、如果步骤2的比较中，实际切削厚度小于最小切削厚度，则建立以耕犁效应为主导的微铣削力预测模型；假设当切削厚度小于最小切削厚度时，刀具对工件的耕犁力大小与耕犁区域过盈体积成比例，则耕犁力的切向分力plow_F_c、径向分力plow_F_r和轴向分力plow_F_a由公式(17)获得；

$\left\{\begin{array}{c}plow\_F_c=K_{cpp}\&CenterDot;A_p\&CenterDot;w\\ plow\_F_r=K_{rpp}\&CenterDot;A_p\&CenterDot;w\\ plow\_F_a=K_{app}\&CenterDot;A_p\&CenterDot;w\end{array}\right.---\left(17\right)$

其中，w为切削宽度，单位为mm；K_cpp，K_rpp和K_rpp分别为切向、径向和轴向耕犁效应力系数，单位为：N/mm³；A_p为耕犁区域过盈面积，单位为：mm²；

耕犁区域过盈面积可根据公式(18)计算；

$\left\{\begin{array}{cc}{A}_p=\frac{1}{2}{r}_e^2({\mathrm{\&alpha;}}_c+{\mathrm{\&alpha;}}_e)+\frac{1}{2}{r}_e\&CenterDot;l_{AB}-\frac{1}{2}{r}_e\&CenterDot;l_{BO}\&CenterDot;\sin ({\mathrm{\&alpha;}}_c+{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_p)& (\mathrm{\&delta;}<t_c<t_{min})\\ A_p=\frac{1}{2}{r}_e^2({\mathrm{\&alpha;}}_c+{\mathrm{\&alpha;}}_e)+\frac{1}{2}{r}_e\&CenterDot;l_{AE}-\frac{1}{2}{r}_e\&CenterDot;l_{EO}\&CenterDot;\sin ({\mathrm{\&alpha;}}_D+{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_{pe})& (\mathrm{\&delta;}<t_c<t_{\min })\end{array}\right.---\left(18\right)$

其中，r_e是刃口圆弧半径，单位为mm；t_c为切削厚度，单位为mm；δ为弹性回复量，单位为mm；式中其他参数如下：

$\begin{array}{cccc}{\mathrm{\&alpha;}}_C={\cos }^{-1}\left(\frac{r_e-t_c}{r_e}\right)& l_{AB}=\frac{\mathrm{\&delta;}-r_e(1-{\mathrm{cos\&alpha;}}_e)}{\sin \mathrm{\&alpha;}}& l_{BO}=\sqrt{{r_e}^2+{l_{AB}}^2}& {\mathrm{\&alpha;}}_P={\tan }^{-1}\left(\frac{l_{AB}}{r_e}\right)\\ {\mathrm{\&alpha;}}_D={\cos }^{-1}\left(\frac{r_e-t_c}{r_e}\right)& l_{AE}=\frac{t_c-r_e(1-{\mathrm{cos\&alpha;}}_e)}{\sin \mathrm{\&alpha;}}& l_{EO}=\sqrt{{r_e}^2+{l_{AE}}^2}& {\mathrm{\&alpha;}}_{Pe}={\tan }^{-1}\left(\frac{l_{AE}}{r_e}\right)\end{array}$

根据公式(19)计算得到K_cpp，K_rpp和K_rpp；

$\left\{\begin{array}{c}{K}_{cpp}\&CenterDot;A_p^{\min \_thickness}\&CenterDot;w=shear\_F_c\left(t_c^{\min \_thickness}\right)\\ K_{rpp}\&CenterDot;A_p^{\min \_thickness}\&CenterDot;w=shear\_F_r\left(t_c^{\min \_thickness}\right)\\ K_{app}\&CenterDot;A_p^{\min \_thickness}\&CenterDot;w=shear\_F_a\left(t_c^{\min \_thickness}\right)\end{array}\right.---\left(19\right)$

其中为当切削厚度为最小切削厚度时耕犁区域的过盈面积，单位为mm²；w为切削宽度，单位为mm； 和分别为当切削厚度为最小切削厚度，通过以剪切效应为主导的微铣削力模型计算得到的微铣削力；

步骤6、将计算得到的切向力F_c、径向力F_r和轴向力F_a分解到分别与X轴、Y轴、Z轴平行的力，得到X方向F_x，Y方向F_y，Z方向F_z，按公式(20)计算：

$\begin{array}{c}{F}_x=a\&CenterDot;(F_c\&CenterDot;\cos \mathrm{\&theta;}+F_r\&CenterDot;\sin \mathrm{\&theta;})\\ F_y=a\&CenterDot;(F_c\&CenterDot;\sin \mathrm{\&theta;}-F_r\&CenterDot;\cos \mathrm{\&theta;})\\ F_z=b\&CenterDot;F_a\end{array}---\left(20\right)$

其中：

θ为齿位角；

公式(20)中a、b为经验系数，由试验确定；

步骤7、建立切削温度预测模型；根据傅里叶定律所述，在导热现象中，单位时间内通过给定截面的热量，正比于垂直于该界面方向上的温度变化率和截面面积，而热量传递的方向则与温度升高的方向相反，按公式(21)计算：

$\frac{dq}{dt}=-KS\frac{\&part;T}{\&part;\stackrel{\&RightArrow;}{n}}---\left(21\right)$

式中：为导热率，单位为W；K为热导率，单位为W/mK；S为空间中垂直于向量的导热面积，单位为m²；为在向量方向的温度梯度，单位为K/m；

假如与x,y,z三个坐标轴的夹角分别为α,β,γ，热源在坐标原点时，温度梯度表示为公式(22)；

$\frac{\&part;T}{\&part;\stackrel{\&RightArrow;}{n}}=\frac{\&part;T}{\&part;x}cos\mathrm{\&alpha;}+\frac{\&part;T}{\&part;y}cos\mathrm{\&beta;}+\frac{\&part;T}{\&part;z}cos\mathrm{\&gamma;}---\left(22\right)$

单位时间通过面积S热量为：

$\frac{dq}{dt}=\underset{s}{\&Integral;\&Integral;}K(\frac{\&part;T}{\&part;x}cos\mathrm{\&alpha;}+\frac{\&part;T}{\&part;y}cos\mathrm{\&beta;}+\frac{\&part;T}{\&part;z}cos\mathrm{\&gamma;})dS=\underset{v}{\&Integral;\&Integral;\&Integral;}K\&lsqb;\frac{{\&part;}^{2}T}{\&part;x^2}+\frac{{\&part;}^{2}T}{\&part;y^2}+\frac{{\&part;}^{2}T}{\&part;z^2}\&rsqb;dv---\left(23\right)$

根据热力学中的理论，单位时间通过面积S的热量还可以表示为：

$\frac{dq}{dt}=\frac{\&part;}{\&part;t}\underset{v}{\&Integral;\&Integral;\&Integral;}c\mathrm{\&rho;}T(x,y,z,t)dv=\underset{v}{\&Integral;\&Integral;\&Integral;}c\mathrm{\&rho;}\frac{\&part;T}{\&part;t}dv---\left(24\right)$

其中c为工件材料的比热容，单位为J/(Kg·K)；ρ为工件材料密度，单位为Kg/m³，t为时间变量，单位为s；

联立公式(23)和公式(24)得：

$\frac{\&part;T}{\&part;t}=\frac{K}{c\mathrm{\&rho;}}(\frac{{\&part;}^{2}T}{\&part;x^2}+\frac{{\&part;}^{2}T}{\&part;y^2}+\frac{{\&part;}^{2}T}{\&part;z^2})=k^{\&prime;}(\frac{{\&part;}^{2}T}{\&part;x^2}+\frac{{\&part;}^{2}T}{\&part;y^2}+\frac{{\&part;}^{2}T}{\&part;z^2})---\left(25\right)$

其中 $k^{\&prime;}=\frac{K}{c\mathrm{\&rho;}};$

设工件中某一点M点的坐标为(x,y,z)，且对公式(25)进行傅里叶变换，并做一系列数学变换，如公式(26)；

$\begin{array}{c}\frac{\&part;F}{\&part;t}=-k^{\&prime;}(x^2+y^2+z^2)F\\ \frac{dF}{F}=-k^{\&prime;}{L_M}^{2}dt\\ \ln F=-k^{\&prime;}{L_M}^{2}t+C\\ F=C^{\&prime;}{e}^{-k^{\&prime;}{L_M}^{2}t}\end{array}---\left(26\right)$

对上述结果进行傅里叶逆变换：

$\mathrm{\&Delta;}T{(x,y,z,t)}_{po\mathrm{int}}=\frac{C^{\&prime;}}{{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{{L_M}^2}{4k^{\&prime;}t}\right)}---\left(27\right)$

其中ΔT(x,y,z,t)_point为点热源在t时刻造成的温升；根据能量守恒定理，得到Q为点热源瞬间发出的热量，单位为W；

当点热源在原点时，造成空间中某一点M(x,y,z)的温升，按公式(28)计算

$\begin{array}{c}\mathrm{\&Delta;}T{(x,y,z,t)}_{po\mathrm{int}}=\frac{Q}{c\mathrm{\&rho;}{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{{L_M}^2}{4k^{\&prime;}t}\right)}\\ =\frac{Q}{c\mathrm{\&rho;}{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{x^2+y^2+z^2}{4k^{\&prime;}t}\right)}\end{array}---\left(28\right)$

将其推广到有限长线热源，假设该线热源与Z轴平行，在XOY平面投影在坐标原点，热源长度为L；该有限长线热源造成空间中某一点的温升，由公式(29)计算得到；

$\begin{array}{c}\mathrm{\&Delta;}T{(x,y,z,t)}_{line}=\frac{Q}{c\mathrm{\&rho;}{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{x^2+y^2}{4k^{\&prime;}t}\right)}\underset{0}{\overset{L}{\&Integral;}}{e}^{-\left(\frac{{(z-z_i)}^2}{4k^{\&prime;}t}\right)}{\mathrm{dz}}_i\\ =\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}t}{e}^{-\left(\frac{x^2+y^2}{4k^{\&prime;}t}\right)}\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}t}}\right)-erf\left(\frac{z-L}{\sqrt{4k^{\&prime;}t}}\right)\&rsqb;\end{array}---\left(29\right)$

其中：

当线热源在XOY平面的投影不在原点时：

$\begin{array}{c}\mathrm{\&Delta;}T{(x,y,z,t)}_{line}=\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}t}{e}^{-\left(\frac{{(x-x_i)}^2+{(y-y_i)}^2}{4k^{\&prime;}t}\right)}\\ \&times;\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}t}}\right)-erf\left(\frac{z-L}{\sqrt{4k^{\&prime;}t}}\right)\&rsqb;\end{array}---\left(30\right)$

其中x_i和y_i分别为线热源在XOY平面投影的横纵坐标，单位为m；微铣削中，将刀齿的每一次切削看作是一个移动的有限长线热源，线热源长度用切削宽度w＝a_p/cosB_h来表示，单位为m，微铣削线热源在XOY平面投影的位置可表示为公式(31)；

$\left\{\begin{array}{c}{x}_i=f\&CenterDot;\frac{\mathrm{\&theta;}}{\mathrm{\&omega;}}+R\&CenterDot;\sin \left(\mathrm{\&theta;}\right)\\ y_i=R\&CenterDot;\cos \left(\mathrm{\&theta;}\right)\end{array}\right.---\left(31\right)$

式中θ为刀齿转过的角度，当刀齿与Y轴正方向重合时，θ＝0；R为刀具半径，单位为m；f为进给量，单位为m/s；ω为刀具转动的角速度，单位为rad/s；

因此根据公式(30)和公式(31)，单个齿一次切削对工件中某一点M(x,y,z)造成的温升ΔT由公式(32)得到；

$\begin{array}{c}\mathrm{\&Delta;}T=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}(\mathrm{\&pi;}-\mathrm{\&theta;})}{e}^{-\left(\frac{{(x-f\frac{\mathrm{\&theta;}}{\mathrm{\&omega;}}-R\sin \mathrm{\&theta;})}^2{(y-R\cos \mathrm{\&theta;})}^2}{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}\right)}\\ \&times;\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)-erf\left(\frac{z-a_p/\cos B}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)\&rsqb;d\mathrm{\&theta;}\end{array}---\left(32\right)$

由于工件中某一点的温升是多齿不断进给切削造成的温升叠加，因此工件中某一点M(x,y,z)的实际总温升total_ΔT由公式(33)得到；

$\begin{array}{c}total\_\mathrm{\&Delta;}T=\underset{i=0}{\overset{\frac{x-R}{f_t}}{\&Sigma;}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}(\mathrm{\&pi;}-\mathrm{\&theta;})}{e}^{-\left(\frac{{(x-i\&CenterDot;f_t-f\frac{\mathrm{\&theta;}}{\mathrm{\&omega;}}-R\sin \mathrm{\&theta;})}^2+{(y-R\cos \mathrm{\&theta;})}^2}{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}\right)}\\ \&times;\&lsqb;ref\left(\frac{z}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)-erf\left(\frac{z-a_p/\cos B_h}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)\&rsqb;d\mathrm{\&theta;}\end{array}---\left(33\right)$

式中；F_c为微铣削力的切向分力，单位为N；J为功热当量；f_t为每齿进给量，单位为m/z；其中v为切削速度，单位为m/s；t_c为切削厚度，单位为m；

γ_s＝cotφ+tan(φ-α)；

切削温度表示为公式(34)：

T_end＝T_room+total_ΔT (34)

其中：T_room为室温；

步骤8、对微铣削力模型和切削温度模型进行热力耦合计算，首先给定切削参数：主轴转速、每齿进给量和切削深度，刀具几何参数：刀具半径、刃口圆弧半径和螺旋角，工件材料的机械物理性能：弹性模量、泊松比以及热力学参数；然后给定温度初值T_start，取室温温度；然后通过微铣削温度模型计算切削温度T_end，比较T_start和T_end，若两者差值的绝对值大于ζ，则对T_start进行修正，重复步骤1到步骤7，直到当T_start和T_end差值小于ζ时，认为热力耦合计算达到平衡，输出此时的微铣削力和温度，实现考虑切削温度的微铣削力计算；

步骤9、由试验确定公式(20)中a、b为经验系数后，带入公式(20)，完成微铣削力和切削温度耦合计算模型。

本发明的有益效果是建模方法基于微铣削瞬时切削厚度模型，通过解析法建立微铣削力和温度预测模型，将微铣削温度的作用加入到微铣削力模型的建立过程中，并通过热力耦合计算，实现微铣削力和温度的准确快速预测。

附图说明

图1微铣削热力耦合模型建立过程图

图2微铣削热力耦合计算流程

图3当实际切削厚度小于等于刃口圆弧半径时有效刀具前角的计算，1-刀具，2-工件，α_e-有效刀具前角，t_c-实际切削厚度，r_e-刃口圆弧半径

图4当实际切削厚度大于刃口圆弧半径时有效刀具前角的计算，1-刀具，2-工件，α_e-有效刀具前角，t_c-实际切削厚度，r_e-刃口圆弧半径

图5 Waldorf滑移线场模型，1-积屑区，2-刀具，3-工件，α_e-有效刀具前角，r_e-刃口圆弧半径，t_c-实际切削厚度，φ-剪切角，η-滑移线和积屑区底部AC之间的角度，γ_w-以A为顶点的扇形圆心角，θ_w-以B为顶点的扇形圆心角，R_w-以A为顶点的扇形和以B为顶点的扇形的半径，ρ₀-切屑与为加工表面之间过渡斜面倾斜角。

图6微铣削热源移动示意图，a)微铣削热源A-A剖视图，b)微铣削热源俯视图，1-工件，2-线热源，3-刀具，F_c-微铣削力的切向分力，F_r-微铣削力的径向分力，F_x-微铣削力沿X轴方向的力，F_y-微铣削力沿Y轴方向的力，θ-刀齿的齿位角。

图7计算得到的切削厚度随切削时间的变化

具体实施方式

下面结合一个计算实例来详细说明该模型的计算过程。

镍基高温合金Inconel718具有良好的综合性能，如优良的高温强度、抗高温氧化、燃气腐蚀的能力以及抗疲劳性能等，可以很好的满足高温环境下具有高强度的微小结构/零件材料的要求。因此本发明以微铣削镍基高温合金Inconel718为例，对建模过程进行详细阐述。

实例所用条件如下：

主轴转速n为60000r/min，每齿进给量f_t为0.9μm/z，切削深度为30μm。铣刀总齿数K为2，刀齿齿尖径向跳动R_t为0.0117mm,跳动初始角为π/3，刀具直径为1mm，刃口圆弧半径为2μm，刀具名义前角为2°。工件材料为Inconel718，材料密度为8190kg/m³，热导率为11.2W/m·K，功热当量为4.1840，比热容为435J/(Kg·K)，材料熔点为1225℃。

建模方法基于微铣削瞬时切削厚度模型，以最小切削厚度为分界点分别建立以剪切效应和以耕犁效应为主导的微切削力解析模型，当切削厚度大于最小切削厚度时，利用以剪切效应为主导的切削力模型进行预测；当切削厚度小于最小切削厚度时，利用以耕犁效应为主导的微铣削力模型进行预测；在剪切效应为主导的微切削力建模过程中，将微铣削力分为两部分，一部分为使材料变形形成切屑的剪切力，另一部分为刃口圆弧对材料的耕犁力，对其分别进行建模；以耕犁效应为主导的微切削力则是基于耕犁力和过盈体积之间的关系建立；从傅里叶定律出发，推导出有限长热源对空间中某一点所造成的温升，将微铣削切削区域简化为移动的有限长热源，建立切削温度模型；最终，进行热力耦合计算，实现对微铣削力和温度的预测；

建模方法的具体步骤如下：

第1步：判断是否发生单齿切削效应。

根据公式(1)可计算得到ξ(f,n_s,K_t,R_t,)＝-0.0108＜0，所以发生了单齿切削效应。

因此利用公式(2)来计算在t时刻刀具切削位置对应的上一次刀具切削到该位置时所对应的时刻t'。然后再根据公式(5)、公式(6)、公式(7)。即可得到实际切削厚度随切削时间变化的图，如图7。该计算过程较为繁杂，因此利用Matlab编程实现。从图中可知，切削厚度随切削时间不断变化，切削厚度的最大值为：1.243μm。

第2步：比较实际切削厚度与最小切削厚度，本发明中最小切削厚度取0.7μm。

当实际切削厚度大于最小切削厚度时，利用步骤3到步骤5中建立的以剪切效应为主导的切削力模型进行预测；当实际切削厚度小于最小切削厚度时，利用步骤6中建立的以耕犁效应为主导的微铣削力模型进行预测。从图6可知，切削厚度在切削过程中从0到最大值不断变化，这里以计算当切削厚度为最大值时的切削力为例进行阐述。切削厚度的最大值为1.243μm，大于最小切削厚度0.7μm。因此采用以剪切效应为主导的切削力预测模型进行预测。

第3步：根据公式(8)计算当剪切效应主导时，刀具对工件的剪切力。公式(8)的参数：刀具螺旋角为B_h为30°，切削宽度w＝a_p/cos30°＝34.64μm＝0.03464mm，t_c在前文已计算得到，切削厚度最大值为1.243μm，即0.001243mm，剪切角φ可根据公式(9)和公式(10)计算得到，为45.22°，摩擦角β_f为21.8°，刀具有效前角α为22.24°，τ_sm根据公式(11)～公式(13)计算可得，其中温度T设置为T_start，初始化T_start为室温20℃。将以上参数带入公式(8)，即可得到剪切效应主导时，刀具对工件的剪切力在切向分力shear_F_sc，径向分力shear_F_sr，轴向分力shear_F_sa。

第4步：根据公式(15)计算以剪切效应为主导的切削过程中刃口圆弧对工件材料的耕犁力。公式(15)中μ取0.4，切屑与为加工表面之间过渡斜面倾斜角ρ₀为30°，其余参数与上文一致。通过Matlab变成计算得到耕犁力的切向分力shear_F_pc、径向分力shear_F_pr和轴向分力shear_F_pa。

第5步：根据公式(16)计算以剪切效应为主导的切削力，其切向分力为shear_F_c、径向分力为shear_F_c、轴向分力为shear_F_c。

第6步：根据公式(20)将上述计算得到的切向力F_c、径向力F_r和轴向力F_a分解到分别与X轴、Y轴、Z轴平行的力，得到X方向F_x，Y方向F_y，Z方向F_z。公式(20)中的经验系数通过试验获得，这里取为a＝0.6409,b＝1.3544。

第7步：计算切削温度。根据推导得到的公式(33)即可计算得到切削温度。其中材料密度为8190Kg/m³，功热当量为4.1840，材料热导率为11.2W/m·K，比热容为435J/(Kg·K)，其余参数与上文一致。即可计算得到切削温度T_end。

第8步：进行热力耦合计算。比较计算得到的T_start和T_end，若其差值的绝对值小于ζ(这里取0.5)。则输出此时计算得到的切削力和切削温度。若其差值的绝对值大于等于ζ，则修正T_start＝(T_start+T_end)/2，重新进行步骤3到步骤8，直到其差值小于ζ，则输出此时计算得到的切削力和切削温度。最终通过耦合运算得到微铣削力F_x为0.308N，F_y为0.3252N，F_z为0.2526N，切削温度为78.6℃。通过试验测得的铣削力F_x为0.304216N，F_y为0.32734N，F_z为0.239411N，切削温度为78.3℃。可以看到该模型具有良好的准确性。

本发明提出的微铣削热力耦合计算模型，可以在一定程度上反映在微铣削过程中，切削力和切削温度相互影响的过程，最终实现微铣削力和温度的准确快速预测，为实际生产中切削参数的优选、表面残余应力以及加工硬化的研究奠定理论基础。

Claims (1)

1.一种微铣削热力耦合解析建模方法，其特征是，建模方法基于微铣削瞬时切削厚度模型，以最小切削厚度为分界点分别建立以剪切效应和以耕犁效应为主导的微切削力解析模型，当切削厚度大于最小切削厚度时，利用以剪切效应为主导的切削力模型进行预测；当切削厚度小于最小切削厚度时，利用以耕犁效应为主导的微铣削力模型进行预测；在耕犁效应为主导的微切削力则是基于耕犁力和过盈体积之间的关系建立；将微铣削切削区域简化为移动的有限长热源，建立切削温度模型；经过热耦合计算，最终，实现对微铣削力和温度的预测；方法的具体步骤如下：

步骤1、计算微铣削瞬时切削厚度

首先，根据公式(1)判断是否发生单齿切削效应；

其中，f为进给率，单位为mm/s；n_s为主轴转速，单位为r/min；K_t为铣刀总齿数；R_t为刀齿齿尖径向跳动，单位为mm；为跳动初始角；f_t为每齿进给量，单位为mm/z；

当切削参数满足时，即出现单齿切削现象；当切削参数满足时，即为多齿交替切削；

然后计算与在t时刻刀具切削位置对应的上一次刀具切削到该位置时所对应的时刻t′；当出现单齿切削现象时，t′采用公式(2)计算；而当多齿交替切削时，t′采用公式(3)计算；

其中：R为刀具半径，单位为mm；t表示时间，单位为s；k为刀齿编号，k＝0,1,2,...,K_t-1；K_t为铣刀总齿数；ω为刀具转动的角速度，单位为rad/s；

采用Newton-Raphson迭代算法求解公式(2)和公式(3)，给定初始迭代值：

t′₁＝t-2π/(ωK) t′_i+1＝t′_i-F(t′_i)/F′(t′_i) (4)

其中F′(t′)为公式(2)和公式(3)的导函数；

求解出t′之后，即可计算实际切削厚度；t时刻切削厚度t_c可表示为公式(5)；

$t_c=\left\{\begin{array}{cc}{t}_c\left(t-2\mathrm{\&pi;}/\left({\mathrm{\&omega;K}}_t\right),k-1\right)+t_c\left(t,k\right)& t_c\left(t-2\mathrm{\&pi;}/\left({\mathrm{\&omega;K}}_t\right),k-1\right)<t_{\min }\\ t_c\left(t,k\right)& t_c\left(t-2\mathrm{\&pi;}/\left({\mathrm{\&omega;K}}_t\right),k-1\right)\&GreaterEqual;t_{\min }\end{array}\right.---\left(5\right)$

其中，t_min为最小切削厚度，单位为mm；t_c(t,k)表示第k齿在t时刻的切削厚度，公式(1)中的t_c(t,k)可根据公式(6)计算：

$t_c(t,k)=R+f_c\sin (\mathrm{\&omega;}t-2k\mathrm{\&pi;}/K_t+{\mathrm{\&omega;}}_0)-\sqrt{R^2-f_c^2{\cos }^2(\mathrm{\&omega;}t-2k\mathrm{\&pi;}/K_t+{\mathrm{\&omega;}}_0)}---\left(6\right)$

其中：R为微径铣刀半径，单位为mm；

步骤2、比较实际切削厚度与最小切削厚度的大小；

当实际切削厚度大于最小切削厚度时，建立以剪切效应为主导的切削过程中，使材料变形形成切屑的剪切力shear_F_s的计算模型；

在形成切屑的剪切力计算过程中，将刀具看作绝对锋利的刀具，单独考虑剪切作用；假设剪切区主剪切面上的剪应力是均匀分布的；以剪切效应为主导的切削过程中，剪切力的切向分力shear_F_sc、径向分力shear_F_sr和轴向分力shear_F_sa可分别表示为公式(8)；

$\begin{array}{c}shear\_F_{sc}=\cos B_h\&CenterDot;\frac{t_c{\mathrm{w\&tau;}}_{sm}\cos \left({\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e\right)}{\sin \mathrm{\&phi;}\cos \left(\mathrm{\&phi;}+{\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e\right)}\\ shear\_F_{sr}=\frac{t_c{\mathrm{w\&tau;}}_{sm}\sin \left({\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e\right)}{\sin \mathrm{\&phi;}\cos \left(\mathrm{\&phi;}+{\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e\right)}\\ shear\_F_{sa}=\sin B_h\&CenterDot;\frac{t_c{\mathrm{w\&tau;}}_{sm}\cos \left({\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e\right)}{\sin \mathrm{\&phi;}\cos \left(\mathrm{\&phi;}+{\mathrm{\&beta;}}_f-{\mathrm{\&alpha;}}_e\right)}\end{array}---\left(8\right)$

其中，B_h为刀具螺旋角；t_c为切削厚度，单位为mm，在步骤1中已计算得到，w为切削宽度，单位为mm，τ_sm是主剪切面上的剪应力，单位为MPa，φ是剪切角，β_f是刀具与切屑的摩擦角，α_e是有效刀具前角；

计算切削宽度w，切削宽度为w＝a_p/cosB_h，a_p为切削深度，单位为mm；

计算剪切角φ，利用Merchant公式计算剪切角，如公式(9)；

$\mathrm{\&phi;}=\frac{\mathrm{\&pi;}}{4}-\frac{{\mathrm{\&beta;}}_f}{2}+\frac{{\mathrm{\&alpha;}}_e}{2}---\left(9\right)$

其中β_f为摩擦角，根据刀具与工件材料之间的摩擦特性获得；α_e为有效刀具前角，有效刀具前角由公式(10)获得；

${\mathrm{\&alpha;}}_e=\left\{\begin{array}{cc}\arcsin \left(\frac{t_c-r_e}{r_e}\right)& t_c\&le;r_e\\ {\mathrm{\&alpha;}}_0& t_c>r_e\end{array}\right.---\left(10\right)$

其中，t_c为切削厚度，单位为mm；r_e为刃口圆弧半径；α₀为刀具名义前角；

主剪切区的剪应力τ_s和主剪切面上的剪应力τ_sm的计算，主剪切区的剪应力τ_s根据Johnson-Cook本构模型进行计算，按公式(11)计算：

${\mathrm{\&tau;}}_s=\frac{1}{\sqrt{3}}\&lsqb;A+B\&CenterDot;{\mathrm{\&epsiv;}}^n\&rsqb;\&CenterDot;\&lsqb;1+C\ln \frac{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}{{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_0}\&rsqb;\&lsqb;1-{\left(\frac{T_{start}-T_r}{T_{melt}-T_r}\right)}^m\&rsqb;---\left(11\right)$

A，B,C,n，m为工件材料的Johnson-Cook本构模型参数；T_start为切削温度，T_r为参考温度，T_melt为工件材料的熔点；主剪切区剪应变率和剪应变ε由公式(12)和公式(13)获得；

$\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}=\left\{\begin{array}{cc}\frac{y{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m}{lh}& (y\&Element;\&lsqb;0,lh\&rsqb;)\\ \frac{h-y}{(1-l)h}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m& (y\&Element;(lh,h\&rsqb;)\end{array}\right.---\left(12\right)$

$\mathrm{\&epsiv;}=\left\{\begin{array}{cc}\frac{{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m{y}^2}{2lhv\sin \mathrm{\&phi;}}& (y\&Element;\&lsqb;0,lh\&rsqb;)\\ \frac{{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m\left(-y^2+2hy-{\mathrm{lh}}^2\right)}{2\left(1-l\right)hv\sin \mathrm{\&phi;}}& (y\&Element;(lh,h\&rsqb;)\end{array}\right.---\left(13\right)$

其中：v为切削速度，单位为mm/s；

$\begin{array}{ccc}l=\frac{cos\mathrm{\&phi;}cos(\mathrm{\&phi;}-{\mathrm{\&alpha;}}_e)}{{\mathrm{cos\&alpha;}}_e}& {\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_m=\frac{2v{\mathrm{cos\&alpha;}}_e}{hcos(\mathrm{\&phi;}-{\mathrm{\&alpha;}}_e)}& h=\frac{t_c}{10\sin \mathrm{\&phi;}}\end{array}---\left(14\right)$

当y＝lh时，带入公式(13)即为主剪切面上应变ε_m；将和ε_m及其余参数带入公式(11)计算得到τ_sm；

将计算得到的切削厚度t_c、切削宽度w、剪切角φ、摩擦角β_f、刀具有效前角α_e以及主剪切面上的剪应力τ_sm带入公式(1)计算得到以剪切效应为主导的微铣削力中形成切屑的剪切力shear_F_s；

步骤3、建立以剪切效应为主导的切削过程中，刃口圆弧对工件材料的耕犁力shear_F_p的计算模型；

应用Waldorf滑移线场理论计算当切削厚度大于最小切削厚度时，刃口圆弧对工件材料的耕犁力的切向分力shear_F_pc、径向分力

shear_F_pr和轴向分力shear_F_pa分别表示为公式(15)；

shear_F_pc＝cos B_h·μ·τ_sm·w·[(1+2θ_w+2γ_w+sin(2η))sin(φ-γ_w+η)

+cos(2η)cos(φ-γ_w+η)]·AC

shear_F_pr＝μ·τ_sm·w·[(1+2θ_w+2γ_w+sin(2η))cos(φ-γ_w+η)

-cos(2η)sin(φ-γ_w+η)]·AC (15)

shear_F_pa＝sin B_h·μ·τ_sm·w·[(1+2θ+2γ_w+sin(2η))sin(φ-γ_w+η)

+cos(2η)cos(φ-γ_w+η)]·AC

其中:

滑移线和积屑区底部AC之间的角度：η＝0.5arccosμ

以A为顶点的扇形圆心角：

以B为顶点的扇形圆心角：

以A为顶点的扇形和以B为顶点的扇形的半径：

$R_w=\sin \mathrm{\&eta;}\sqrt{{\&lsqb;r_{e}tan(\frac{\mathrm{\&pi;}}{4}+\frac{{\mathrm{\&alpha;}}_e}{2})+\frac{\sqrt{2}{R}_w{\mathrm{sin\&rho;}}_0}{tan(\frac{\mathrm{\&pi;}}{4}+\frac{{\mathrm{\&alpha;}}_e}{2})}\&rsqb;}^2+2{\left(R_w{\mathrm{sin\&rho;}}_0\right)}^2}$

积屑区的底线长度AC则为：

上式中，μ为刀具与切屑之间的摩擦因子；ρ₀为切屑与为加工表面之间过渡斜面倾斜角；φ为剪切角；

步骤4、叠加步骤2和步骤3中计算得到的剪切力shear_F_s和耕犁力shear_F_p，获得以剪切效应为主导的微铣削力shear_F，其切向分力shear_F_c、径向分力shear_F_r和轴向分力shear_F_a按公式(16)计算：

shear_F_c＝shear_F_sc+shear_F_pc

shear_F_r＝shear_F_sr+shear_F_pr (16)

shear_F_a＝shear_F_sa+shear_F_pa

步骤5、如果步骤2的比较中，实际切削厚度小于最小切削厚度，则建立以耕犁效应为主导的微铣削力预测模型；假设当切削厚度小于最小切削厚度时，刀具对工件的耕犁力大小与耕犁区域过盈体积成比例，则耕犁力的切向分力plow_F_c、径向分力plow_F_r和轴向分力plow_F_a由公式(17)获得；

$\left\{\begin{array}{c}plow\_F_c=K_{cpp}\&CenterDot;A_p\&CenterDot;w\\ plow\_F_r=K_{rpp}\&CenterDot;A_p\&CenterDot;w\\ plow\_F_a=K_{app}\&CenterDot;A_p\&CenterDot;w\end{array}\right.---\left(17\right)$

其中，w为切削宽度，单位为mm；Kcpp，Krpp和Krpp分别为切向、径向和轴向耕犁效应力系数，单位为：N/mm³；A_p为耕犁区域过盈面积，单位为：mm²；

耕犁区域过盈面积可根据公式(18)计算；

$\left\{\begin{array}{cc}{A}_p=\frac{1}{2}{r}_e^2\left({\mathrm{\&alpha;}}_c+{\mathrm{\&alpha;}}_e\right)+\frac{1}{2}{r}_e\&CenterDot;l_{AB}-\frac{1}{2}{r}_e\&CenterDot;l_{BO}\&CenterDot;\sin \left({\mathrm{\&alpha;}}_c+{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_p\right)& \left(\mathrm{\&delta;}<t_c<t_{\min }\right)\\ A_p=\frac{1}{2}{r}_e^2\left({\mathrm{\&alpha;}}_c+{\mathrm{\&alpha;}}_e\right)+\frac{1}{2}{r}_e\&CenterDot;l_{AE}-\frac{1}{2}{r}_e\&CenterDot;l_{EO}\&CenterDot;\sin \left({\mathrm{\&alpha;}}_D+{\mathrm{\&alpha;}}_0+{\mathrm{\&alpha;}}_{pe}\right)& \left(\mathrm{\&delta;}<t_c<t_{\min }\right)\end{array}\right.---\left(18\right)$

其中，r_e是刃口圆弧半径，单位为mm；t_c为切削厚度，单位为mm；δ为弹性回复量，单位为mm；式中其他参数如下：

$\begin{array}{cccc}{\mathrm{\&alpha;}}_C={\cos }^{-1}\left(\frac{r_e-t_c}{r_e}\right)& l_{AB}=\frac{\mathrm{\&delta;}-r_e(1-{\mathrm{cos\&alpha;}}_e)}{\sin \mathrm{\&alpha;}}& l_{BO}=\sqrt{{r_e}^2+{l_{AB}}^2}& {\mathrm{\&alpha;}}_P={\tan }^{-1}\left(\frac{l_{AB}}{r_e}\right)\end{array}$

$\begin{array}{cccc}{\mathrm{\&alpha;}}_D={\cos }^{-1}\left(\frac{r_e-t_c}{r_e}\right)& l_{AE}=\frac{t_c-r_e(1-{\mathrm{cos\&alpha;}}_e)}{\sin \mathrm{\&alpha;}}& l_{EO}=\sqrt{{r_e}^2+{l_{AE}}^2}& {\mathrm{\&alpha;}}_{Pe}={\tan }^{-1}\left(\frac{l_{AE}}{r_e}\right)\end{array}$

根据公式(19)计算得到K_cpp，K_rpp和Kr_pp；

$\left\{\begin{array}{c}{K}_{cpp}\&CenterDot;A_p^{\min \_thickness}\&CenterDot;w=shear\_F_c\left(t_c^{\min \_thickness}\right)\\ K_{rpp}\&CenterDot;A_p^{\min \_thickness}\&CenterDot;w=shear\_F_r\left(t_c^{\min \_thickness}\right)\\ K_{app}\&CenterDot;A_p^{\min \_thickness}\&CenterDot;w=shear\_F_a\left(t_c^{\min \_thickness}\right)\end{array}\right.---\left(19\right)$

其中为当切削厚度为最小切削厚度时耕犁区域的过盈面积，单位为mm²；w为切削宽度，单位为mm； 和分别为当切削厚度为最小切削厚度，通过以剪切效应为主导的微铣削力模型计算得到的微铣削力；

步骤6、将计算得到的切向力F_c、径向力F_r和轴向力F_a分解到分别与X轴、Y轴、Z轴平行的力，得到X方向F_x，Y方向F_y，Z方向F_z，按公式(20)计算：

F_x＝a·(F_c·cosθ+F_r·sinθ)

F_y＝a·(F_c·sinθ-F_r·cosθ) (20)

F_z＝b·F_a

其中：

θ为齿位角；

公式(20)中a、b为经验系数，由试验确定；

步骤7、建立切削温度预测模型；根据傅里叶定律所述，在导热现象中，单位时间内通过给定截面的热量，正比于垂直于该界面方向上的温度变化率和截面面积，而热量传递的方向则与温度升高的方向相反，按公式(21)计算：

$\frac{dq}{dt}=-KS\frac{\&part;T}{\&part;\stackrel{\&RightArrow;}{n}}---\left(21\right)$

式中：为导热率，单位为W；K为热导率，单位为W/mK；S为空间中垂直于向量的导热面积，单位为m²；为在向量方向的温度梯度，单位为K/m；

假如与x,y,z三个坐标轴的夹角分别为α,β,γ，热源在坐标原点时，温度梯度表示为公式(22)；

$\frac{\&part;T}{\&part;\stackrel{\&RightArrow;}{n}}=\frac{\&part;T}{\&part;x}cos\mathrm{\&alpha;}+\frac{\&part;T}{\&part;y}cos\mathrm{\&beta;}+\frac{\&part;T}{\&part;z}cos\mathrm{\&gamma;}---\left(22\right)$

单位时间通过面积S热量为：

$\frac{dq}{dt}=\underset{s}{\&Integral;\&Integral;}K\left(\frac{\&part;T}{\&part;x}\cos \mathrm{\&alpha;}+\frac{\&part;T}{\&part;y}\cos \mathrm{\&beta;}+\frac{\&part;T}{\&part;z}\cos \mathrm{\&gamma;}\right)dS=\underset{v}{\&Integral;\&Integral;\&Integral;}K\&lsqb;\frac{{\&part;}^{2}T}{\&part;x^2}+\frac{{\&part;}^{2}T}{\&part;y^2}+\frac{{\&part;}^{2}T}{\&part;z^2}\&rsqb;dv---\left(23\right)$

根据热力学中的理论，单位时间通过面积S的热量还可以表示为：

$\frac{dq}{dt}=\frac{\&part;}{\&part;t}\underset{v}{\&Integral;\&Integral;\&Integral;}c\mathrm{\&rho;}T(x,y,z,t)dv=\underset{v}{\&Integral;\&Integral;\&Integral;}v\mathrm{\&rho;}\frac{\&part;T}{\&part;t}dv---\left(24\right)$

其中c为工件材料的比热容，单位为J/(Kg·K)；ρ为工件材料密度，单位为Kg/m³，t为时间变量，单位为s；

联立公式(23)和公式(24)得：

$\frac{\&part;T}{\&part;t}=\frac{K}{c\mathrm{\&rho;}}(\frac{{\&part;}^{2}T}{\&part;x^2}+\frac{{\&part;}^{2}T}{\&part;y^2}+\frac{{\&part;}^{2}T}{\&part;z^2})=k^{\&prime;}(\frac{{\&part;}^{2}T}{\&part;x^2}+\frac{{\&part;}^{2}T}{\&part;y^2}+\frac{{\&part;}^{2}T}{\&part;z^2})---\left(25\right)$

其中

设工件中某一点M点的坐标为(x,y,z)，且

对公式(25)进行傅里叶变换，并做一系列数学变换，如公式(26)；

$\begin{array}{c}\frac{\&part;F}{\&part;t}=-k^{\&prime;}\left(x^2+y^2+z^2\right)F\\ \frac{dF}{F}=-k^{\&prime;}{L_M}^{2}dt\\ \ln F=-k^{\&prime;}{L_M}^{2}t+C\\ F=C^{\&prime;}{e}^{-k^{\&prime;}{L_M}^{2}t}\end{array}---\left(26\right)$

对上述结果进行傅里叶逆变换：

$\mathrm{\&Delta;}T{(x,y,z,t)}_{po\mathrm{int}}=\frac{C^{\&prime;}}{{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{{L_M}^2}{4k^{\&prime;}t}\right)}---\left(27\right)$

其中ΔT(x,y,z,t)_point为点热源在t时刻造成的温升；根据能量守恒定理，得到Q为点热源瞬间发出的热量，单位为W；

当点热源在原点时，造成空间中某一点M(x,y,z)的温升，按公式(28)计算

$\begin{array}{c}\mathrm{\&Delta;}T{(x,y,z,t)}_{po\mathrm{int}}=\frac{Q}{c\mathrm{\&rho;}{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{{L_M}^2}{4k^{\&prime;}t}\right)}\\ =\frac{Q}{c\mathrm{\&rho;}{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{x^2+y^2+z^2}{4k^{\&prime;}t}\right)}\end{array}---\left(28\right)$

将其推广到有限长线热源，假设该线热源与Z轴平行，在XOY平面投影在坐标原点，热源长度为L；该有限长线热源造成空间中某一点的温升，由公式(29)计算得到；

$\begin{array}{c}\mathrm{\&Delta;}T{\left(x,y,z,t\right)}_{line}=\frac{Q}{c\mathrm{\&rho;}{\left(4{\mathrm{\&pi;k}}^{\&prime;}t\right)}^{3/2}}{e}^{-\left(\frac{x^2+y^2}{4k^{\&prime;}t}\right)}\underset{0}{\overset{L}{\&Integral;}}{e}^{-\left(\frac{{\left(z-z_i\right)}^2}{4k^{\&prime;}t}\right)}{\mathrm{dz}}_i\\ =\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}t}{e}^{-\left(\frac{x^2+y^2}{4k^{\&prime;}t}\right)}\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}t}}\right)-erf\left(\frac{z-L}{\sqrt{4k^{\&prime;}t}}\right)\&rsqb;\end{array}---\left(29\right)$

其中：

当线热源在XOY平面的投影不在原点时：

$\begin{array}{c}\mathrm{\&Delta;}T{\left(x,y,z,t\right)}_{line}=\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}t}{e}^{-\left(\frac{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2}{4k^{\&prime;}t}\right)}\\ \&times;\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}t}}\right)-erf\left(\frac{z-L}{\sqrt{4k^{\&prime;}t}}\right)\&rsqb;\end{array}---\left(30\right)$

其中x_i和y_i分别为线热源在XOY平面投影的横纵坐标，单位为m；微铣削中，将刀齿的每一次切削看作是一个移动的有限长线热源，线热源长度用切削宽度w＝a_p/cosB_h来表示，单位为m，微铣削线热源在XOY平面投影的位置可表示为公式(31)；

$\left\{\begin{array}{c}{x}_i=f\&CenterDot;\frac{\mathrm{\&theta;}}{\mathrm{\&omega;}}+R\&CenterDot;\sin \left(\mathrm{\&theta;}\right)\\ y_i=R\&CenterDot;\cos \left(\mathrm{\&theta;}\right)\end{array}\right.---\left(31\right)$

式中θ为刀齿转过的角度，当刀齿与Y轴正方向重合时，θ＝0；R为刀具半径，单位为m；f为进给量，单位为m/s；ω为刀具转动的角速度，单位为rad/s；

因此根据公式(30)和公式(31)，单个齿一次切削对工件中某一点M(x,y,z)造成的温升ΔT由公式(32)得到；

$\begin{array}{c}\mathrm{\&Delta;}T=\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}\left(\mathrm{\&pi;}-\mathrm{\&theta;}\right)}{e}^{-\left(\frac{{\left(x-f\frac{\mathrm{\&theta;}}{\mathrm{\&omega;}}-R\sin \mathrm{\&theta;}\right)}^2+{\left(y-R\cos \mathrm{\&theta;}\right)}^2}{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}\right)}\\ \&times;\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)-erf\left(\frac{z-a_p/\cos B}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)\&rsqb;d\mathrm{\&theta;}\end{array}---\left(32\right)$

由于工件中某一点的温升是多齿不断进给切削造成的温升叠加，因此工件中某一点M(x,y,z)的实际总温升total_ΔT由公式(33)得到；

$\begin{array}{c}total\_\mathrm{\&Delta;}T=\underset{i=0}{\overset{\frac{x-R}{f_t}}{\mathrm{\&Sigma;}}}\underset{0}{\overset{\mathrm{\&pi;}}{\&Integral;}}\frac{Q}{8{\mathrm{\&pi;c\&rho;k}}^{\&prime;}\left(\mathrm{\&pi;}-\mathrm{\&theta;}\right)}{e}^{-\left(\frac{{\left(x-i\&CenterDot;f_t-f\frac{\mathrm{\&theta;}}{\mathrm{\&omega;}}-R\sin \mathrm{\&theta;}\right)}^2+{\left(y-R\cos \mathrm{\&theta;}\right)}^2}{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}\right)}\\ \&times;\&lsqb;erf\left(\frac{z}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)-erf\left(\frac{z-a_p/\cos B_h}{\sqrt{4k^{\&prime;}\frac{\mathrm{\&pi;}-\mathrm{\&theta;}}{\mathrm{\&omega;}}}}\right)\&rsqb;d\mathrm{\&theta;}\end{array}---\left(33\right)$

式中F_c为微铣削力的切向分力，单位为N；J为功热当量；f_t为每齿进给量，单位为m/z；其中v为切削速度，单位为m/s；t_c为切削厚度，单位为m；

γ_s＝cotφ+tan(φ-α)；

切削温度表示为公式(34)：

T_end＝T_room+total_ΔT (34)

其中：T_room为室温；

步骤8、对微铣削力模型和切削温度模型进行热力耦合计算，首先给定切削参数：主轴转速、每齿进给量和切削深度，刀具几何参数：刀具半径、刃口圆弧半径和螺旋角，工件材料的机械物理性能：弹性模量、泊松比以及热力学参数；然后给定温度初值T_start，取室温温度；然后通过微铣削温度模型计算切削温度T_end，比较T_start和T_end，若两者差值的绝对值大于ζ，则对T_start进行修正，重复步骤1到步骤7，直到当T_start和T_end差值小于ζ时，认为热力耦合计算达到平衡，输出此时的微铣削力和温度，实现考虑切削温度的微铣削力计算；

步骤9、由试验确定公式(20)中a、b为经验系数后，带入公式(20)，完成微铣削力和切削温度耦合计算模型。