黄贤振¹, 臧云飞¹, 宋增旺², 矫宝龙¹
1. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819;
2. 北方重工集团有限公司, 辽宁 沈阳 110141
收稿日期：2018-04-18
基金项目：国家自然科学基金资助项目(51575094)；国家博士后基金资助项目(2017M611244)；中央高校基本科研业务费专项资金资助项目(N160304004)。
作者简介：黄贤振(1982-), 男, 山东定陶人, 东北大学教授, 博士。

摘要：利用单频率法建立薄壁构件铣削稳定性及表面位置误差模型, 采用Kriging方法对薄壁构件铣削稳定性及表面位置误差进行了可靠性灵敏度分析, 评价铣削参数对薄壁构件铣削稳定性及表面位置误差的影响程度.研究结果表明:随薄壁构件铣削系统的阻尼比及y向刚度的增加, 系统的可靠度增加; 随切向力系数、固有频率、径向侵入比及轴向切深的增加, 系统可靠度降低; 径向力系数和每齿进给量对系统的影响较小.研究结果能为薄壁构件高速铣削加工提供合理的理论依据.
关键词：薄壁构件Kriging铣削稳定性表面位置误差可靠性灵敏度
Reliability Sensitivity Analysis of Thin-Walled Components Milling Machining
HUANG Xian-zhen¹, ZANG Yun-fei¹, SONG Zeng-wang², JIAO Bao-long¹
1. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China;
2. Northern Heavy Industries Group Co. Ltd., Shenyang 110141, China
Corresponding author: HUANG Xian-zhen, E-mail: xzhhuang@mail.neu.edu.cn
Abstract: A stability and surface position error model for the thin-walled components milling was established by using a single frequency method. The Kriging method was used to analyze the reliability sensitivity of the stability and surface position error of thin-walled components milling, and to evaluate the influence of milling parameters on the stability and surface position error of thin-walled components milling. The results show that the system reliability increases with the increase of damping ratio and y direction stiffness of the thin-walled components milling system. The system reliability decreases with the increase of tangential force coefficient, natural frequency, radial invasion ratio and axial cutting depth. The radial force factor and feed per tooth have less impact on the system. The research results may provide reasonable theoretical basis for high speed milling of thin-walled components.
Key words: thin-walled componentsKrigingmilling stabilitysurface position errorreliability sensitivity
薄壁构件刚性较差, 在加工过程中易发生振动和变形, 因此加工质量往往难以得到保证^[1].加工振动主要包括自激振动和强迫振动.自激振动根据系统内部反馈方式的不同分为摩擦颤振、再生颤振及振型耦合型颤振等^[2].强迫振动主要是由系统外部的周期激励引起的振动, 是产生表面位置误差的主要因素^[3].
Jemielniak等^[4]提出发生颤振的主要原因是动态切削厚度的再生效应.Li等^[5]基于时域法讨论了薄壁件铣削稳定性.Thevenot等^[6]考虑薄壁构件加工过程中刚度和固有频率的变化绘制三维稳定性图.Tang等^[7]绘制出基于最大去除率的三维稳定性图.Bravo等^[8]绘制出刀具和工件相对位置变化时的三维稳定性叶瓣图.Seguy^[9]绘制出薄壁件不同加工路径下的三维稳定性叶瓣图.上述文献中的铣削参数被认为是确定的, 但实际加工中往往都是随机不确定的, 尤其对于考虑薄壁构件表面位置误差与稳定性结合的研究不多.
1 铣削稳定性及表面位置误差模型1.1 铣削稳定性模型的建立常见的薄壁构件铣削加工动力学模型如图 1所示.其中Ω为主轴转速(r/min); j为第j刀齿; ?_j为第j刀齿的瞬时切削角(rad); F_tj和F_rj分别为第j刀齿受到的切向切削力(N)和径向切削力(N); f为工件相对刀具的进给方向; a_r为径向切深(mm); k_y和c_y分别为工件在y方向的等效刚度(N/m)和等效阻尼.本文将采用目前计算效率最高的单频率法对薄壁构件进行稳定性求解^[3].
图 1(Fig. 1)

