李朕均, 赵春雨, 闻邦椿, 卢泽宸
东北大学 机械工程与自动化学院, 辽宁 沈阳 110819
收稿日期：2018-09-18
基金项目：国家自然科学基金资助项目(51775094)。
作者简介：李朕均(1992-)，男，吉林四平人，东北大学博士研究生;
赵春雨(1963-)，男，辽宁黑山人，东北大学教授，博士生导师;
闻邦椿(1930-), 男，浙江温岭人，东北大学教授，博士生导师，中国科学院院士。
通信作者：ZHAO Chun-yu, E-mail: chyzhao@mail.neu.edu.cn。

摘要：数控机床进给系统具有热源多且各热源发热率时变性强的特点.但以往研究对各热源发热率的时变性关注较小, 因此很难实现丝杠温升和热误差的精确预测.本文将Monte Carlo(MC)模拟集成FEM方法捕获到进给系统各热源的发热率, 进而确定了各热源发热率占系统总发热率的比例——分配比随系统运行时间的变化规律.然后, 基于进给系统各热源的发热率分配比, 提出利用进给系统伺服电机的力矩电流计算各热源发热率的方法.最后, 基于一维热传导方程的有限差分法, 提出了丝杠温升和热误差的数值预测模型, 并利用试验验证了模型预测的精确性.
关键词：数控机床滚珠丝杠进给系统温度场热传导热误差
Heat Source Rate Identification and Thermal Error Predictions of Ball Screw Feed Drive System for CNC Machine Tools
LI Zhen-jun, ZHAO Chun-yu, WEN Bang-chun, LU Ze-chen
School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Abstract: The ball screw feed drive system of CNC(computerized numerical control)machine tools is characterized as multi-heat sources with strong time-varying heat rate. But the past work was scarely focused on the time-varying of the heat sources. This fact results in inaccuracy of predicting temperature rise and thermal error of the screw. In this paper, the heat rates of the multi-heat sources were captured by virtue of the FEM integrated with Monte Carlo method. Consequently, relationship of the ratio of each heat source to the total heat rate of the system with the running time was obtained. Then, based on heat ratios of the heat sources, the heat rates of the heat sources were determined using torque current of the servo motor in the ball screw feed drive system. Finally, numeric prediction model of temperature rise and thermal error of the screw was put forward based on finite difference method of one-dimension heat conduction equation, and the accuracy of the prediction method was verified by experimental results.
Key words: CNC machine toolball screw feed drive systemtemperature fieldheat conductionthermal error
滚珠丝杠进给驱动系统作为数控机床精密传动和定位的关键部件, 对机床定位精度起着至关重要的作用.然而在工作过程中, 由于各运动副接触面间的摩擦导致温升, 引起结构热变形影响刀具定位精度^[1-2].以往研究表明, 进给系统热变形引起的加工误差占总误差的40%到70%^[3]，因此, 进给系统热误差问题的研究具有重要意义.
进给系统的热误差预测及补偿控制或采用附加措施减小丝杠温升是降低热误差及提高加工精度的重要手段^{[1, 2, 4]}.进给系统热源识别及温度分布预测是热误差预测的基础, 也是本领域国内外学者的研究热点.如基于系统大量试验数据, 利用回归分析和神经网络实现温度分布与热误差的模型预测^[5-6]; 基于集中热容量模型, Yun等利用GEA算法实现了进给系统温升与热误差的预测^[3]; 基于丝杠一维热传导方程的解析解, 陈诚等提出了通过辨识模型参数实现丝杠温度场预测方法^[7].但是, 由于进给系统热源较多, 且随进给速度和环境的变化而变化^[1], 导致各热源发热率及系统温升分布具有较大的时变性.因此, 给进给系统的温升及热误差的精确预测带来较大的困难.
针对上述问题, 本文通过集成Monte Carlo(MC)模拟和有限元的方法确定进给系统各热源发热率的分配比, 进而利用伺服电机的力矩电流辨识进给系统各热源发热率, 提出了一种全新的基于伺服电机力矩电流的丝杠热误差预测模型, 为数控机床进给系统热误差的实时补偿提供依据.
1 进给系统热特性试验试验系统如图 1所示, 试验对象为HTC2050i数控车床的x轴滚珠丝杠进给系统.数控系统为FANUC 0i Mate-TD.工作台质量为500 kg, 该滚珠丝杠型号为3210, 公称直径为32 mm, 导程为10 mm, 行程为220 mm, 导轨形式为滚动导轨.丝杠由2个轴承支撑, 其中靠近伺服电机的轴承为固定支撑, 另一端为支持轴承.PC104计算机通过网线接口与数控系统相连.FOCAS函数接口程序运行在PC104工控机上, 实现进给系统伺服电机力矩电流信号的读取和存储.试验过程中进给系统空载运行, 其进给速度分别为5, 10, 15 m/min; 试验时间为1 h.
图 1(Fig. 1)

