基于热特性分析的几何与热复合的定位误差建模及补偿

冯文龙，南江红，杨建国

(1. 上海航天电子技术研究所，上海 201109； 2. 上海交通大学 机械与动力工程学院，上海 200240)

0 引言

为能生产高精度的工业产品，减小加工过程中的各项误差是一种最为有效的途径。影响加工精度的误差包括各部件的几何误差、热误差、切削力误差、伺服控制误差及插补误差等[1]。其中，热误差的比重达到了40%～70%[2]。对于丝杠螺母系统，丝杠受内部加工热源的影响，其温度发生较大变化，从而受热膨胀产生热误差。内部加工热源包括丝杠螺母摩擦热源、轴承热源、变速箱热源、电机热源等[3]。在热误差的影响下，机床加工误差往往会超出容许误差10倍以上，故控制机床热误差越来越受到人们重视[4]。

减小热误差的方法一般可分为误差补偿法和误差防止法[5-7]。误差防止法是通过优化机床结构设计和提高各部件加工精度的方法来提高机床的整体性能。但因存在多种物理限制，此种方法对机床精度的提升非常有限，且应用此方法消耗的成本较大[8]。误差补偿法是人为制造一项误差用于补偿，该误差与原始误差相比，大小相等，方向相反，可抵消原始误差[9]。基于现有先进的误差识别技术的支持，机床原始误差能被激光干涉仪等测量设备测得。相较于误差防止法，误差补偿不乏是一种经济有效的精度提升方案。文献[10]中提出了基于正交多项式的热误差建模方法，补偿后，在既定运动轨迹下，几何误差及其热误差减小了80%以上，但其使用单个温度点来表征整个丝杠的热状态具有先天性缺陷，一旦丝杠的运动范围改变，该热误差模型将不再适用于变化后的加工环境，误差补偿效果也将大打折扣。

本文考虑复杂的加工状况，即丝杠并非处于热均匀状态。把丝杆轴细分成多个细小的热均匀单元，并将每个细小热均匀单元应用到热均匀丝杠定位误差模型，计算可得任意时刻每个细小单元的几何与热复合的定位误差值，继而把这些细小单元的定位误差模型叠加，建立整个丝杠轴的几何与热复合的定位误差模型。通过实验验证，结果表明：补偿后，丝杠螺母传动系统的热误差明显减小，机床加工精度显著提高。

1 丝杠热特性分析

为研究丝杠的热误差，利用立式加工中心(VMC850E)展开实验。

以立式加工中心的Y轴为例,丝杠螺母传动系统机构及其传动机理如图1所示，在传动过程中，系统根据需要移动的位移量，依据其丝杠和螺母的螺距值，可精确计算得到伺服电机需要旋转的角度值。然而丝杠本身存在螺距误差，且当外界环境变化时，其螺距误差也发生变化，是一项时变误差。虽然该误差导致的丝杠转动圈数相同，但其带来的位移变化却并不相同。故对丝杠螺母系统定位误差的研究，可转化为对丝杠螺距误差的研究，即对丝杠伸长量的研究。丝杠轴正常工作时，会受到电机热源、前后轴承热源、环境温度变化以及丝杠和螺母摩擦热源的影响，从而引起丝杠本身的热温升，继而发生热伸长，导致丝杠螺距发生变化。可在电机外壳、前后轴承座、丝杠螺母以及床身等多处布置温度传感器，并实时监测这些位置的温度，来评估其对丝杠伸长的影响。

在对丝杠展开研究之前，为便于理解，作如下假设：

1) 丝杠被认为是一个实心的圆柱体，其热伸长或者冷缩现象均被认为是线性的。

2) 丝杠表面与内部(截面上)的温度一致。

3) 以空气热传导方式影响丝杠温度的热源，实验证明其影响较小，如电机热源、前后轴承热源等,均不在本文考虑范围之内。本文认为丝杠和螺母的摩擦热源是导致丝杠伸长的主要因素。

4) 在同一进给速度下，认为丝杠和螺母摩擦产生的热量是相同的。

5) 将丝杠等分成多个细小的单元，并假设这些单元是热均匀的，即在任何时候，独立单元的每个位置均处于同一温度下，且任何一个独立单元，只要其处在同一温度下，该单元的热状态就认为是相同的，与历史过程无关，即每个单元的温度值与热伸长量一一对应。

丝杠受热源的影响，会发生热胀冷缩现象，其长度变化量可表示为

ΔLscr=ΔTscr·L·α