图 1 薄壁构件单自由度铣削动力学模型Fig.1 Dynamic model of single-degree-of-freedom milling for thin-walled components

对薄壁构件的动力学模型分析后得

(1)

式中: y(t)为薄壁构件在y方向的位移; m_y为薄壁构件在y方向的等效质量; F_y(t)为薄壁构件在y方向的铣削力.
对式(1)进行拉氏变换后得

(2)

式中:s为拉普拉斯算子; 为铣削系统的固有频率; ξ_y=c_y/(2m_yω_N)为铣削系统的阻尼比.
由自动控制理论知, 当系统传递函数特征方程根s=σ+iω的实部σ=0时, 系统处于稳定和不稳定的临界状态.将s=iω代入式(2):

(3)

式中:ω为沿y方向的颤振频率; r=ω/ω_N.
通过力学模型可以得到铣削力^[6]:

(4)

式中:N为刀齿数; a_ax为铣削系统的轴向切深，m; α_yy为y方向系数; k_tc为铣削系统的切向力系数，N/m²; Δy为y向位移, m.
利用谐波函数, 将式(4)转化为

(5)

将欧拉公式e^-iωT=cosωT-isinωT代入式(5)并与式(3)进行对比后得到系统临界轴向切深计算公式:

(6)

(7)

将代入式(7), 可得到相应的主轴转速:

(8)

式中, k=0, 1, 2….
1.2 表面位置误差模型的建立在薄壁构件的铣削加工过程中, 产生表面位置误差的主要原因是强迫振动的发生^[10].强迫振动造成理论加工表面与实际加工表面不符, 影响加工工件的质量及精度.强迫振动导致的表面位置误差原理如图 2所示.
图 2(Fig. 2)

图 2 表面位置误差形成图Fig.2 Surface position error forming diagram

由理论可得铣削过程中y方向的切削力:

(9)

将F_y(φ)用傅里叶级数的形式展开后得

(10)

式中:a₀=[F_y(φ)]_mⁿ; a_n和b_n为傅里叶级数的系数.当系统为顺铣时，m=φ₁, n=π; 当系统为逆铣时, m=0, n=φ₁.
将F_y(φ)用离散傅里叶变换转化为F_y(ω)并与薄壁构件的频响函数H_y(ω)相乘后得Y(ω), 将Y(ω)通过傅里叶逆变换即可得到铣削系统的时域位移y(t).
2 基于Kriging方法的可靠性灵敏度分析Kriging方法作为一种改良线性回归分析技术, 它能利用已知点的实验数据模拟出未知点的实验数据^[11].Kriging模型包括线性组合和随机分布^[12].具体模型为

(11)

式中:f(x)是多项式函数; β是回归系数; f(x)提供模型的全局近似; z(x)提供局部近似.z(x)服从正态分布且协方差非零, 协方差为

(12)

式中:σ为标准差; R(θ, w, x)作为核函数, 对Kriging模型拟合的准确性影响非常大.目前常采用高斯函数作为核函数^[11].
如果想通过已知点求出未知点的响应值, 可以利用已知点响应值的线性组合来估计未知点的响应值, 即

(13)

式中:c为线性组合向量; Y为响应值向量.
真值与估计值之间的偏差为^[13-18]

(14)

当偏差均值为零时, 可以保证上述过程的无偏性, 即满足F^Tc-f(x)=0.
预测的均方误差为

(15)

式中: R为相关函数; r为相关向量.要求得φ(x)的最小值, 即要求解如下最优化问题^[19-23]:

(16)

最终得到估计值为

(17)

式中: β^*=(F^TR^-1F)^-1F^TR^-1Y.
由式(17)可知, 如果模型的参数β及核函数的表达式已知, 就可以得到未知点的响应值.
因为用数值模拟法计算动态表面位置误差需要耗费大量时间, 所以采用Kriging模型代替铣削动力学模型, 提高计算效率, 节约成本.
用Kriging方法进行可靠度灵敏度分析的具体步骤如图 3所示.
图 3(Fig. 3)

图 3 Kriging可靠性灵敏度分析流程图Fig.3 Kriging reliability sensitivity analysis flow chart

