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

摘要

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

主权项

一种微铣削热力耦合解析建模方法，其特征是，建模方法基于微铣削瞬时切削厚度模型，以最小切削厚度为分界点分别建立以剪切效应和以耕犁效应为主导的微切削力解析模型，当切削厚度大于最小切削厚度时，利用以剪切效应为主导的切削力模型进行预测；当切削厚度小于最小切削厚度时，利用以耕犁效应为主导的微铣削力模型进行预测；在耕犁效应为主导的微切削力则是基于耕犁力和过盈体积之间的关系建立；将微铣削切削区域简化为移动的有限长热源，建立切削温度模型；经过热耦合计算，最终，实现对微铣削力和温度的预测；方法的具体步骤如下：步骤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\pi /\left({\mathrm{\omega K}}_t\right),k-1\right)+t_c\left(t,k\right)& t_c\left(t-2\pi /\left({\mathrm{\omega K}}_t\right),k-1\right)<t_{\min }\\ t_c\left(t,k\right)& t_c\left(t-2\pi /\left({\mathrm{\omega K}}_t\right),k-1\right)\ge 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 (\omega t-2k\pi /K_t+{\omega }_0)-\sqrt{R^2-f_c^2{\cos }^2(\omega t-2k\pi /K_t+{\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·\frac{t_c{\mathrm{w\tau }}_{sm}\cos \left({\beta }_f-{\alpha }_e\right)}{\sin \phi \cos \left(\phi +{\beta }_f-{\alpha }_e\right)}\\ shear\_F_{sr}=\frac{t_c{\mathrm{w\tau }}_{sm}\sin \left({\beta }_f-{\alpha }_e\right)}{\sin \phi \cos \left(\phi +{\beta }_f-{\alpha }_e\right)}\\ shear\_F_{sa}=\sin B_h·\frac{t_c{\mathrm{w\tau }}_{sm}\cos \left({\beta }_f-{\alpha }_e\right)}{\sin \phi \cos \left(\phi +{\beta }_f-{\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)；$\phi =\frac{\pi }{4}-\frac{{\beta }_f}{2}+\frac{{\alpha }_e}{2}---\left(9\right)$其中β_f为摩擦角，根据刀具与工件材料之间的摩擦特性获得；α_e为有效刀具前角，有效刀具前角由公式(10)获得；${\alpha }_e=\left\{\begin{array}{cc}\arcsin \left(\frac{t_c-r_e}{r_e}\right)& t_c\le r_e\\ {\alpha }_0& t_c>r_e\end{array}\right.---\left(10\right)$其中，t_c为切削厚度，单位为mm；r_e为刃口圆弧半径；α₀为刀具名义前角；主剪切区的剪应力τ_s和主剪切面上的剪应力τ_sm的计算，主剪切区的剪应力τ_s根据Johnson‑Cook本构模型进行计算，按公式(11)计算：${\tau }_s=\frac{1}{\sqrt{3}}[A+B·{ϵ}^n]·[1+C\ln \frac{\stackrel{·}{ϵ}}{{\stackrel{·}{ϵ}}_0}][1-{\left(\frac{T_{start}-T_r}{T_{melt}-T_r}\right)}^m]---\left(11\right)$A，B,C,n，m为工件材料的Johnson‑Cook本构模型参数；T_start为切削温度，T_r为参考温度，T_melt为工件材料的熔点；主剪切区剪应变率和剪应变ε由公式(12)和公式(13)获得；$\stackrel{·}{ϵ}=\left\{\begin{array}{cc}\frac{y{\stackrel{·}{ϵ}}_m}{lh}& (y\in [0,lh\left]\right)\\ \frac{h-y}{(1-l)h}{\stackrel{·}{ϵ}}_m& (y\in (lh,h\left]\right)\end{array}\right.---\left(12\right)$$ϵ=\left\{\begin{array}{cc}\frac{{\stackrel{·}{ϵ}}_m{y}^2}{2lhv\sin \phi }& (y\in [0,lh\left]\right)\\ \frac{{\stackrel{·}{ϵ}}_m\left(-y^2+2hy-{\mathrm{lh}}^2\right)}{2\left(1-l\right)hv\sin \phi }& (y\in (lh,h\left]\right)\end{array}\right.---\left(13\right)$其中：v为切削速度，单位为mm/s；$\begin{array}{ccc}l=\frac{cos\phi cos(\phi -{\alpha }_e)}{{\mathrm{cos\alpha }}_e}& {\stackrel{·}{ϵ}}_m=\frac{2v{\mathrm{cos\alpha }}_e}{hcos(\phi -{\alpha }_e)}& h=\frac{t_c}{10\sin \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+η)]·ACshear_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 \eta \sqrt{{[r_{e}tan(\frac{\pi }{4}+\frac{{\alpha }_e}{2})+\frac{\sqrt{2}{R}_w{\mathrm{sin\rho }}_0}{tan(\frac{\pi }{4}+\frac{{\alpha }_e}{2})}]}^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_pcshear_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}·A_p·w\\ plow\_F_r=K_{rpp}·A_p·w\\ plow\_F_a=K_{app}·A_p·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({\alpha }_c+{\alpha }_e\right)+\frac{1}{2}{r}_e·l_{AB}-\frac{1}{2}{r}_e·l_{BO}·\sin \left({\alpha }_c+{\alpha }_0+{\alpha }_p\right)& \left(\delta <t_c<t_{\min }\right)\\ A_p=\frac{1}{2}{r}_e^2\left({\alpha }_c+{\alpha }_e\right)+\frac{1}{2}{r}_e·l_{AE}-\frac{1}{2}{r}_e·l_{EO}·\sin \left({\alpha }_D+{\alpha }_0+{\alpha }_{pe}\right)& \left(\delta <t_c<t_{\min }\right)\end{array}\right.---\left(18\right)$其中，r_e是刃口圆弧半径，单位为mm；t_c为切削厚度，单位为mm；δ为弹性回复量，单位为mm；式中其他参数如下：$\begin{array}{cccc}{\alpha }_C={\cos }^{-1}\left(\frac{r_e-t_c}{r_e}\right)& l_{AB}=\frac{\delta -r_e(1-{\mathrm{cos\alpha }}_e)}{\sin \alpha }& l_{BO}=\sqrt{{r_e}^2+{l_{AB}}^2}& {\alpha }_P={\tan }^{-1}\left(\frac{l_{AB}}{r_e}\right)\end{array}$$\begin{array}{cccc}{\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 \alpha }& l_{EO}=\sqrt{{r_e}^2+{l_{AE}}^2}& {\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}·A_p^{\min \_thickness}·w=shear\_F_c\left(t_c^{\min \_thickness}\right)\\ K_{rpp}·A_p^{\min \_thickness}·w=shear\_F_r\left(t_c^{\min \_thickness}\right)\\ K_{app}·A_p^{\min \_thickness}·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{\partial T}{\partial \overrightarrow{n}}---\left(21\right)$式中：为导热率，单位为W；K为热导率，单位为W/mK；S为空间中垂直于向量的导热面积，单位为m²；为在向量方向的温度梯度，单位为K/m；假如与x,y,z三个坐标轴的夹角分别为α,β,γ，热源在坐标原点时，温度梯度表示为公式(22)；$\frac{\partial T}{\partial \overrightarrow{n}}=\frac{\partial T}{\partial x}cos\alpha +\frac{\partial T}{\partial y}cos\beta +\frac{\partial T}{\partial z}cos\gamma ---\left(22\right)$单位时间通过面积S热量为：$\frac{dq}{dt}=\underset{s}{\int \int }K\left(\frac{\partial T}{\partial x}\cos \alpha +\frac{\partial T}{\partial y}\cos \beta +\frac{\partial T}{\partial z}\cos \gamma \right)dS=\underset{v}{\int \int \int }K\left[\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right]dv---\left(23\right)$根据热力学中的理论，单位时间通过面积S的热量还可以表示为：$\frac{dq}{dt}=\frac{\partial }{\partial t}\underset{v}{\int \int \int }c\rho T(x,y,z,t)dv=\underset{v}{\int \int \int }v\rho \frac{\partial T}{\partial t}dv---\left(24\right)$其中c为工件材料的比热容，单位为J/(Kg·K)；ρ为工件材料密度，单位为Kg/m³，t为时间变量，单位为s；联立公式(23)和公式(24)得：$\frac{\partial T}{\partial t}=\frac{K}{c\rho }(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2})=k^{\prime }(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2})---\left(25\right)$其中设工件中某一点M点的坐标为(x,y,z)，且对公式(25)进行傅里叶变换，并做一系列数学变换，如公式(26)；$\begin{array}{c}\frac{\partial F}{\partial 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)$对上述结果进行傅里叶逆变换：$\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}\Delta T{(x,y,z,t)}_{po\mathrm{int}}=\frac{Q}{c\rho {\left(4{\mathrm{\pi k}}^{\prime }t\right)}^{3/2}}{e}^{-\left(\frac{{L_M}^2}{4k^{\prime }t}\right)}\\ =\frac{Q}{c\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}\Delta T{\left(x,y,z,t\right)}_{line}=\frac{Q}{c\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}{\int }}{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)}\left[erf\left(\frac{z}{\sqrt{4k^{\prime }t}}\right)-erf\left(\frac{z-L}{\sqrt{4k^{\prime }t}}\right)\right]\end{array}---\left(29\right)$其中：当线热源在XOY平面的投影不在原点时：$\begin{array}{c}\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 \left[erf\left(\frac{z}{\sqrt{4k^{\prime }t}}\right)-erf\left(\frac{z-L}{\sqrt{4k^{\prime }t}}\right)\right]\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·\frac{\theta }{\omega }+R·\sin \left(\theta \right)\\ y_i=R·\cos \left(\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}\Delta T=\underset{0}{\overset{\pi }{\int }}\frac{Q}{8{\mathrm{\pi c\rho k}}^{\prime }\left(\pi -\theta \right)}{e}^{-\left(\frac{{\left(x-f\frac{\theta }{\omega }-R\sin \theta \right)}^2+{\left(y-R\cos \theta \right)}^2}{4k^{\prime }\frac{\pi -\theta }{\omega }}\right)}\\ \times \left[erf\left(\frac{z}{\sqrt{4k^{\prime }\frac{\pi -\theta }{\omega }}}\right)-erf\left(\frac{z-a_p/\cos B}{\sqrt{4k^{\prime }\frac{\pi -\theta }{\omega }}}\right)\right]d\theta \end{array}---\left(32\right)$由于工件中某一点的温升是多齿不断进给切削造成的温升叠加，因此工件中某一点M(x,y,z)的实际总温升total_ΔT由公式(33)得到；$\begin{array}{c}total\_\Delta T=\underset{i=0}{\overset{\frac{x-R}{f_t}}{\Sigma }}\underset{0}{\overset{\pi }{\int }}\frac{Q}{8{\mathrm{\pi c\rho k}}^{\prime }\left(\pi -\theta \right)}{e}^{-\left(\frac{{\left(x-i·f_t-f\frac{\theta }{\omega }-R\sin \theta \right)}^2+{\left(y-R\cos \theta \right)}^2}{4k^{\prime }\frac{\pi -\theta }{\omega }}\right)}\\ \times \left[erf\left(\frac{z}{\sqrt{4k^{\prime }\frac{\pi -\theta }{\omega }}}\right)-erf\left(\frac{z-a_p/\cos B_h}{\sqrt{4k^{\prime }\frac{\pi -\theta }{\omega }}}\right)\right]d\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)，完成微铣削力和切削温度耦合计算模型。