2. 沈阳化工大学 装备可靠性研究所, 辽宁 沈阳 110142
2. Equipment Reliability Institute, Shenyang University of Chemical Technology, Shenyang 110142, China.
现代加工对切削精度和效率的要求越来越高, 且高性能切削都以稳定切削为前提, 而车削颤振是影响车削稳定性的重要因素.车削过程中的刀具磨损问题是不可避免的, 由于刀具磨损使得切削力不断增大, 切削系数不断增加[1-5], 影响车削颤振稳定性的预测, 因此在研究车削颤振稳定性时考虑刀具磨损非常必要.Campocassoa等[1]证明了刀刃几何形状和切削角度随刀具磨损的变化, 使切削力系数随之变化.Li等[2]发现刀具磨损使切削力在各方向上均有增加.Malakooti等[3]研究了刀具磨损与刀具寿命的关系, 将刀具磨损分为初级磨损、正常磨损和急剧磨损.Naves等[4]研究了AISI316不锈钢在加工过程中刀具磨损量与加工时间的关系.Cascón等[5]发现在车削过程中的切削力随着切削时间的增加而逐渐增大.
由于受到加工系统环境的变化及加工材料表面分布不均匀等因素的影响, 加工系统的结构参数和材料切削性能参数等具有不确定性, 无法准确预测分析车削颤振稳定性, 因此建立车削加工过程中的颤振时变可靠性概率模型.在时变可靠性分析方面已有相当多的研究, 并且广泛应用于工程实际[6-11].Schueller[6]和Sudret等[7]通过离散随机过程, 采用了Monte Carlo方法(MCS)计算时变可靠度.Singh等[8]提出了采用重要抽样方法近似计算超越率(即失效率)的方法.González-Fernández[9]基于MCS方法, 提出了序列交叉熵MCS方法评估时变可靠性指标.Andrieu-Renaud等[10]提出了PHI2方法, 通过离散时间区间的可靠度模型计算时变可靠度.Hu等[11]提出了first order simulation approach(FOSA)方法, 采用高斯随机过程表示响应的随机过程, 进而计算时变可靠度.本文考虑到车削加工系统颤振时变可靠度模型的非线性高且失效概率小的特点, 采用代理模型与方差缩减这两类相结合的方法进行可靠性分析.
代理模型类可靠性分析方法能减少计算可靠度中调用非线性程度高的极限状态函数的次数, 而方差缩减类可靠性分析方法能在保证精度的前提下减小样本数, 有效解决小失效概率问题, 因此代理模型与方差缩减这两类方法的结合可以解决车削加工颤振时变可靠度的计算问题.近年来, 代理模型与方差缩减相结合的方法[12-14]也逐渐增多.Bourinet等[12]提出了基于支持向量机的子集模拟方法, Papadopoulos等[13]提出了神经网络与子集模拟相结合的方法, Echard等[14]将重要抽样方法与Kriging模型相结合, 提出了AK-IS方法.这些方法对于小失效概率问题的求解效率皆有极大程度上的提高, 但对于解决时变可靠度问题, 依然没有较准确高效的计算方法.
本文针对考虑刀具磨损的车削颤振时变可靠性问题, 采用Gamma过程描述了切削力系数的变化过程, 建立了车削加工系统的颤振时变可靠性模型.综合代理模型和方差缩减的优点, 并结合首次超越方法, 提出了主动学习的Kriging模型与基于首次超越的子集模拟方法相结合的时变可靠性分析方法——AK-SSFP.根据主动学习原则, 采用学习函数U选择最佳样本点优化Kriging模型, 代替真实的极限状态函数, 再应用基于首次超越的子集模拟时变可靠性方法计算车削加工系统的时变可靠度.
1 考虑刀具磨损的车削颤振稳定性分析 1.1 包含时变切削力系数的车削加工动力学模型由刀具磨损引起的切削力系数的变化是一个时变过程, 且切削力系数具有随机性.本文采用Gamma过程描述切削力系数随时间的变化过程.定义ΔKs(t)为0时刻到t时刻切削力系数的增量.合力切削力系数的表达式为
(1) |
式中:Ks0为刀具磨损前合力的切削力系数, 由切向切削力系数Kt和径向切削力Kr组成, 其表达式为
根据Cascón等[5]的实验数据拟合描述切削力系数增量变化的Gamma过程参数为b=0.615 8, c=0.921 8, β=0.027 41, 即Gamma过程的形状参数和尺寸参数分别为
(2) |
图 1为车削加工系统的单自由度(SDoF)再生型颤振动力学模型简图.由于刀具系统为车削加工系统的主振系统, 并且切削宽度的动态变化仅由再生型颤振引起, 车削加工系统的动力学微分方程为
(3) |
式中:m, c和k分别为车削系统的等效质量kg、等效阻尼(N · s/m)和等效刚度(N/m); ΔF(t)为动态切削力合力; φ为切削力合力与刀具振动方向的夹角; α为主振动方向与刀具振动方向的夹角.
动态切削力合力由切削宽度的动态变化决定[15], 因此动态切削力合力的表达式为
(4) |
式中:ap为切削宽度(mm); h(t)为切削厚度变化量(mm); η为车削加工重叠系数, 0 < η < 1;y(t)为当前时刻车削振动位移(mm); y(t-T)为上一转时车削振动位移(mm); T为主轴旋转周期(s), T=60/Ω,Ω为主轴转速(r/min).
由图 1可知, 刀具振动方向上振动位移与主振动方向上振动位移的关系为y(t)=u(t) · cosα, 并假设车削加工振动位移表示为y(t)=Asin(ωt+φ0), 且初始条件均为0, 将考虑刀具磨损的车削加工动力学微分方程进行Laplace变换:
(5) |
式中:s=σ+iω为特征方程的根;K=-ωn2Ksapu/k; ωn为车削加工系统的固有频率, ωn2=k/m; ξ为车削加工系统的阻尼比, ξ=c/2mωn; u为方向系数, u=cos(φ-α)cosα.
1.2 车削颤振时变稳定性根据Nyquist稳定判据, 式(5)所示特征方程根的实部σ=0时, 系统处于稳定与不稳定的临界状态[16], 即s=iω, 解得临界状态下的时变极限切削宽度为
(6) |
对应的主轴转速为
(7) |
式中:λ=ω/ωn, 由于系统振动频率ω略大于系统固有频率ωn, λ>1;j=0, 1, 2, …, 为叶瓣数.
2 车削颤振时变可靠性分析 2.1 车削颤振时变可靠性模型在实际车削加工过程中, 由于测量误差和环境等因素, 车削颤振稳定性分析需要考虑参数的随机性, 车削颤振可靠度为车削加工系统不发生颤振的概率, 即车削加工极限切宽大于给定切宽的概率.
车削颤振时变可靠性极限状态函数为
(8) |
式中:ap为给定切削宽度; S (t)=[X, Y (t)]为车削加工系统参数的样本点, X为关于车削加工系统参数x的随机变量, x=[k, c, m, φ, α, η]T, Y (t)=(Ks(t), Ω(t))为切削力系数和主轴转速随机过程.
车削加工系统从0到t时刻内发生颤振可表示为
(9) |
因此0到t时刻的系统累积失效概率为PF(0, t)=Prob(E), 则车削颤振时变可靠性模型为
(10) |
在分析车削颤振时变可靠性时, 将时间区间[0, t]平均分为nt-1个时间段, 各时间节点可表示为ti=(i-1)Δt, nt为时间节点数, Δt=t/(nt-1)为时间间隔.目前, MCS为解决时变可靠度问题最精确的方法, 若失效概率在10-k数量级上, 为达到相应的精度, 则需要10k+2~10k+3个随机样本[17].采用MCS方法分析车削颤振时变可靠性问题时, 由于车削颤振时变可靠性极限状态函数为非线性程度较高的函数, MCS方法所需的计算量相当大, 而采用主动学习的Kriging模型来代替真实的极限状态函数能够减少极限状态函数的调用次数, 提高计算效率.由于车削加工系统各时间段的失效概率较小, 采用子集模拟方法高效解决小失效概率问题.针对车削加工时变系统, 本文在计算累积失效概率时, 采用了结合首次超越方法[10]和子集模拟的时变可靠性方法.综合子集模拟时变可靠性方法和主动学习的Kriging模型提出了AK-SSFP方法, 详细计算流程如下:
1) 将时间区间[0, t]平均分为nt-1个时间段, 各时间节点为ti, i=1, 2, …, nt, 生成N个离散在各个时间点上的车削参数候选样本S.
2) 采用Latin hypercube sampling(LHS)方法生成m个样本点D, 即初始车削参数DoE.计算各样本点Dim在各时间点上的极限状态函数值Gim.本文取m为12, 通过主动学习过程更新DoE.
3) 根据时间点ti处的Dti和Gti构建该时间节点上的Kriging模型, 相关模型选为Gaussian函数, 回归多项式取常数1.计算S中各样本点SiN在时间ti处的预测值ĜiN, ti=ĝ(SiN, ti)及学习函数值U(SiN, ti).
4) 在时间点ti处, 对于所有样本点有min(U(SiN, ti))≥2, 则停止主动学习; 否则, 将最佳样本点Sti*和极限状态函数值g(Sti*)加入到Dti和Gti中, 返回步骤3), 改进时间点ti处的Kriging模型.以此类推, 构建所有时间点处的Kriging模型.
5) 基于上述时间点ti-1和ti处的Kriging模型, 采用子集模拟方法计算[ti-1, ti]时间段的失效概率Pf,ti.
6) 若ti=t, 输出所有时间段的失效概率Pf, 并计算其变异系数δpf; 否则, i=i+1, 返回步骤5).
7) 预设[δ]=5 % [18], 如果所有时间段失效概率的变异系数都不大于[δ], 即max(δPf)≤[δ], 算法结束, 输出失效概率Pf; 反之, 增加车削参数候选样本的数量, 返回步骤3).
8) 由失效概率Pf计算车削加工系统各时间段[0, ti]的累积失效概率PF, 进而得到可靠度R.
AK-SSFP结合了子集模拟方法和主动学习Kriging模型的优势.对于车削加工系统颤振时变可靠性问题, 主动学习Kriging模型通过少量的样本建立代理模型, 代替真实的高非线性极限状态函数, 进而采用子集模拟时变可靠性方法将小失效概率转化为较大概率的乘积, 从而达到了在保证高精度的前提下降低计算量、提高运算效率的目的.
3 算例分析 3.1 车削加工系统在车削加工系统[19]中, 以50CrMo4硬化钢为工件材料[5], 切向切削力系数Kt=2 560 N/mm2, 径向切削力系数Kr=1 625 N/mm2, 车削加工系统的性能参数分布特征如表 1所示.主轴转速在车削加工时也不是恒定的, 可以用随机过程描述, 均值为2 500 r/min, 自相关函数为
(11) |
式中, 相关系数l=60.
3.2 颤振时变稳定性叶瓣图将切向切削力系数Kt和径向切削力系数Kr代入式(1), 得到合力切削力系数Ks(t)的拟合结果如图 2所示.将表 1中的参数代入式(6)和式(7), 可得车削加工系统的切削颤振时变稳定性叶瓣图, 如图 3所示.
根据2.1节中车削颤振时变可靠性模型, 采用本文提出的AK-SSFP方法, 并与MCS和AK-MCS方法进行对比.AK-SSFP方法的初始样本数取为103, AK-MCS方法的初始样本数取为104, MCS方法的样本数取为106.本文中的计算均在Intel i5处理器下运行, 并且运行5次取平均值以降低计算结果的不确定性.切削时长为35 min, 时间间隔为1 min, 在给定实际切宽ap=1.6 mm, 主轴转速的均值为2 500 r/min的情况下, 采用几种方法计算车削加工系统颤振时变可靠度的结果如图 4所示.
由图 4可知, 车削加工系统颤振时变可靠度在0时刻为1, 随着切削时间的增加, 可靠度逐渐减小.而且AK-MCS方法和AK-SSFP方法的计算结果与MCS方法计算结果的一致性都很好.
表 2为在主轴转速为2 500 r/min的情况下, 采用3种时变可靠性方法计算车削加工系统颤振时变可靠度, 计算结果之间的相对误差, 以及运算时调用真实极限状态函数的次数NMCS, NAK-MCS, NAK-SSFP和运算时长tAK-MCS, tAK-SSFP.
由表 2可知, 采用AK-MCS方法计算的车削颤振时变可靠度与MCS方法计算结果之间的平均相对误差εAK-MCS=0.548 62 %, 在工程允许范围内, 而AK-MCS方法所需的NAK-MCS=17 942, 明显少于MCS方法所需的NMCS=3.6×107, 大幅提高了计算效率, 但AK-MCS方法依然需要大量的随机样本作为候选样本, 所以其运算时长较长.虽然AK-SSFP方法的平均相对误差εAK-SSFP=2.166 %,与AK-MCS方法相比有所增大, 但依然在工程允许范围内, 而AK-SSFP方法的运算时长tAK-SSFP= 9776.78 s, 与AK-MCS方法的运算时长tAK-MCS=30 254.27 s相比, 有很大程度上的降低.因此, 采用AK-SSFP方法解决车削颤振时变可靠度问题, 在保证计算精度的前提下大幅提高了运算效率.
4 结论1) 本文针对考虑刀具磨损的车削加工系统, 采用Gamma过程描述切削力系数的变化, 分析颤振时变稳定性, 进而建立了车削加工系统的颤振时变可靠性模型.
2) 本文结合子集模拟方法和主动学习Kriging模型提出用AK-SSFP时变可靠性分析方法计算车削颤振时变可靠度.
3) 通过对车削加工实例进行分析, 验证了AK-SSFP方法与MCS方法计算结果的一致性, 且在减少调用真实极限状态函数次数的基础上, 进一步降低了所需样本数, 缩减了运算时间, 极大地提高了计算效率, 有效地解决了车削颤振这类高非线性极限状态函数且失效概率较小的实际工程时变可靠度计算问题.
[1] |
Campocasso S, Costes J P, Fromentin G, et al.
A generalised geometrical model of turning operations for cutting force modeling using edge discretisation[J]. Applied Mathematical Modeling, 2015, 39: 6612–6630.
DOI:10.1016/j.apm.2015.02.008 |
[2] |
Li H, Lai X, Li C, et al.
Modeling and experimental analysis of the effects of tool wear, minimum chip thickness and micro tool geometry on the surface roughness in micro-end-milling[J]. Journal of Micromechanics and Microengineering, 2007, 18(2): 025006.
DOI:10.1088/0960-1317/18/2/025006 |
[3] |
Malakooti B, Wang J, Tandler E C.
A sensor-based accelerated approach for multi-attribute machinability and tool life evaluation[J]. The International Journal of Production Research, 1990, 28(12): 2373–2392.
DOI:10.1080/00207549008942872 |
[4] |
Naves V T G, Da Silva M B, Da Silva F J.
Evaluation of the effect of application of cutting fluid at high pressure on tool wear during turning operation of AISI 316 austenitic stainless steel[J]. Wear, 2013, 302(1): 1201–1208.
|
[5] |
Cascón I, Sarasua J A.
Mechanistic model for prediction of cutting forces in turning of non-axisymmetric parts[J]. Procedia CIRP, 2015, 31: 435–440.
DOI:10.1016/j.procir.2015.03.084 |
[6] |
Schueller G I.
A state-of-the-art report on computational stochastic mechanics[J]. Probabilistic Engineering Mechanics, 1997, 12(4): 197–321.
DOI:10.1016/S0266-8920(97)00003-9 |
[7] |
Sudret B, Der Kiureghian A.
Stochastic finite element methods and reliability:a state-of-the-art report[J]. University of California Berkeley, 2000, 11: 189.
|
[8] |
Singh A, Mourelatos Z, Nikolaidis E.
Time-dependent reliability of random dynamic systems using time-series modeling and importance sampling[J]. SAE International Journal of Materials and Manufacturing, 2011, 4(1): 929–946.
DOI:10.4271/2011-01-0728 |
[9] |
González-Fernández R A.
Reliability assessment of time-dependent systems via sequential cross-entropy Monte Carlo simulation[J]. IEEE Transactions on Power Systems, 2011, 26(4): 2381–2389.
DOI:10.1109/TPWRS.2011.2112785 |
[10] |
Andrieu-Renaud C, Sudret B, Lemaire M.
The PHI2 method:a way to compute time-variant reliability[J]. Reliability Engineering and System Safety, 2004, 84(1): 75–86.
DOI:10.1016/j.ress.2003.10.005 |
[11] |
Hu Z, Du X.
First order reliability method for time-variant problems using series expansions[J]. Structural and Multidisciplinary Optimization, 2015, 51(1): 1–21.
|
[12] |
Bourinet J M, Deheeger F, Lemaire M.
Assessing small failure probabilities by combined subset simulation and support vector machines[J]. Structural Safety, 2011, 33(6): 343–353.
DOI:10.1016/j.strusafe.2011.06.001 |
[13] |
Papadopoulos V, Giovanis D G, Lagaros N D, et al.
Accelerated subset simulation with neural networks for reliability analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 223: 70–80.
|
[14] |
Echard B, Gayton N, Lemaire M, et al.
A combined importance sampling and Kriging reliability method for small failure probabilities with time-demanding numerical models[J]. Reliability Engineering & System Safety, 2013, 111: 232–240.
|
[15] |
Fu X, Zheng S.
New approach in dynamics of regenerative chatter research of turning[J]. Communications in Nonlinear Science & Numerical Simulation, 2014, 19(11): 4013–4023.
|
[16] |
胡寿松.
自动控制原理[M]. 北京: 科学出版社, 2007.
( Hu Shou-song. Principle of automatic control[M]. Beijing: Science Press, 2007. ) |
[17] |
Baumgärtner A, Binder K.
Applications of the Monte Carlo method in statistical physics[M]. Berlin: Springer-Verlag, 1987.
|
[18] |
Bichon B J, Eldred M S, Swiler L P, et al.
Efficient global reliability analysis for nonlinear implicit performance functions[J]. AIAA Journal, 2008, 46(10): 2459–2468.
DOI:10.2514/1.34321 |
[19] |
Liu Y, Li T X, Liu K, et al.
Chatter reliability prediction of turning process system with uncertainties[J]. Mechanical Systems and Signal Processing, 2016, 66/67: 232–247.
DOI:10.1016/j.ymssp.2015.06.030 |