3 数值算例3.1 基于Kriging方法的铣削稳定性可靠性灵敏度分析仿真参数^[24]:铝合金2A12, 切向切削力系数k_tc=682.7 N/mm², 径向切削力系数k_rc=62 N/mm², 切向犁切力系数k_tp=40.4 N/mm, 径向犁切力系数k_rp=27 N/mm, 刀具齿数N=3, 刀具直径D=10 mm, 模态刚度k_y=1.7×10⁶ N/m, 固有频率ω_N=1 073 Hz, 阻尼比ζ=0.041, 径向侵入比r_ra=0.6, 逆铣.
将实际轴向切深为0.6 mm代入构造的模型中, 采用蒙特卡洛和Kriging方法分别计算得到的可靠度如图 4所示.
图 4(Fig. 4)

图 4 Kriging和MCS结果对比Fig.4 Comparison of Kriging and MCS results

由图 4可知, 选取n=6 000 r/min, 临界轴向切深为0.6 mm得到铣削系统的可靠度为0.842 2.基于Kriging方法的失效概率对输入变量进行均值灵敏度分析和方差灵敏度分析并与利用MCS方法进行可靠性灵敏度分析的结果进行对比, 如图 5和图 6所示.
图 5(Fig. 5)

图 5 Kriging和MCS均值灵敏度分析Fig.5 Kriging and MCS mean sensitivity analysis

图 6(Fig. 6)

图 6 Kriging和MCS方差灵敏度分析Fig.6 Kriging and MCS variance sensitivity analysis

由图 5和图 6可知, 薄壁构件铣削系统的阻尼比、y方向的刚度的增加会增加系统的可靠度, 切向力系数的增加会降低系统的可靠度, 固有频率和径向力系数对系统的影响较小.
3.2 基于Kriging方法的铣削表面位置误差可靠性灵敏度分析仿真参数除3.1节外还包括轴向切深a_ax=1 mm, 每齿进给量f_t=0.01 mm.
选取n=11 000 r/min及表面位置误差许用值SLE=65 μm, 得到铣削系统的可靠度为0.817 8.基于Kriging方法的失效概率对输入变量进行均值灵敏度分析和方差灵敏度分析，并与利用MCS方法进行可靠性灵敏度分析的结果进行对比, 如图 7和图 8所示.
图 7(Fig. 7)

图 7 Kriging和MCS均值灵敏度分析Fig.7 Kriging and MCS mean sensitivity analysis

图 8(Fig. 8)

图 8 Kriging和MCS方差灵敏度分析Fig.8 Kriging and MCS variance sensitivity analysis

由图 7和图 8可知, 薄壁构件铣削系统的阻尼比、y方向刚度的增加会增加系统的可靠度, 切向力系数、固有频率、径向侵入比以及轴向切深的增加会降低系统的可靠度, 径向力系数和每齿进给量对系统的影响较小.
4 结论随着薄壁构件铣削系统的阻尼比及y向刚度的增加, 系统的可靠度增加; 随着切向力系数、固有频率、径向侵入比及轴向切深的增加, 系统的可靠度降低; 固有频率、径向力系数和每齿进给量对系统的影响较小.
参考文献