(1)

式中：ΔLscr表示丝杠长度的变化量，单位为μm，正值表示伸长，负值表示缩短；ΔTscr表示丝杠温度的变化量，单位为℃，正值表示升高，负值表示降低；L表示丝杠的总长度，单位为m；α表示丝杠的线性热膨胀系数，即丝杠改变1 ℃长度的变化量，取值为11.7×10-3，单位为℃-1。

丝杠的热伸长量等同于丝杠的定位热误差量，使用激光干涉仪，对丝杠的定位误差进行测量，测量结果如图2所示，每个时间点的定位误差曲线形状基本保持不变，每条曲线可由0 min这条曲线绕着起始点旋转一定角度得到。故所测得的定位误差数据可拆分成两部分，一部分是固定的几何定位误差，只跟坐标位置信息相关，也即0 min这条曲线；另一部分是热复合定位误差，跟旋转角度及坐标信息相关，将旋转角度对应的斜率变化定义为此时丝杠的热状态参数。实测定位误差采用误差分离技术，表示为

δp(y,T)=δpG(y)+δpT(T)·y

(2)

式中：δp(y,T)表示跟Y轴坐标位置及温度相关的定位误差；δpG(y)为只跟Y轴位置相关的几何定位误差；δpT(T)表示热状态参数，即跟旋转角相关的斜率变化值；y表示对应的Y轴坐标值。

图2 Y轴几何与热复合的定位误差Fig.2 Geometrical error and thermally induced positioning error of Y-axis

丝杠定位热误差等同于丝杠的热伸长量，对比式(1)和式(2)，可得到热状态参数的表达式

δpT(T)=ΔTscr·α

(3)

由式(3)可知，热状态参数δpT(T)跟丝杠的温升量呈线性关系，只要能准确地测量出丝杠的实时温度，即可求解热状态参数，继而获得定位热误差，完成对其的建模过程。然而丝杠在工作条件下将一直处于旋转状态，无法在其表面直接吸附传感器用以测量其温度值，可通过测量丝杠螺母的实时温度，来近似得到丝杠的温度。在丝杠螺母摩擦的过程中，丝杠会受丝杠螺母摩擦热源的影响，温度升高，同时，丝杠也会向低于其温度的外界空气释放出热量，温度降低。丝杠在温升阶段的热状态可由以下瞬态热平衡方程表示

(4)

式中：Qnut为丝杠和螺母摩擦产生的热量；hw为强制对流热传导系数；A为热扩散面积；Tscr、Tamb分别为丝杠温度和环境温度；ρ为丝杠材料密度；c为丝杠材料比热容；V为丝杠体积；tup为丝杠温升的累积时间。

为得到丝杠瞬态温度Tscr的表达式，对式(4)作解微分方程运算。丝杠瞬态温度可表示为

(5)

而当丝杠停止转动，丝杠和螺母停止摩擦，机床处于停止状态时，丝杠不再受摩擦热源的影响，而只是向低于其温度的外界空气散热，此时丝杠处于瞬态热平衡，热平衡方程为

(6)

式中：hc表示自由对流热传导系数；tdown表示丝杠降温累积时间。

对式(6)作解微分方程运算，得到降温状态时丝杠瞬态温度表达式，即

(7)

丝杠被认为是一个实心的圆柱体，故其体积V和散热面积A都可用丝杠直径d来表示，综合式(2)、(3)、(5)和式(7)，可得到温升和冷却阶段丝杠定位误差求解表达式：

(8)

由式(8)可知，丝杠的定位误差与丝杠和螺母的摩擦热量、丝杠直径、散热面积、强制对流系数、自由对流热传导系数、温升时间和冷却时间等参数相关。其中散热面积和运动范围相关，丝杠和螺母的摩擦热源与进给速度相关。当加工参数发生变化时，丝杠温度预测方程式也可通过调整模型系数，来适应环境的变化，以此保证丝杠温度函数对丝杠温度的准确预测。摩擦热求解表达式为

(9)

式中：Qnut为丝杠螺母摩擦产生的发热量，单位为W；M为丝杠螺母的摩擦力矩，单位为N·mm；n为丝杠的转速，单位为r/min；a为系数；F为机床进给速度，单位为mm/min；S为丝杠的导程，单位为mm。

2 丝杠定位误差建模

2.1 丝杠定位误差模型理论

分析丝杠在工作环境中的热状态可得丝杠在温升和冷却阶段的定位误差模型，但其表达式中仍有很多未知的参数，不能直接用于定位误差的计算，本节将求解这些未知系数。

