摘 要：为了降低加工过程的热误差，提高数控机床加工精度，基于时序相关分析理论与数值计算方法，建立了一种以温度场分布及加工参数为输入的新型机床主轴热误差建模方法.所建模型由热误差模型、主轴动压轴承热特性模型以及主轴热传递模型三部分组成.该方法首先根据时序相关理论建立热误差与温度测点之间的相关模型，再通过灰色相关理论完成关键温度测点位置与数量的优化，同时，基于数值计算与热传导理论，建立了动压主轴系统热特性模型.以一台大型龙门导轨磨床为实验对象，建立了磨床主轴箱热误差预测模型.实验结果表明，所建立的热误差模型具有良好的热误差辨识性能.
关键词：时序相关分析；雷诺方程；有限差分法； 热误差
中图分类号：TH161 文献标识码：A
机床热误差是精密加工误差的主要来源之一，对其控制和补偿是提高机床加工精度的关键技术＼[1＼].机床热误差补偿涉及热误差检测、热误差预测模型和热误差实时补偿三类问题，其中热误差预测模型的建立最为关键＼[2＼].
机床热误差控制的最新技术是基于热误差预测模型的误差实时补偿技术（RTEC）＼[3＼].由于机床热误差可以看作是关于机床热分布、位置以及运动方向的确定性函数＼[4＼]，故而，可以建立以机床温度场分布和运动方向为自变量的热误差模型.同时，机床温度场分布特性可以用若干温度测点处的数据来表征＼[5＼].
本文运用时间序列分析理论拟合热误差与相关温度测点数据之间的关系，采用数值计算方法建立某高精度磨床主轴热特性模型，提高了此类热误差模型的预测精度和鲁棒性.
1 机床热误差模型的理论基础
1.1 时序分析建模方法
时序分析是对观测所得的有序随机数据进行分析和处理的一种方法，适用于输入不可测或不确定测量数据的描述，常被用于机床热误差建模.
根据时序分析理论＼[6＼]，对一组测点进行第i次测量的热误差χi可以分解为两部分，如式（1）和式（2）所示，一部分是确定性部分，完全取决于各温度测点处温度序列Tr，i的线性组合fi；另一部分是随机性部分，完全独立于Tr，i，由残量εi确定（如仪器的测量误差）.同时，假定εi是零均值的独立序列，即白噪声序列，因而，对于不同i值，εi是相互独立的.
χi=fi+εi.（1）
fi=β0+β1T1，i+β2T2，i+…+βRTR，i.（2）
式中：i=1，2，…，I，βo为常数，βi为系数，r为温度测点序号，R为温度测点总数.
基于热误差χi和各温度测点处温度Tr，i的实测样本序列，采用最小二乘法获取系数β0和βi的估计值0和i，如式（3）所示：
=TTT-1TTχ.（3）
式中，χ=χ1，χ2，…，χIT，
=o，1，2，…，RT，
T=1T1，1T2，1…TR，11T1，2T2，2…TR，21T1，IT2，I…TR，I.
将系数估计值代入式（2）和（1），得到残量序列εi，如式（4）所示：
εi=χi-o+1T1，i+2T2，i，…，RTR，i.（4）
采用一阶自回归模型描述残差，如式（5）所示：
εi=φ1εi-1+ai.（5）
采用最小二乘估计求解上述模型的参数，如式（6）所示.同时，ai是均值为0的正态分布序列：
1=∑Ii=2εiεi-1/∑Ii=2ε2i.（6）
综上所述，由一组测点确定的热误差计算模型如式（7）所示：
χi=0+1T1，i+2T2，i+…+RTR，i+10+1T1，i-1+2T2，i-1+…+RTR，i-1+ai.（7）
1.2 机床温度测点优化方法
机床各部分温度的分布特性是热误差产生的内在原因，因而，温度测点数量和位置的选取将直接影响式（7）计算误差的准确性.基于已有的测试数据，采用灰色关联度分析方法进行温度测点位置和数量的优化选取，可以将测点数目控制在合理范围内的同时，提高式（7）的计算精度.
灰色关联度分析是按照系统中各特征参量序列之间的相似程度进行系统分析的一种方法[7-8]，保证各特征参量间的等效性和同序性，是进行关联度分析的前提.
采用极差变换法的一种改进形式，对热误差与各温度测点的实测结果进行无量纲化处理，如式（8）所示：
x（k）=max x（k）-x（k）max x（k）-min x（k），
k=1，2，…，n.（8）
设对应于热误差实测值的x0为参考序列，对应于各温度测点的温度实测值的xi为比较序列，则x0对xi在第k点的灰关联系数如式（9）所示：
ξ0i（k）=min imin kΔ0i（k）+ρmax imax kΔ0i（k）Δ0i（k）+ρmax imax kΔ0i（k）.（9）
其中，Δ0i（k）为第k点x0与xi差值的绝对值，ρ为分辨系数，min imin kΔ0i（k）为Δ0i（k）的两级最小值，max imax kΔ0i（k）为Δ0i（k）的两级最大值.两序列之间的灰关联度为：
γ0i=1m∑mk=1ξ0i（k）.（10）
式中：γ0i为比较序列xi与参考序列x0的关联度，m为比较序列的数据个数.
将各个比较序列对同一参考序列的关联度按大小顺序排列起来，即组成了反映各个比较序列对同一参考序列的“主次”关系的关联序.依据关联序即可对磨床主轴各温度测点位置分布及测点数量进行优化处理.
1.3 动压轴承的热特性建模方法
精密磨床运行时，主轴轴承的摩擦热是热误差的主要来源，通常也是主要的温度测点.在机床的实际运行中，高速旋转和结构复杂的特点，难以用传感器对轴承处的温度进行准确测量.通过建立主轴的热特性描述模型，计算不同工况下主轴轴承处的温度参数用于热误差建模，对于提高模型的精度和适用性都有主要意义.
建立主轴动压滑动轴承油膜热特性计算模型的思路是通过对雷诺方程、温粘方程等的联立求解，得到动压轴承达到平衡状态时的平均温度.计算流程如图1所示.
由于刚性径向滑动轴承的油膜压强分布满足雷诺方程，基于相似理论，将坐标原点置于轴承轴线上，且位于轴承宽度的1/2处，对雷诺方程进行无量纲化后得到＼[9＼]：
φ3φ+DB2γh-3γ=6φ.（11）
其中， =ψ2p/6μVtB，ψ为相对间隙，μ为润滑油动力粘度，p为油膜压强，Vt为轴颈上任意点的切向速度，φ=x/R，x为轴颈上任意点的周向坐标，R为轴瓦半径，为油膜无量纲厚度，=1+εcos φ，ε为偏心率，γ=2y/B，y为轴向坐标，B为轴承宽度，D为轴瓦直径.
式（11）基于如下假设：
1）润滑油重量及惯性力忽略不计；
2）油膜为牛顿流体，且为层流；
3）润滑油不可压缩；
4）润滑油为定常流体；
5）载荷方向为垂直于轴向.
对于任意动压轴承油膜，若求解区域为φ1≤φ≤φ2，0≤γ≤1，则将求解区域网格化，且令i表示沿φ方向的列数，j表示沿γ方向的行数.根据有限差分原理，将式（11）化为有限差分方程求解.
1.3.1 雷诺方程的有限差分化
采用半步长插入点上压强的中差商来表示任意节点上压强的一阶导数，相邻半步长插入点上压强的一阶导数的中差商来表示压强的二阶导数，则式（11）可表示为：
求解上述方程时，先确定承载区四周边界的压强，进而求解整个求解区域上的压强分布.
1.3.2 油膜偏位角的修正
由于已知载荷方向，为了确保油膜压力分布求解的准确性，需要进行油膜偏位角的校核，从而保证油膜承载力的合力与载荷方向相反.
已知载荷方向为垂直于轴向，则轴承的偏位角应保证油与膜合力方向相反，采用迭代法求解稳定状态下轴承的偏位角.修正方程为：
1.3.3 温粘方程的求解
油膜温升的求解需要考虑温度对润滑油粘度的影响.为了克服采用能量方程求解粘度过程的复杂且耗时的缺点，采用润滑油标准粘度η50（t=50℃）来计算任意润滑油温度下的粘度＼[10＼]：
η=ηxexp 400lg η50ηx1t+95-0.001.（17）
其中，t为油膜温度，ηx=0.18 mPa·s.
将修正后的粘度值代入雷诺方程进行迭代，直到收敛，此时，所得压强分布即为精确解.
1.3.4 油膜温升的计算
已知轴颈表面的摩擦阻力为：
Ff=vωR2Bψ∑ζa=1∫10∫φa2φa11+311φdφd.（18）
轴颈表面的摩擦系数为f=f/，其中， =x2+z2.
基于流量连续原理和差分原理，轴承泄油量为：
Q=∑ζa=12R3ωψBΔφΔγ∑m-103i，n-1-i，n.（19）
基于热平衡公式，假设由摩擦功转化的热量全部由润滑油及热扩散带走，则热平衡方程式为＼[11＼]：
Δt=f·FfψDB/cρQψUDB+πksψU.（20）
则油膜温度为：
t=t0+f·FfψDB/cρQψUDB+πksψU.（21）
其中，t0为润滑油初始温度； c为润滑油比热容；ρ为润滑油密度；Q为轴承泄油量，ks为散热系数.
1.4 主轴热传递模型
已知主轴运转过程中，动压轴承油膜温度如式（21）所示，下面讨论求解主轴热特性分布的另一个关键条件，即主轴热传递模型.
采用圆柱坐标系描述主轴温度分布，则不含内热源的非稳态导热微分方程为＼[12＼]：
2Tr2+1rTr+2Tz2=1αTt.（22）
其中，T为温度，r为半径，z为轴向坐标，t为时间，α为主轴散热系数.
式（22）还基于如下假设：
1）主轴材质均匀且各向同性；
2）主轴材料的特性系数为常数；
3）物体内无温度跃变或物体相变；
4）主轴温度不随角度φ的变化而变化.
采用分离变量法，假定式（22）有如下变量分离形式的解：
2 热误差模型精度分析
下面以某高精度龙门导轨磨床主轴箱系统（如图2所示）为研究对象，依据上述研究结果进行测试和建模，并对热误差模型精度进行分析验证.
首先根据经验和初步分析确定12个测试点：主轴后端皮带轮附近布置2个测点（1～2号测点）；主轴箱体上布置4个测点（3～6号测点）；主轴前段布置3个测点（7～9号测点）；电机上布置1个测点（10号测点）；安装座上布置2个测点（11～12号测点），如图2所示.进行温度测试采样.测试条件为，磨头主轴以1 040 r/min和1 500 r/min的恒定转速各空载运行2.5 h，每间隔5 min用非接触式红外温度测量仪测量12个测点的温度各一次，与此同时，用千分表测量主轴在X，Y，Z三个方向上对应的热漂移，总共获取2组温度和热误差数据.
2.1 主轴箱系统温度测点优化
采用灰色关联度分析法，将所测得的主轴在X，Y，Z三个方向上的热误差作为参考序列，所测得的12组温度测点数据作为比较序列.利用式（8）对各数据序列进行处理，通过式（9）和式（10）求得各比较序列与参考序列间的灰关联度，如图3所示.
经过关联度大小的对比，排除关联度小的温度测点后，确定影响主轴三个方向热漂移的关键温度测点总数由12个降低为8个，分别是：1号、3号、4号、5号、6号、7号、11号及12号等8个温度测点.其中，1号与7号测点位于主轴轴承处，需采用2.3节与2.4节建立的热特性模型对该两个测点处的温度进行预测.
2.2 主轴热特性建模与精度分析
实验中所采用的磨床的主轴箱动压轴承及其润滑油的参数包括：宽径比0.8，初始偏位角90°，偏心率0.5，轴颈半径0.120 m，轴瓦半径0.167 m，内切圆半径0.121 m，轴承张角147°，初始润滑油动力粘度0.02 Pa·s，初始润滑油温度28 ℃，润滑油比热1 888 J/（kg·℃），润滑油密度867 kg/m3，轴承散热系数80 W/m2·℃，主轴转速为1 470 r/min或1 500 r/min.
采用主轴热特性模型和主轴热传递模型，对1号与7号测点处温度时间序列进行预测，同时，将预测值序列i与两个测点处的温度实测值序列xi进行对比，并计算其相对误差δ，结果如表1所示.相对误差的计算公式如式（27）所示：
δ=xi-ixi.（27）
由表1可知：预测值序列与实测值序列的相对误差较小，相对误差平均值为14.2%.其中，1号测点相对误差最大值13.3%，7号测点相对误差36.4%.
因此，主轴热特性辨识模型的预测数据与实测数据基本吻合，另外考虑到1号和7号测点位于主轴外露的轴颈处，测试所得温度与轴承处的温度本身存在一定的误差，所以，本文建立的模型能有效地描述主轴的热特性.
2.3 热误差模型实验验证
综合主轴热特性辨识数据和温度测点数据，利用式（7）进行热误差模型的验证.
首先，基于第1组数据中各温度测点数据序列和热误差数据序列，建立主轴X，Y，Z三个方向上的热误差预测模型.然后，选取未使用过的第2组数据中的温度数据，代入热误差预测模型中，输出X，Y，Z方向的热误差预测值；最后，将预测值与热误差实测值比较.同时，预测过程中，第1号和7号测点的温度数据采用主轴箱热特性模型的输出值.模型的热误差预测值与实测值，以及两者之间的残差如图4，图5和图6所示.
3 结 论
1）基于空载实验的精密导轨磨床热误差建模方法，研究的热误差不受其它加工因素的耦合影响，降低了热误差建模的复杂程度，对于切削力不大的精密磨削具有实用意义.
2）基于时序相关分析理论与数值计算理论的热误差模型能有效地对机床热误差进行辨识.
3）实验验证了该建模方法的有效性与可行性，为机床精密加工精度的提高与热误差实时补偿技术提供了理论借鉴.
参考文献
[1] 林伟青， 傅建中， 许亚洲， 等. 基于最小二乘向量机的数控机床热误差预测[J]. 浙江大学学报， 2008， 42（6）： 905-908.
LIN Weiqing， FU Jianzhong， XU Yazhou， et al. Thermal error prediction of numerical control machine tools based on least squares support vector machines [J]. Journal of Zhejiang University， 2008， 42（6）： 905-908.（In Chinese）
[2] 李永祥， 杨建国.灰色系统模型在机床热误差建模中的应用[J]. 中国机械工程， 2006， 17（23）： 2439-2442.
LI Yongxiang， YANG Jianguo. Application of grey system model to thermal error modeling on machine tools [J]. China Mechanical Engineering， 2006， 17（23）： 2439-2442.（In Chinese）
[3] YANG Hong， NI Jun. Dynamic neural network modeling for nonlinear， nonstationary machine tool thermally induced error[J]. International Journal of Machine Tool& Manufacture， 2005， 45（4）： 455-465.
[4] JOSEF M， JERZY J， ECKART U， et al. Thermal issues in machine tools [J]. Annals of CIRP， 2012， 61（2）： 771-791.
[5] WECKA M， MCKEOWNB P， BONSEC R， et al. Reduction and compensation of thermal errors in machine tools [J]. Annals of CIRP， 1995， 44（2）： 589-598.
[6] 杨叔子， 吴雅， 轩建平， 等.时间序列分析的工程应用[M]. 武汉：华中科技大学出版社， 2007.
YANG Shuzi， WU Ya， XUAN Jianping， et al. Engineering application of time series analysis[M]. Wuhan： Huazhong University of Science and Technology Press， 2007.（In Chinese）
[7] 刘思峰， 党耀国， 方志耕.灰色系统理论及其应用[M]. 北京：科学出版社， 2005.
LIU Sifeng， DANG Yaoguo， FANG Zhigeng. Theory and application of grey system[M]. Beijing： Science Press， 2005.（In Chinese）
[8] 闫嘉钰， 张宏韬， 刘国良，等. 基于灰色综合关联度的数控机床温度测点分组优化[J]. 湖南大学学报： 自然科学版，2008， 35（4）： 37-41.
YAN Jiayu， ZHANG Hongtao， LIU Guoliang， et al. Optimization of measuring points for machine tool thermal error modeling based on grouping of synthetic grey correlation method [J]. Journal of Hunan University： Natural Science， 2008， 35（4）： 17-41.（In Chinese）
[9] 杨沛然. 流体润滑数值分析[M]. 北京：国防工业出版社， 2005.
YANG Peiran. Numerical analysis of fluid lubrication[M]. Beijing： National Defense Industry Press， 2005.（In Chinese）
[10]李中华. 数控机床主轴轴承部件温度及热力学建模[J]. 机械制造， 2010， 48（6）： 30-33.
LI Zhonghua. Temperature of spindle bearings of CNC machine tools and thermodynamic modeling[J]. Mechanical Manufacturing， 2010， 48（6）： 30-33.（In Chinese）
[11]孙军， 桂长林. 内燃机曲轴强度研究的现状、讨论与展望[J]. 内燃机学报， 2002， 20（2）： 179-184.
SUN Jun， GUI Changlin. A review of crankshaft strength analysis for internal combustion engines [J]. Journal of Combustion Engine， 2002， 20（2）： 179-184.（In Chinese）
[12]俞昌铭. 热传导[M]. 北京：高等教育出版社， 1983.
YU Changming. Heat conduction [M]. Beijing： Higher Education Press， 1983.（In Chinese）