图 1 试验设备Fig.1 Experimental equipment

为了分析进给系统热源发热特性, 应合理地选择系统温度和定位误差检测点.温度检测传感器配置尽量靠近热源, 以准确地反映各热源温度的变化.系统温度和定位误差的测点分布如图 2所示.
图 2(Fig. 2)

图 2 试验测点分布Fig.2 Distribution of testing points

1) 使用磁性吸附式热电偶检测温度, 其检测点为:电机端轴承座的温度T₁、下端轴承座温度T₂、螺母法兰温度T₃和导轨滑块温度T₄.采样间隔设定为48 ms.
2) 丝杠表面点的温度使用红外热成像仪进行全行程检测, 每往复运动10 min检测一次.热成像的数据处理过程按等间距拾取丝杠表面温度点, 如图 2中P₁~P₆, 其点间距离为44 mm.
3) 使用雷尼绍XL80激光干涉仪测量定位误差, 其测量点仍选定为P₁~P₆点, 工作台每往复运动10 min检测一次.
2 集成Monte Carlo法的有限元法发热率分配关系辨识2.1 集成Monte Carlo法的有限元仿真计算选择x轴进给系统为研究对象, 利用ANSYS软件建立其有限元模型.其结构由z轴滑板和x轴的滑板、轴承系统及丝杠组成.单元类型选用Solid90, 各运动副接触面选为TARGE170和CONTA174接触单元, 如图 3所示.
图 3(Fig. 3)

图 3 进给系统有限元模型Fig.3 Finite element model of feed system

本文用x₁, x₂, x₃和x₄分别表示螺母、电机端轴承、下端轴承、滑块发热率.为了求得机床进给系统的发热率, 建立目标函数如下:

(1)

其中:T_ij是第j个采样时间步长的给定点i的温度; 上标EM和MC分别代表试验测量值和有限元模拟仿真值.
蒙特卡洛模拟确定各热源的热流密度发热率x_m.若设x_m的下界和上界分别为x_Lm和x_Hm, 则利用随机函数产生k组热流密度样本值:

(2)

分别将每组热流密度值代入有限元模型进行模拟计算, 提取温度检测点与检测时间对应的温度值T_lij^MC.分别将每组模拟温度T_lij^MC和试验检测结果T_ij^EM代入式(1)计算目标函数F_l, 找到k次模拟最小目标函数对应的热流密度, 记为x′₁, x′₂, x′₃, x′₄.如果满足计算精度min(F_l)≤ε, 输出计算结果, 结束模拟仿真计算; 否则按二分法更改上下限值:若x′_m－x_Lm < 0.5δ, 则x_Hm=x_Hm－0.5δ; 否则x_Lm=x_Lm+0.5δ, m=1, 2, 3, 4.然后, 重复上述模拟过程.
由于数控机床进给系统的滚珠丝杠螺母和滑块在工作过程中属于往复运动状态, 为了真实反映实际温度场, 本文对螺母和滑块移动热源, 采用移动热源加载方式, 即丝杠螺母和导轨热源位置随工作台位置变化而变化, 其求解过程使用APDL命令流编程实现.
2.2 各接触面发热率的辨识及分配比利用上述蒙特卡洛模拟对进给系统的各热源热流密度进行辨识.其中, 样本数量k=70;计算精度ε=0.01;各热源热流密度下限为0, 上限取进给系统消耗的总功率P_总.图 4为进给速度5 m/min下各热源发热率的计算结果.由图 4可以看出, 各热源发热率在开始时急速下降, 随系统运行时间增加逐渐达到稳态值, 这种变化与系统温度热平衡过程相似^[6].其原因是初始阶段各接触面温度较低, 润滑油膜黏度大, 接触面间摩擦热流密度较高; 随接触面温度升高, 润滑油膜黏度减小致使接触面间的发热率下降.文献[8]利用龚帕兹模型描述各热源的发热率随进给时间的变化规律, 与本文的计算结果吻合.
图 4(Fig. 4)