以2 000 mm/min的进给速度来回移动Y轴工作台，移动范围为整个行程。丝杠温度会伴随着移动而上升，在Y轴螺母处布置温度传感器，用以测量和螺母接触处丝杠的温度。丝杠和螺母摩擦共同产生摩擦热，且因丝杠和螺母紧密连接，故可认为在丝杠和螺母接触处丝杠和螺母的温度是相同的。与此同时，在靠近丝杠螺母传动系统后轴承的床身处布置温度传感器，监控丝杠温升量基准，因整个实验过程耗时较长，很难保证整个过程中环境温度不发生变化。而丝杠螺母传动系统中，后轴承处受到丝杠螺母摩擦热源影响较小，其温度仅受环境温度变化的影响，能有效测量出丝杠在环境温度影响中的温度变化，螺母处温度传感器值和此处温度的差值即丝杠仅受丝杠螺母摩擦热源影响时的温度变化量。移动Y轴工作台直到丝杠到达热平衡状态，继而停止工作台，进入丝杠冷却阶段，在这个过程中，实时测量并记录螺母温度和丝杠基准温度。螺母位置和轴承端位置的温度传感器布置情况如图3所示。

图3 温度传感器安装图Fig.3 Temperature sensors installation on machine tool

图4 温升和冷却阶段丝杠温度变化量Fig.4 Temperature variations during warming and cooling phases

指数拟合测得的温度曲线，可得丝杠温度的预测模型。测量得到的丝杠温度数据拟合情况如图4所示。拟合得到的丝杠温度的预测模型可表示为

(10)

由式(10)可知，冷却阶段的指数项系数与温升阶段不同，究其原因：理论模型的前提是温升到热平衡后再降温，其降温的起点是热平衡状态，但在实际实验中，温升阶段历时140 min，看起来曲线变化趋势已经逼近热平衡状态，但其终究还是没有达到热平衡，若在此点降温，其相距热平衡状态必有一个时间间隔。对比式(10)和式(5)、(7)，可得

(11)

式中：丝杠密度ρ=7.85×103kg/m3；丝杠比热容c=0.46×103J/(kg·℃)；丝杠直径d=40 mm；丝杠长度L=850 mm；丝杠散热面积A=0.11 m2。

根据上述已知参数，求解式(11)，可以得到：hw=13.4 W/(m2·℃),Qnut=16.5 W,hc=12.0 W/(m2·℃)。

把求得的参数代入式(8)，其中线性热膨胀系数取值α=11.7×10-3μm/mm，可得丝杠定位误差预测模型，其表达式为

δp(y,T)=

(12)

在实际加工过程中，在温升阶段，依据定位误差模型，实时定位误差可易得；在热平衡后冷却阶段，实时定位误差也可轻松获得。而在实际的加工中，很少出现理想的温升到热平衡再冷却的情况，一般都是温升一段时间就开始冷却，冷却一段时间又开始温升，温升和冷却阶段随机交替出现。故下文将研究上述温升和冷却阶段随机交替出现时定位误差的计算方法。

如图5所示，当丝杠从冷态开始移动，一直处于温升状态，热状态参数也随着丝杠温度的上升而变大，但在A点，丝杠停止运动，开始冷却，此时的丝杠可看作一个均匀的整体，且丝杠各点温度基本相同，找到与其热状态相同的一个点，且处于自由冷却中，复制其变化趋势就能得到丝杠中间过程冷却的热状态参数变化趋势。作一条平行于时间轴的直线，交冷却阶段曲线于B点，平移曲线BD至A、B点重合即可得到曲线AC，丝杠中间过程冷却后的热状态参数变化趋势如曲线AC所示。

图5 中间过程冷却情况下的热状态参数变化趋势Fig.5 Thermal parameter variation when cooling in the midway of warming up

丝杠在温升过程中，中途进行冷却后，丝杠处于冷却过程中，此时，丝杠重新开始移动，丝杠温度必定又会上升，丝杠的热状态参数变化趋势如图6所示，丝杠在经历第1次温升段OA后，在A点冷却，冷却至F点后，再度温升，如曲线段FG，直至热平衡。过F点作时间轴的平行线，交第1次温升曲线OH于E点，平移曲线EH直至E点和F点重合，即可获得2次温升曲线段FG。在丝杠温升阶段FG过程中，如果再次出现冷却和冷却后再温升的情况，可通过类似的方法计算丝杠的热状态参数，即当丝杠冷却阶段和温升阶段随机交替时丝杠热状态参数都可通过上述计算方法计算得到。