[1]	卢晓红, 王凤晨, 王华, 等. 铣削过程颤振稳定性分析的研究进展[J].振动与冲击, 2016, 35(1): 74–82. ( Lu Xiao-hong, Wang Feng-chen, Wang Hua, et al. Research progress on flutter stability analysis in milling process[J].Vibration and Impact, 2016, 35(1): 74–82.)
[2]	Wu G, Luan W. Self-excited vibration of driveline for vehicle launch[J].Journal of Beijing Institute of Technology, 2013(3): 330–336.
[3]	Ding H, Ding Y. On time-domain methods for milling stability analysis[J].Chinese Science Bulletin, 2012(33): 4336–4345.
[4]	Jemielniak K, Widota A. Suppression of self-excited vibration by the spindle speed variation method[J].International Journal of Machine Tool Design & Research, 1984, 24(3): 207–214.
[5]	Li H, Shin Y C. A comprehensive dynamic end milling simulation model[J].Journal of Manufacturing Science & Engineering, 2006, 128(1): 86–95.
[6]	Thevenot V, Arnaud L, Dessein G, et al. Influence of material removal on the dynamic behavior of thin-walled structures in peripheral milling[J].Machining Science and Technology, 2006, 10(3): 276–283.
[7]	Tang A, Liu Z. Three-dimensional stability lobe and maximum material removal rate in end milling of thin-walled plate[J].The International Journal of Advanced Manufacturing Technology, 2009(1/2): 33–39.
[8]	Bravo U, Altuzarra O. Stability limits of milling considering the flexibility of the workpiece and the machine[J].International Journal of Machine Tools & Manufacture, 2005, 45(15): 1669–1680.
[9]	Seguy S. Toolpath dependent stability lobes for the milling of thin-walled parts[J].International Journal of Machining & Machinability of Materials, 2008, 4(3/4): 261–274.
[10]	Kaymaz I. Application of Kriging method to structural reliability problems[J].Structural Safety, 2005, 27(2): 133–151.DOI:10.1016/j.strusafe.2004.09.001
[11]	Sun Z, Wang J, Li R, et al. LIF:a new Kriging based learning function and its application to structural reliability analysis[J].Reliability Engineering & System Safety, 2017, 157: 152–165.
[12]	Schmitz T L, Mann B P. Closed-form solutions for surface location error in milling[J].International Journal of Machine Tools and Manufacture, 2006, 46(12/13): 1369–1377.
[13]	Echard B, Gayton N, Lemaire M. AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[J].Structural Safety, 2011, 33(2): 145–154.DOI:10.1016/j.strusafe.2011.01.002
[14]	Zhang W, Dai B, Liu Z, et al. Unconfined seepage analysis using moving Kriging mesh-free method with Monte Carlo integration[J].Transport in Porous Media, 2016, 116(1): 1–18.
[15]	Ma J, Ren Z, Zhao G, et al. A new reliability analysis method combining adaptive Kriging with weight index Monte Carlo simulation[J].IEEE Transactions on Magnetics, 2018, 54(3): 1–4.
[16]	Gong Q, Zhang J G, Su D. Reliability simulation combining Kriging and Monte Carlo radius-outside importance sampling in space structure latch[J].Applied Mechanics and Materials, 2012, 166/167/168/169: 1872–1878.
[17]	Emery X. Two ordinary Kriging approaches to predicting block grade distributions[J].Mathematical Geology, 2006, 38(7): 801–819.
[18]	Ginsbourger D, Rosspopoff B, Pirot G, et al. Distance-based Kriging relying on proxy simulations for inverse conditioning[J].Advances in Water Resources, 2013, 52: 275–291.DOI:10.1016/j.advwatres.2012.11.019
[19]	Emery X. Simple and ordinary multigaussian Kriging for estimating recoverable reserves[J].Mathematical Geology, 2005, 37(3): 295–319.
[20]	Janusevskis J, Riche R L. Simultaneous Kriging-based estimation and optimization of mean response[J].Journal of Global Optimization, 2013, 55(2): 313–336.DOI:10.1007/s10898-011-9836-5
[21]	Kleijnen W, Jack P C. Simulation-optimization via Kriging and bootstrapping: a survey[J].Journal of Simulation, 2014, 8(4): 241–250.DOI:10.1057/jos.2014.4
[22]	Zhang L, Lu Z, Wang P. Efficient structural reliability analysis method based on advanced Kriging model[J].Applied Mathematical Modeling, 2015, 39(2): 781–793.
[23]	Ng S H, Yin J. Bayesian Kriging analysis and design for stochastic simulations[J].ACM Transactions on Modeling and Computer Simulation, 2012, 22(3): 1–26.
[24]	刘彪.薄壁件铣削稳定性与动态加工误差研究[D].北京: 北京理工大学, 2015. ( Liu Biao.Study on stability and dynamic machining error of thin-walled parts[D]. Beijing: Beijing Institute of Technology, 2015.http://cdmd.cnki.com.cn/Article/CDMD-10007-1016716946.htm)