基于ANSYS Workbench的电主轴热误差预测

Thermal Error Prediction of Motorized Spindle Based on ANSYS Workbench

摘要: 针对电主轴受热变形影响加工精度的问题，本文先运用Solidworks软件对电主轴进行几何建模，然后通过有限元分析软件ANSYS Workbench对电主轴进行瞬态热–结构耦合分析，从而预测出电主轴热误差随时间变化的趋势，最后以实验室的TC-E650型数控雕铣机为研究对象验证仿真结果。实验结果表明：电主轴在运行的前1500 s内热误差变化幅度大，短时间内迅速增长到40 μm，而后缓慢增长，在运行约4500 s时达到热稳态，主轴末端最大轴向热误差为52 μm，实验结果有效验证了电主轴热误差的有限元瞬态预测模型的精度，为电主轴的热平衡设计及热误差的补偿提供了参考。

Abstract: In order to solve the problem that the thermal deformation of motorized spindle affects the machining accuracy, this paper first uses Solidworks software for geometric modeling of motorized spindle, and then carries out transient thermal-structure coupling analysis of motorized spindle by ANSYS Workbench finite element analysis software, so as to predict the trend of thermal error of motorized spindle with time change. Finally, the tC-E650 CNC engraving and milling machine in the laboratory is taken as the research object to verify the simulation results. The experimental results show that: The thermal error of the motorized spindle changes greatly in the first 1500 s of operation, rapidly increases to 40 μm in a short time, then slowly increases, and reaches thermal steady state at about 4500 s. The maximum axial thermal error of the motorized spindle is 52 μm. The experimental results effectively verify the accuracy of the finite element transient prediction model of thermal error of the motorized spindle. It provides reference for thermal balance design and thermal error compensation of motorized spindle.

1. 引言

随着现代制造技术的发展，对机床加工精度的要求越来越高，改善机床的热变形成为提高加工精度需要解决的迫切问题。热误差占机床总误差的40%~70% [1]，电主轴作为机床的核心部件，在高速运动中产生的热变形严重影响了加工精度，因此电主轴的热变形预测及补偿研究至关重要。

热误差建模预测的主要方法有经验热误差建模和理论热误差建模，根据是否考虑热弹性效应，经验热误差模型又分为准静态模型和动态模型 [2]。理论热误差建模的方法主要包括有限元法、热网络法和集中质量法，计算机科学的快速发展使有限元法成为工程分析中广泛应用的数值分析方法。因此，通过Ansys仿真分析软件可以对电主轴的整体热态特性进行深入研究。邓小雷利用有限元软件建立了机床主轴–立柱系统耦合分析模型，对机床的热态特性进行研究 [3]。张丽秀将有限元模型与实验数据相结合，建立了电主轴的温升预测模型 [4]。Andreas Naumann提出了一种多速率高阶时间步进法，对机床的有限元模型进行实时精确的瞬态热变形仿真 [5]。Liu Zhe等人利用Ansys Workbench建立了主轴的有限元模型，并对不同转速下的温度、热变形及旋转精度进行了研究 [6]。刘俊龙提出了一种基于最小二乘支持向量回归理论和遗传算法的有限元模型热学边界条件辨识方法，实现电主轴有限元模型的热学参数辨识计算 [7]。李郝林基于机床热误差有限元分析方法与信息论，对机床温度测量点的位置进行优化 [8]。范开国采用有限元方法、理论方法和实验方法，建立在线监测应用程序对机床热变形进行实时预测 [9]。

上述论文基于Ansys Workbench从瞬态时间步的设置、热边界条件的修正及温度测点的优化布置等角度对机床电主轴的热误差进行建模预测，本文根据前人的各项研究成果，以实验室的TC-E650型数控雕铣机的电主轴为研究对象，结合了传热学及有限元的相关理论对电主轴的各项热学参数和边界条件进行了分析计算，然后借助数值模拟技术对电主轴的热误差进行理论建模，从而预测出电主轴从冷态到稳态整个运行过程中的热误差随时间变化的趋势，为后续热误差的补偿研究奠定基础。

2. 电主轴的生热、传热计算

2.1. 角接触球轴承的生热计算

轴承在高速旋转时，滚珠与内外圈滚道之间摩擦生热，生热量主要与转速和摩擦力矩有关，可根据下式计算：

$Q=1.047\times {10}^{-4}Mn$ (1)