图 4 各热源发热率Fig.4 Heating rate of each heat source

文献[9]分析了进给速率对滚珠丝杠温度场的影响, 证明热源发热率随进给速率线性增加.因此, 可利用不同进给速度下的试验数据对各热源发热率辨识, 并对辨识结果进行非线性回归分析, 确定各热源发热率随进给速度和时间的变化关系：

(3)

其中:v为进给速度; t为进给时间; K_i和b_i为各热源的非线性回归系数.
对不同进给速度下各个热源发热率的辨识结果P_i进行分配比计算, 即各热源发热率与总发热率的比:

(4)

其中, w_i为各热源的分配比.
将式(3)代入式(4), 得

(5)

利用式(5)可以得到不同进给速度工况下各热源发热率分配比随时间变化趋势.尽管各热源发热率随进给速度线性增加, 但是其分配比不受进给速度影响, 即不同进给速度下发热率的分配关系不随进给速度发生变化, 与确定的相关系数和加工时间有关.
2.3 基于伺服电机力矩电流的发热率辨识在试验过程中, 利用FOCAS函数接口程序获得伺服电机的力矩电流.伺服电机输出转矩与力矩电流的关系为

(6)

式中:i为驱动电机电流; K_t为驱动电机的转矩系数; T为驱动电机力矩.
滚珠丝杠进给系统消耗总功率为

(7)

其中:P_总为总摩擦热功率; ω为角速度; m_t为刀架的质量; g为重力加速度; α₀为滚珠丝杠驱动系统与水平面相对的倾斜角; v为刀架的进给速度.
将蒙特卡洛模拟辨识的各热源发热率相加求得进给系统的总摩擦热功率, 图 5为进给速度5 m/min的有限元辨识与由力矩电流计算总发热率的对比, 可以看出两者十分吻合, 随时间的变化趋势也相同.因此, 利用式(5)所得的发热率分配关系与试验测得的伺服电机力矩电流, 可以计算各个热源发热率:
图 5(Fig. 5)

图 5 发热率对比Fig.5 Comparison of heating rate

(8)

3 热误差的预测模型滚珠丝杠的轴向热变形是导致机床加工误差的主要原因.在此忽略丝杠的径向变形, 滚珠丝杠下端固定, 上端自由, 根据传热学理论可以简化为一维杆导热体^[10].如图 6所示, 使用空间尺度S把丝杠轴离散为N段, 即L=j·S, j=1, 2, …, N; 同时用时间尺度τ划分系统运行时间为M段, 即t=i·τ, i=1, 2, …, M, 进行区域离散, 建立节点离散方程.
图 6(Fig. 6)

图 6 滚珠丝杠离散化示意图Fig.6 Schematic diagram of discretization of ball screw

在系统运行过程中, 两端轴承热源向丝杠注入热量进行热传导; 同时螺母运行至丝杠某接触点时产生摩擦热; 另外, 丝杠与空气接触段S_j产生对流换热, 与紧领的S_j-1和S_j+1两段进行热传导.热平衡方程为

(9)

其中:ΔT_jⁱ为S_j段在第i时段的温度变化; Q_1－1ⁱ为下端轴承产生的总热量; Q_2－Nⁱ为电机端轴承在第i时段产生的总热量; Q_3－jⁱ为移动螺母在第i时段产生的总摩擦热; Q_L－jⁱ是S_j段沿轴向热传导量; Q_out－jⁱ为S_j段与空气换热量; c为丝杠材料比热容; m为每段的质量.
热源发热量Q_1－1ⁱ, Q_2－Nⁱ, Q_3－jⁱ为在时间尺度τ内产生的总热量, 热源产生的总热量为

(10)

Q_L－jⁱ是在时间尺度τ内进行的热传导, 计算式如下:

(11)

其中:λ为导热率; A_c为丝杠截面积; T_jⁱ为第j段在i时刻温度; T_j－1ⁱ为j-1段在i时刻温度; T_j+1ⁱ为j+1段在i时刻温度.
滚珠丝杠与工作环境的热交换是散热方式, 在时间尺度τ内, 丝杠与空气换热量Q_out－jⁱ为

(12)

其中:h_v为丝杠与空气对流换热系数; A为每段丝杠与空气接触表面积; T_a为工作环境温度.
对ΔT_jⁱ进行求和计算, 为第j段丝杠在系统工作期间内温度变化量, 则热伸长量为