图6 丝杠二次温升对应的热状态参数变化趋势Fig.6 Thermal parameter variation in double warming situation

2.2 分段指数模型理论

2.1节建立了几何与热复合的定位误差模型，兼容温升及冷却交替出现情况。该模型利用丝杠全行程来回移动模拟加工运动，把丝杠等效为一个均匀的整体，继而研究丝杠的热状态。但在实际的加工中，基本不存在全行程来回移动的情况，实际的加工状况是围绕工件，在全行程的某个区间内来回移动。实际的加工条件和2.1节中建立的误差模型前提不相符，这也预示着建立的误差模型不能直接应用于实际加工的误差补偿中。

在实际的加工中，加工工件不同，运行的G代码程序不同。机床的运动轨迹是无规律及不可预测的，这势必带来一个严重问题，即丝杠各点受到摩擦热源的影响将不再相同，丝杠各点的热状态将不再相等。使用一个热状态参数已经不能表征整根丝杠的热状态。本文采用对丝杠进行细分的方法，把丝杠等分成若干足够小的单元，因其足够小，其内部就不存在温度梯度，近似可认为其绝对热均匀，对每个等分单元应用2.1节建立的误差模型，得到每个等分单元的几何与热复合的定位误差模型，再将这些等分单元的误差模型叠加，即可得到整根丝杠的误差模型。

为验证这一方法的可行性，实验研究机床的Y向运动。为模拟实际工件的加工状况，使丝杠各点处于不同的热状态。本实验采用各区域依次热机的方法，热机轨迹为：

1) 螺母工作台在位置坐标为200 mm和600 mm之间以进给速度2 000 mm/min来回移动10 min；

2) 螺母工作台再在位置坐标为0 mm和400 mm之间以进给速度2 000 mm/min来回移动10 min；

3) 螺母工作台在位置坐标为400 mm和800 mm之间以进给速度2 000 mm/min来回移动10 min。

运动过程中，实时记录机床的位置坐标值，如图7所示。丝杠在这30 min温升阶段的热状态参数模型将不同于上述全行程运动温升过程，因其运动范围只是全行程的一半，故其散热面积是全行程温升情况的一半，用A2=A/2表示。上述30 min温升过程中，丝杠温升和冷却时热状态参数计算方法可表示为

图7 温升阶段工作台位置坐标Fig.7 Recorded location of nut during warming phase

经过实验验证，选定最小的热均匀单元长度为40 mm，将Y轴等分为20份，对这20个细小热均匀单元分别建模。建模结果如表1所示。

表1 20个单元的热状态参数值

20个独立单元的热状态参数分别用δpT1,δpT2,δpT3,…,δpT20表示，整根丝杠的定位误差可表示为

(14)

式中：δpG(y)表示Y轴的几何定位误差；p为丝杠的位置，单位为mm。

几何定位误差可通过最小二乘法建模，建模结果可表示为

几何定位误差和热状态参数建模结束后，依据式(14)，即可预测丝杠任意时刻的定位误差。为评估定位误差模型的模型精度，在上述30 min热机过程中，每隔10 min测量一次丝杠的定位误差，并对比测量误差值和模型预测误差值之间的差异，用以验证模型的准确性。对比结果如图8所示，结果表明：上文所建立的定位误差模型能准确预测上述30 min实验过程中的定位误差，该误差在5 μm以内，预测精度大于80%。

此定位误差建模方法对运动轨迹没有诸多限制，其适用于任何加工轨迹后的定位误差预测，可广泛应用于实际生产加工过程中，只需要同时对每个独立细分单元做热状态分析，准确记录每个单元的温升和冷却时间，并依据模型算法计算每个单元各个时刻的热状态，根据式(14)叠加，即可获得丝杠定位误差模型。用户可根据自己的要求来细分丝杠，细分单元越小，模型精度越高，定位误差建模流程图如图9所示。

图8 定位误差模型预测值与实测值对比图Fig.8 Comparison of positioning errors between the measured values and the predicted values

图9 定位误差建模流程图Fig.9 Algorithm flow chart of positioning error mode

3 补偿结果分析