式中： $Q$ ——轴承的生热量，W； $M$ ——轴承的总摩擦力矩， $\text{N}\cdot \text{mm}$ ； $n$ ——轴承转速，r/min。

Palmgren总结出轴承摩擦力矩计算的经验公式，把轴承摩擦力矩分为外加载荷引起的摩擦力矩和润滑剂黏性引起的摩擦力矩两大部分 [10]，计算公式如下：

$\{\begin{array}{l}M=M_0+M_1\\ M_0={10}^{-7}{f}_0{\left(vn\right)}^{2/3}{d}_m^3\left(vn\ge 2000\right)\\ M_0=160\times {10}^{-7}{f}_0{d}_m^3\left(vn<2000\right)\\ M_1=f_1{P}_1{d}_m\end{array}$ (2)

式中： $M_0$ ——与轴承类型有关的摩擦力矩， $\text{N}\cdot \text{mm}$ ； $M_1$ ——与负载有关的摩擦力矩， $\text{N}\cdot \text{mm}$ ； $f_1$ ——与轴承类型相关的参数； $P_1$ ——轴承摩擦力矩负荷， $\text{N}$ ； $\nu$ ——工作温度下润滑剂的运动粘度， ${\text{mm}}^2\cdot {\text{s}}^{-1}$ ； $d_m$ ——轴承节圆直径，mm； $f_0$ ——与轴承类型相关的参数，可查表选取。

2.2. 定子、转子生热计算

电机的损耗发热主要由机械损耗、磁损耗和电损耗组成的，附加损耗占大约1%~5%，各部分损耗有着对应的理论计算公式，但用理论公式对电主轴的功率损耗进行求解需要查找大量的计算参数，因此本文采用了一种在工程实践中简便易行的功率转化法估算电主轴的功率损耗：

$\{\begin{array}{l}{P}_{loss}=P\cdot \left(1-\eta \right)\\ P_{定子}=2/3P_{\text{loss}}\\ P_{转子}=1/3P_{\text{loss}}\end{array}$ (3)

其中， $P$ ——电机的额定功率； $\eta$ ——电机效率。

电主轴的损耗发热是以生热率的形式加载在模型上，生热率的计算公式为：

$q=\frac{Q}{V}$ (4)

式中， $Q$ ——产生的热量； $V$ ——热源体积。

2.3. 电主轴传热的分析与计算

2.3.1. 定子和冷却油之间的对流换热

电主轴定子和冷却油之间换热属于流体强制对流换热，冷却油的流动状态不同，其换热系数也不同。因此需要根据雷诺数判断冷却水的流动状态，对流换热系数可由下式计算：

当Re < 2200时，流体处于层流状态

$a=\frac{\lambda }{L}{\left(1.86RePr\frac{D_H}{L}\right)}^{\frac{1}{3}}$ (5)

当2200 < Re < 10,000时，流体处于层流向紊流过渡的状态

$a=\frac{\lambda }{L}\left(0.023R{e}^{0.8}P{r}^{0.4}\right)$ (6)

当10,000 < Re时，流体处于紊流状态

$a=\frac{\lambda }{L}\left(R{e}^{\frac{2}{3}}-125\right)P{r}^{\frac{1}{3}}{\left(1+\frac{D_H}{L}\right)}^{\frac{2}{3}}$ (7)

其中：Re——雷诺常数，Pr——普朗特常数，D_H——冷却液当量直径，m； $\lambda$ ——流体的导热系数 ${\text{W/(m}}^{\text{2}}\cdot \text{k)}$ ；L——特征长度，m。

2.3.2. 电主轴定子与转子之间气隙的对流换热

转子通过气隙把部分热量传送到定子，需要根据雷诺数计算定转子间气隙的对流换热系数。当气隙处于纯层流状态，此时热量将会以纯导热的形式从一个表面传送到另一个表面。当气隙处于层流和紊流状态时，对流换热系数可由下式计算：

$a=\frac{\lambda }{h}0.23{\left(\frac{\delta }{r_1}\right)}^{0.25}R{e}^{0.5}$ (8)

式中， $\delta$ ——定转子间的气隙厚度； $\lambda$ ——流体的导热系数； $r_1$ ——转子外圈半径； $h$ ——气隙厚度。

2.3.3. 电主轴端部和周围空气的对流换热

电主轴端部与周围空气之间的对流换热系数可以看成一个用空气流量和主轴转速来描述的多项式函数：

$a=\left(c_0+c_1{u}^{c_2}\right)$ (9)

式中，c₁、c₂、c₃——通过试验总结的常数；u——电主轴的速度，m/s。

2.3.4. 轴承与冷却油间的对流换热

在轴承对流换热系数的计算中，目前常用的确定方法为：比拟法、分析法、数值法和分析法 [11] [12]。Harris在层流假设的基础上提出了轴承对流换热系数的计算模型，轴承和冷却油之间的强迫对流换热系数计算公式为：

$a=0.332\frac{k}{x}P{r}^{\frac{1}{3}}R{e}^{\frac{1}{2}}$ (10)

式中，k——冷却油的热导率；x——轴承节圆直径；Pr——冷却油的普朗特数；Re——雷诺数。

3. 电主轴热特性的有限元仿真分析

3.1. 电主轴的有限元网格模型

本文的研究对象是TC-E650型数控雕铣机，主轴两端是7008C/7006C-2RZ-P4角接触球轴承，安装方式为双联背靠背。为了提高计算效率，在进行网格划分前，先对几何模型进行适合的简化处理：忽略电主轴上的螺纹孔、通油孔、倒角及其它的细小结构；忽略电主轴前端的拉刀结构和后端气动装置；定子和转子简化为厚壁圆筒。把简化后的几何模型导入Ansys Workbench，定义好零件各部分的材料属性，并对电主轴内部相邻的零件之间建立接触对，其中轴承内外圈滚道与滚珠之间设置为摩擦接触，摩擦系数设定为0.2，其它不发生相对运动的零件设置为绑定接触。在网格划分时，由于定子和转子的几何形状比较规则，采用六面体主导网格划分，其它零件采用自动网格划分的方法。此外，为了平衡计算效率与有限元分析的精度，设定整体网格单元尺寸为3 mm，但对主轴、定转子的网格单元进行加密处理，设定网格单元尺寸为2 mm。划分后的网格如图1所示。

在Ansys Workbench的网格统计中可以看出，网格划分的节点数是536,412，单元数为379,599，网格单元的平均质量为0.82729，大于规定的最小值0.7，因此符合要求，能够用于进行接下来的热仿真。

Figure 1. Finite element meshing model of motorized spindle

图1. 电主轴有限元网格模型

3.2. 电主轴的瞬态热分析

在Ansys中建立电主轴的瞬态热分析系统，初始条件设定为：电主轴转速为10,000 r/min，实验环境温度为22℃，冷却液温度为24℃。根据上文的公式计算的电主轴生热及散热参数汇总如表1所示。

Table 1. Statistical table of heat dissipation boundary condition parameters of motorized spindle

表1. 电主轴的热边界条件参数统计表

瞬态分析设置的仿真时间3小时，时间步长为100 s，时间子步为108个，模拟电主轴从冷态到稳态整个过程的工作状态。由于瞬态热–结构耦合分析对网格质量的要求高，整个模型的网格精细复杂，需要处理及存储的数据量巨大，为保证有限元分析软件的正确运行，外接大容量的移动硬盘存储数据。进行热–结构耦合分析前需要先仿真出电主轴每一时刻的瞬态温度场，图2为电主轴在运行第5000 s时的温度场仿真结果。

从图3中可以看出，电主轴在运行第5000 s时的温度最高点为51.37℃，位于前端角接触球轴承的内圈，造成这样的原因一方面可能是轴承的内圈与滚珠相互摩擦生热量大，另一方面轴承的内圈相对封闭，散热能力较差；温度最低点为27.77℃，位于定子外面的螺旋冷却管道上，这是由于冷却液的循环流动，不断吸收定子产生的热量，使得此点温度明显低于电主轴其它部分的温度。

将全部时刻的温度场导入结构分析模块并添加相应的约束条件，求解后可以得到对应时刻下电主轴的热变形云图，图3所示为电主轴在运行第5000 s时的整体热变形云图，图4、图5、图6所示分别为在X、Y、Z方向上的热变形云图。

Figure 2. Temperature field after running 5000 seconds

图2. t = 5000 s时的温度场

Figure 3. Thermal deformation after running 5000 seconds

图3. 电主轴在t = 5000 s时的整体热变形云图

Figure 4. Thermal deformation in X direction after running 5000 seconds

图4. 电主轴在t = 5000 s时的X方向热变形云图

Figure 5. Thermal deformation in Y direction after running 5000 seconds

图5. 电主轴在t = 5000 s时的Y方向热变形云图

Figure 6. Thermal deformation in Z direction after running 5000 seconds

图6. 电主轴在t = 5000 s时的Z方向热变形云图

通过选取电主轴内零部件为分析对象，重新分析后可以更直观的看到电主轴内部各位置的热变形情况，如图7所示为主轴–轴承系统的热变形云图，图8所示为转子的热变形云图。

Figure 7. Thermal deformation of spindle-bearing system after running 5000 seconds

图7. 主轴–轴承系统在t = 5000 s时的热变形云图

Figure 8. Thermal deformation of rotor after running 5000 seconds

图8. 转子在t = 5000 s时的热变形云图

上图有限元模型中的坐标系继承的是Solidworks中的坐标系，因此其X轴方向对应电主轴在实际加工中的轴向方向。根据有限元仿真分析的结果可以看出，在电主轴运行第5000 s时的最大热误差为40.52 μm，最大热误差出现在主轴的末端位置，轴向热误差最大为40.19 μm，其它两个方向的热误差为9.79 μm和9.59 μm。此外，定子的热误差主要集中在定子两端的位置，且定子前端的热变形比定子后端的热变形大，最大为17.347 μm，原因可能是转子两端距电主轴的前后轴承组较近，而前端轴承组的发热量又高于后端轴承组。

利用Ansys Workbench的后处理功能可以得到电主轴从冷态到稳态整个过程中每个位置处的热变形数据，提取电主轴主要零部件及主轴末端的相关数据进行处理后，结果如图9、图10所示。

根据图9的仿真预测结果可知，电主轴的加工主轴热变形程度最大，先在前1500 s时迅速增加至35 μm，然后在约4500 s达到热稳态时的最大热变形为40 μm；从图10可以看出，电主轴的热变形主要集中在轴向上，与整体热误差保持一致，达到热稳态时的最大热变形为52 μm，X、Y方向上的热变形较小，达到热稳态时大约都是16 μm。由于电主轴的结构在X和Y方向是对称，在相应的约束条件下，热误差累积表现在了轴向方向上。

Figure 9. Thermal deformation of main components of motorized spindle

图9. 电主轴主要部件的热变形

Figure 10. Total and anisotropic thermal deformation of motorized spindle end

图10. 电主轴末端整体及各向热变形

4. 实验验证

根据上述的仿真结果，本次实验主要对主轴末端的轴向热误差进行测量。测量前先把千分表固定在雕铣机的工作台上，并调整千分表的测杆轴线垂直于电主轴的轴线，然后调零位，使千分表测头与主轴前端面处于刚好接触的状态，并转动刻度盘使零线与大指针对齐。具体安装如图11所示。

Figure 11. Spindle thermal error experimental measurement

图11. 电主轴热误差实验测量现场

实验开始时，先把千分表调零位置时机床的机械坐标设置为参考坐标，然后控制电主轴运动的G代码首先令电主轴沿Z轴正方向上移5 mm，使主轴前端面与千分表测头分离并保持安全距离，然后以10,000 r/min的转速运动5 min，停转后沿Z轴负方向下移5 mm回到参考坐标的原点位置，电主轴静止10秒并读出此时千分表的数值，即为机床运行5分钟的热位移值，然后程序使电主轴再次沿Z轴正方向上移5 mm，如此循环往复，直到测得的热变形数据保持稳定，即机床达到热稳态为止。测得的实验结果及上文的仿真结果如图12所示。

Figure 12. Simulation and experimental results of thermal error

图12. 主轴末端热误差的仿真与实验结果

5. 结论

本文运用Ansys Workbench有限元软件对数控雕铣机的电主轴单元进行了瞬态热–结构耦合分析，构建了高速电主轴热误差预测有限元模型，并通过实验室的TC-E650型数控雕铣机对仿真结果进行验证，得到的主要结论有：

1) 电主轴在开始运行的前1500 s内热误差变化幅度较大，在短时间内迅速增长到40 μm，而后缓慢增长，在约4500 s达到热稳态，此时的主轴末端的最大轴向热误差为52 μm。

2) 电主轴的热误差主要集中在轴向上，X、Y方向的热误差几乎相同，达到稳态后大约为15 μm，相比轴向热误差要小很多。转子稳态时的热误差大约为20 μm，定子由于受到螺旋冷却管道中冷却液的影响，热变形最小为3 μm。

3) 实验有效验证了电主轴热误差有限元模型有良好的预测精度，为后续的热误差的补偿及热平衡设计提供了参考。

NOTES

^*通讯作者。