(13)

其中:ΔL_j为丝杠j段的热伸长量; ΔT_j为丝杠j段的温度变化量; α为丝杠轴热膨胀系数.
在进给系统工作过程中, 工作台运行至位置j时的热误差为

(14)

4 结果分析在对进给系统发热率辨识完成后, 应用本文预测数学模型计算滚珠丝杠温升与红外热成像仪检测结果的对比, 如图 7所示.检测点间距44 mm, 检测周期为10 min.由图 7可知, 模型预测与试验检测温升吻合得较好, 证明了利用滚珠丝杠进给系统热源强度时变特性辨识热源发热率是可行的, 继而证明了数学模型预测实际滚珠丝杠温度场的有效性.
图 7(Fig. 7)

图 7 丝杠温升随时间的变化Fig.7 Variations of screw temperature rise with time

利用式(13)和式(14), 可以得到工作台处于丝杠给定位置的热误差, 其计算结果与试验检测热误差对比如图 8所示.由图可以看出, 热误差预测模型的计算结果与实际试验结果残差在1 μm之内, 证明了该热误差数学模型能够很好预测实际热误差.
图 8(Fig. 8)

图 8 热误差预测与试验结果比较Fig.8 Comparison of thermal error between predicting results and experiment data

5 结论1) 利用进给系统温度检测数据和二分搜索算法, 提出了集成蒙特卡洛法和有限元法的热源发热率辨识方法, 并通过各热源发热总功率与伺服电机输出功率的对比, 验证了辨识结果的有效性.
2) 基于进给系统各热源发热率与系统运行时间和进给速度的函数关系, 推导系统工作过程各热源的分配比关系式, 进而提出了利用伺服电机的力矩电流预测各热源发热率的方法.
3) 基于丝杠一维推导方程的有限差分法, 提出了利用进给系统的各热源发热率预测丝杠温升和热误差分布的方法, 并利用试验检测结果验证预测方法的精确性.
参考文献

[1]	Ramesh R, Mannan M A, Poo A N. Error compensation in machine tools—a review.part Ⅱ:thermal errors[J].International Journal of Machine Tools & Manufacture, 2000, 40(9): 1257–1284.
[2]	Bryan J. International status of thermal error research[J].CRIP Annals Manufacturing Technology, 1990, 39(2): 645–656.DOI:10.1016/S0007-8506(07)63001-7
[3]	Yun W S, Kim S K, Dong W C. Thermal error analysis for a CNC lathe feed drive system[J].International Journal of Machine Tools & Manufacture, 1999, 39(7): 1087–1101.
[4]	Xu Z Z, Liu X J, Kim H K, et al. Thermal error forecast and performance evaluation for an air-cooling ball screw system[J].International Journal of Machine Tools & Manufacture, 2011, 51(7): 605–611.
[5]	Chen J S, Yuan J, Ni J. Thermal error modelling for real-time error compensation[J].International Journal of Advanced Manufacturing Technology, 1996, 12(4): 266–275.DOI:10.1007/BF01239613
[6]	Li Y, Zhao W, Wu W, et al. Thermal error modeling of the spindle based on multiple variables for the precision machine tool[J].International Journal of Advanced Manufacturing Technology, 2014, 72(9/10/11/12): 1415–1427.
[7]	陈诚, 裘祖荣, 李醒飞, 等. 伺服系统中滚珠丝杠的温度场模型[J].光学精密工程, 2011, 19(5): 1151–1158. ( Chen Cheng, Qiu Zu-rong, Li Xing-fei, et al. Temperature field model of ball screws used in servo systems[J].Optics and Precision Engineering, 2011, 19(5): 1151–1158.)
[8]	Takabi J, Khonsari M M. Experimental testing and thermal analysis of ball bearings[J].Tribology International, 2013, 60(7): 93–103.
[9]	Fang B, Gu T, Ye D, et al. An improved thermo-mechanical model for vertical machining center[J].International Journal of Advanced Manufacturing Technology, 2016, 87(9/10/11/12): 1–12.
[10]	刘彦丰, 高正阳, 梁秀俊. 传热学[M]. 北京: 中国电力出版社, 2015: 173-182. ( Liu Yan-feng, Gao Zheng-yang, Liang Xiu-jun. Heat transfering[M]. Beijing: China Electric Power Press, 2015: 173-182.)