为验证模型的鲁棒性和有效性，在立式加工中心VMC850E上进行一系列验证实验。首先，机床空载运行，以2 000 mm/min的进给速度在位置为300 mm和500 mm之间来回移动30 min，分别测量机床Y轴补偿前和补偿后的定位误差；继而，机床停机10 min，自由冷却，然后再分别测量机床Y轴补偿前和补偿后的定位误差，对比这两组补偿前后的定位误差量，评估补偿的效果。对比结果如图10所示，在300～500 mm段定位误差在温升后有明显增加，冷却后又有明显下降。经过指数模型补偿后，定位误差从30 μm减小至8 μm，然而多项式模型补偿后残差更大，这是因为用单个温度传感器来表征整个丝杠的热状态对复杂的运动情况已不再适用。

图10 补偿结果对比Fig.10 Comparison of compensation results

接下来进行加工实验，对6个试件进行定位孔加工，前2个试件分别在有补偿和无补偿的条件下进行，因从冷态开始加工，这2个试件的加工条件认为是初始状态；再加工2个试件，这段过程被认为是温升过程；最后加工2个试件，分别在有补偿和无补偿的条件下进行，这2个试件的加工条件认为是温升后状态。加工过程和试件如图11所示。

对加工完成试件的4个孔距送交第三方检测机构“上海市计量测试技术研究院”检测，测试结果如表2所示，孔距从最大10 μm减小为3 μm，补偿效果显著。

图11 试件加工过程Fig.11 Machining process of test pieces

mm

4 结束语

针对旋转编码器定位的丝杠螺母传动系统，在温升和冷却阶段分析丝杠轴的热特性，建立丝杠轴瞬态温度和时间的指数关系表达式，根据温升或冷却时间计算丝杠实时瞬态温度。基于热膨胀理论，根据丝杠轴温度，建立热均匀丝杠几何与热复合的定位误差模型。考虑到复杂的加工状况，丝杠并非处于热均匀状态，把丝杆轴细分成多个细小热均匀单元，对每个细小热均匀单元应用热均匀丝杠定位误差模型，可计算得到任意时刻每个细小单元的几何与热复合的定位误差值，继而把这些细小单元的定位误差模型叠加，即可建立整个丝杠轴的几何与热复合的定位误差模型。应用该误差模型对多个试件进行加工，经过补偿后，试件的孔距精度明显提高。研究结果表明：应用该误差补偿方法，可有效地减小加工热源影响时批量生产零件的尺寸偏差，相比于基于单温度点的误差补偿方法，该补偿方法更有效。不仅如此，因其是对丝杠螺母传动系统进行建模，故其对机床结构约束较小，可适用于大部分现有机床结构。针对大批量的零件加工任务，可实现各加工机床之间组网，并上载显示各机床的加工状态，实现对各机床加工状态的实时监控，此项工作可作为本文的后续研究方向。

[1] LEE D S, CHOI J Y, CHOI D H. ICA based thermal source extraction and thermal distortion compensation method for a machine tool[J]. International Journal of Machine Tools and Manufacture, 2003, 43(6): 589-597.

[2] RAMESH R, MANNAN M A, POO A N. Error compensation in machine tools: a review Part II: thermal errors[J]. International Journal of Machine Tools and Manufacture, 2000, 40(9): 1257-1284.

[3] XUZ 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 and Manufacture, 2011, 51 (7/8): 605-611.

[4] RAMESH R, MANNAN M A, POO A N. Support vector machines model for classification of thermal error in machine tools[J]. The International Journal of Advanced Manufacturing Technology , 2002, 20(2): 114-120.

[5] NI J. CNC machine accuracy enhancement through real-time error compensation[J]. Journal of Manufacturing Science and Engineering, 1997, 119(4B): 717-725.

[6] 杨建国. 数控机床误差综合补偿技术及应用[D]. 上海: 上海交通大学, 1998.

[7] 冯文龙. 大型数控机床多误差元素建模及综合补偿[D]. 上海: 上海交通大学, 2016.

[8] WANG W, ZHANG Y, YANG J G, et al. Geometric and thermal error compensation for CNC milling machines based on Newton interpolation method[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2013, 227(4): 771-778.

[9] HAN Z Y, JIN H Y, LIU Y L, et al. A Review of Geometric Error Modeling and Error Detection for CNC Machine Tool[J]. Applied Mechanics and Materials, 2013, 303/304/305/306: 627-631.

[10] FAN K G, YANG J G, YANG L Y. Orthogonal polynomials-based thermally induced spindle and geometric error modeling and compensation[J]. The International Journal of Advanced Manufacturing Technology, 2013, 65(9/10/11/12): 1791-1800.