上海理工大学机械工程学院，上海

收稿日期：2023年12月18日；录用日期：2024年3月12日；发布日期：2024年3月19日

摘要

针对电主轴热误差仿真建模中忽略冷却液自身的温度变化对冷却效果的影响，导致预测精度不准确的问题。建立流–热–固多物理场耦合有限元模型，在热–结构耦合模型中引入流体域，通过流–固交界面将冷却液的温度场作为热分析中螺旋冷却管道的精确热边界条件、压力场作为结构分析中的力约束，更为真实地反映冷却液的实际换热能力和螺旋冷却管道的实际受力情况。采用基于热流密度的载荷施加方式与流–热–固多物理场耦合方法相结合的ANSYS仿真模型，仿真结果表明电主轴最大轴向热变形为51 μm，与实验测得的值相比，残差在5.2 μm左右。可以有效地实现高速电主轴热误差预测。

关键词

电主轴，热误差，流–热–固耦合，热流密度，ANSYS

Research on Thermal Error Modeling of Electric Spindle Based on Fluid Thermal Solid Multiphysics Coupling

Yushan Sun, Guangsheng Chen, Haoshuo Wang

School of Mechanical Engineering, University of Shanghai for Science and Technology, Shanghai

Received: Dec. 18^th, 2023; accepted: Mar. 12^th, 2024; published: Mar. 19^th, 2024

ABSTRACT

In the simulation modeling of thermal errors in electric spindles, the influence of the temperature change of the coolant itself on the cooling effect is ignored, resulting in inaccurate prediction accuracy. Establish a fluid thermal solid multi-physics coupling finite element model, introduce a fluid domain into the thermal structure coupling model, and use the fluid solid interface to treat the temperature field of the coolant as the accurate thermal boundary condition for the spiral cooling pipeline in thermal analysis, and the pressure field as the force constraint in structural analysis. This more accurately reflects the actual heat transfer capacity of the coolant and the actual stress situation of the spiral cooling pipeline. The ANSYS simulation model combines the load application method based on heat flux density with the fluid thermal solid multi-physics coupling method. The simulation results show that the maximum axial thermal deformation of the electric spindle is 51 μm. Compared with the experimentally measured values, the residual is 5.2 μm around, it can more accurately predict the thermal error of the electric spindle.

Keywords:Electric Spindle, Thermal Error, Flow Thermal Solid Coupling, Heat Flux, ANSYS

Copyright © 2024 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

数控机床作为工业母机，代表着一个国家制造业的总体技术水平 [1] 。在高速高精度的加工条件下，热变形误差是影响精密机床加工精度的关键因素，由热变形造成的误差大约占到机床总误差的40%~70% [2] ，综合热误差的比重和建模有效性考虑，针对热误差进行建模与补偿对于减小机床的误差、提高机床的精度有着重要意义。

电主轴热误差主要采用基于数值模拟方法的有限元方法和采用数据驱动两种方法进行建模 [3] 。对于数据驱动方法，Chen研究了隐含层的节点数量和学习速度对神经网络模型的影响，并分析了学习速度与模型收敛性之间的具体关系 [4] 。李彬使用遗传算法对神经网络作出优化，从而建立电主轴的热误差补偿模型 [5] 。姜磊采用了粒子群算法优化BP神经网络，构造机床热误差优化目标函数从而提高模型的预测精度 [6] 。但是数据驱动模型需要大量的数据进行训练，且所建立的模型只在采集数据的机床上适用。有限元方法以热力学的相关理论为基础，具有较强的通用性。Fang等人将有限元法与深度学习组合应用，对复杂工况下的电主轴热变形进行预测 [7] 。Yang等人基于逆热传导理论实现对流换热系数的精确求解，在Ansys Workbench中对主轴的温度场和热误差进行了预测 [8] 。Wang等通过ANSYS Workbench软件对电主轴进行瞬态热–结构耦合分析，从而预测出电主轴热误差随时间变化的趋势 [9] [10] 。Andreas Naumann提出了一种多速率高阶时间步进法，对机床的有限元模型进行实时精确的瞬态热变形仿真 [11] 。

本文将重点以热流密度的方式下的流热固多物理场展开理论推导及建模。考虑了轴承内部的实际温度梯度，对轴承总的热量进行精细化分配，计算出每个接触面的热流密度。通过建立流–热–固耦合模型得到较为准确的边界条件，从而提高有限元模型的预测精度。

2. 电主轴边界条件的确定

2.1. 定、转子生热与轴承热流密度计算

2.1.1. 定、转子生热率计算

由于电机损耗发热的理论公式需要大量的参数，且较多参数不易确定，因此本文采用功率转换法来近似计算电主轴的功率损耗：

$\{\begin{array}{l}{P}_{\text{loss}}=P\cdot \left(1-\eta \right)\\ q=Q/V\end{array}$ (1)

式中P表示电机的额定功率；η表示电机效率，Q表示产生的热量；V表示热源体积。

2.1.2. 轴承热流密度计算

Palmgren通过实验总结出轴承摩擦力矩计算的经验公式，如式2所示。

$\{\begin{array}{l}Q=1.047\times {10}^{-4}Mn\\ M=M_0+M_1\end{array}$ (2)

式中，Q表示轴承的生热量，W；M表示轴承的总摩擦力矩，N∙mm；n表示轴承转速，r/min； $M_0$ 表示与轴承类型有关的黏性摩擦力矩，N∙mm； $M_1$ 表示与负载有关的负载摩擦力矩，N∙mm。

通过查阅资料和相关参数求解出角接触球轴承的总发热量后，先对轴承总的热量进行精细化分配，再计算出每个接触面的热流密度。以前端角接触球轴承为例，其内部各零件热流密度的理论计算公式如下：

$\{\begin{array}{l}{Q}_1=Q_2=\frac{Q}{4}\\ Q_3=\frac{Q}{2Z}\\ q=\frac{Q}{S}\end{array}$ (3)

其中， $Q_1$ 、 $Q_2$ 为轴承内圈滚道及外圈滚道分配的发热量，W； $Q_3$ 为轴承单个滚珠分配的发热量，W；Z为轴承的滚珠个数。q为热流密度，W/m²，S为面积，m²。

2.2. 电主轴对流换热计算

高速电主轴在传热途径上主要包括固体与固体之间的热传导、固体和流体之间的热对流、以及空间中互不接触的两物体之间的热辐射三种基本方式。热量在传递过程通常需要遵循傅里叶定律、牛顿冷却公式、以及斯忒藩–波尔兹曼定律三大定律。

本文所要确定的对流换热系数主要有定子与冷却液之间的对流换热、定子与转子之间气隙的对流换热、轴承与润滑脂之间的对流换热、电主轴端部与周围空气的对流换热。具体的计算公式和步骤可在文献 [9] 中得到，本文对此不再赘述。

3. 多物理场耦合热模型原理分析

3.1. 多物理场耦合基本方程

有限元法是将求解的物理温度或变形随时间、空间连续分布的问题，离散为在时间域和空间域一定数量的单元求解，再对每一个单元选择插值函数求解，使有限元结果收敛于实际的精确解。以下分别为温度场微分方程、流体流动的运动方程和结构分析方程。

温度场微分方程

导热体中离散出来的微元六面体，边长分别为 $\text{d}x$ 、 $\text{d}y$ 、 $\text{d}z$ ，微元体内发生的热量是关于局部温度t的函数，而导热率λ、比热c和密度ρ都是常量。根据傅里叶定律可以求出热流量在三个方向分量上的热量：

$\{\begin{array}{l}\text{d}{Q}_x=-\lambda \frac{\partial t}{\partial x}\text{d}y\text{d}z\text{d}\tau \\ \text{d}{Q}_y=-\lambda \frac{\partial t}{\partial y}\text{d}x\text{d}z\text{d}\tau \\ \text{d}{Q}_z=-\lambda \frac{\partial t}{\partial z}\text{d}x\text{d}z\text{d}\tau \end{array}$ (4)

在 $x=\text{d}x$ ， $y=\text{d}y$ ， $z=\text{d}z$ 三个表面上流出微元体的热量简化后可以得到：

$F\left(x+\text{d}x,\tau \right)=F\left(x,\tau \right)+\frac{\partial F}{\partial x}\text{d}x=\lambda \frac{\partial t}{\partial x}+\frac{\partial }{\partial x}\left(\lambda \frac{\partial t}{\partial x}\right)\text{d}x$ (5)

根据能量守恒定律，可以最终得到下面的平衡方程式：

$\frac{{\partial }^{2}t}{\partial x^2}+\frac{{\partial }^{2}t}{\partial y^2}+\frac{{\partial }^{2}t}{\partial z^2}+q_V^m=\frac{1}{\alpha }\frac{\partial t}{\partial \tau }$ (6)

式中， $\alpha =\lambda /c\rho$ 表示热扩散率，其物理意义是热量传递的速度，m²/s。c表示比热容，J/(Kg∙˚C)；ρ表示密度，Kg/m³， $q_V^m$ 表示单位面积时间内热源的发热强度，单位为J/m³。

式(6)为均匀介质中导热微分方程的普遍表达形式，可根据系统中是否存在内热源作进一步的简化。热分析的有限元方程求得：

$C\stackrel{\dot{}}{T}+K_{T}T=Q$ (7)

式中，C表示比热矩阵；K_T表示热传导矩阵；T表示节点温度向量； $\stackrel{\dot{}}{T}$ 表示节点温度变化率向量；Q表示热通量向量。

3.2. 电主轴流热固多物理场耦合的求解方法

耦合问题在求解时也有两种方法，第一种是基于单向耦合的交叉迭代求解，第二种是基于双向耦合的直接求解。双向耦合问题的求解步骤复杂，对有限元网格模型的精度要求高，求解占用的时间久，而单向耦合相对而言简便易行，求解时间适中。因此，本文将采用单向耦合的方式建立多物理场耦合有限元模型，实现电主轴热变形的预测。

机床电主轴在高速旋转过程中，产生的热量通过电主轴的固体域传递到冷却液与冷却套之间的流固交界面以及通过电主轴表面传递到周围空气中，两种对流方式分别为强迫对流和自然对流换热。电主轴的流热固耦合分析中的热量流动如图1所示。

图1. 电主轴流热固耦合分析中的热量流动图

Ansys Workbench软件集成了结构静力学、稳瞬态热分析模块、Fluent及CFD流体分析模块，能够进行多种物理场的耦合仿真分析。如图2所示为电主轴流热固耦合分析的具体流程。

图2. 电主轴流热固耦合分析流程图

4. 电主轴流热固多物理场耦合有限元分析

4.1. 电主轴有限元模型的建立

在运用Solidworks三维建模软件建立电主轴的几何模型的过程中，为了平衡有限元模型求解的精度与效率，对模型进行了忽略和简化处理。忽略电主轴上的螺纹孔、通油孔、倒角等细小结构；忽略电主轴前后端的连接装置；将冷却液槽简化为螺旋环形结构、定、转子简化为厚壁圆筒。将处理好的电主轴模型导入ANSYS Workbench软件中，并设置电主轴系统的主要零部件材料属性。

利用Fluent中的填充功能根据电主轴的模型生成冷却液的模型，将其设定为流体域。为了使仿真更好地模拟实际工况，还需要设置边界层，将其总厚度设置为1 mm，层数为5，膨胀率为1.1。之后对冷却液的进口处和出口处进行定义，并对流体域进行网格划分，由于Fluent中对流体域的网格质量要求高于固体域的网格质量，因此，对初步划分的网格单元进一步细化，将流体域的单元尺寸加密为1 mm，划分后含流体域的电主轴有限元网格模型如图3所示。

图3. 含流体域的电主轴有限元网格模型

在ANSYS Workbench的流–固耦合分析中，流体分析所得到的数据可以用温度和对流两种方式加载在热–结构耦合分析模块中，本文选取温度的形式作为电主轴热分析中的精确边界条件代替由经验公式计算的固定值进行分析。冷却液温度边界条件的导入如下图4所示：

图4. 冷却液温度边界条件的导入

为了模拟冷却液的流动，需要确定合适的湍流模型。本文选择RNG $k\text{-}\epsilon$ 模型对冷却液进行仿真，打开能量方程并定义冷却油的密度、热传导率、比热容及粘度等，冷却油的各项物理参数如表1所示。

表1. 冷却油各项物理参数的设置

在Ansys中建立电主轴的瞬态热分析系统，机床电主轴以10,000 rpm运行3小时左右，室内环境温度保持在22℃，冷却油温度设定为24℃。计算的电主轴热流密度及散热参数如表2和表3所示：

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

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

将上表汇总的电主轴的各部分对流换热系数添加到相应的位置，按热流密度载荷施加到前后端的角接触球轴承上，定、转子都按生热率施加，具体如下图5：

图5. 热边界条件施加示意图

在Ansys Workbench中建立电主轴的瞬态热分析系统，将稳态热分析的网格模型及材料参数传递到瞬态分析中。本次瞬态分析设置的仿真时间14,400 s，时间步长为100 s，时间子步为144个，模拟电主轴从冷态到热稳态的整个过程，图6为电主轴在t = 3000 s时的瞬态温度场。冷态开始，电主轴因摩擦、电机生热等原因使得电主轴温度急速上升，主要集中在轴承、定转子处。随着时间的增加，由于各部件间的对流换热，电主轴内部零部件之间的温度梯度逐渐减小，整体温度场趋于统一、平稳；在3000 s时最高温度为70.7℃，最高温度主要集中在滚珠位置，在热流密度施加方式中热量主要集中在滚珠处导致此点温度最高。

将热分析中的144个时间节点的温度场全部导入ANSYS Workbench中的结构分析模块中并设置适当的约束条件，求解后得到如图7和图8所示的热流密度载荷施加方式下电主轴在t = 3000 s时整体及X方向上的热变形云图。由于建立的Solidworks三维模型导入ANSYS中时坐标系选定有所不同，使得实际的机床电主轴的Z方向在ANSYS中是X方向，即电主轴的轴向方向对应着云图中的X方向，云图中的Z方向对应着实际的X方向。由热变形云图可知，热流密度施加方式下在t = 3000 s时的电主轴整体最大瞬态热变形为50.86 μm，在Z方向上的最大瞬态热变形是50.34 μm。轴承、主轴、转子三者紧密接触，因此，轴承热源载荷的不同施加方式对这三者的温度及热变形影响最大。由于X，Y方向热变形相对较小，本文只对电主轴Z方向的热变形进行研究。

图6. t = 3000 s时电主轴瞬态温度场

图7. t = 3000 s时电主轴整体热变形

电主轴瞬态热变形随时间变化的曲线如图9所示，基于热流密度的多物理场耦合有限元热变形模型在前3000 s内各部件温升明显，导致热变形变化幅度非常大，当达到约4000 s时各零部件基本达到稳态，主轴达到稳态时的轴向热变形预测结果是51 μm。整体热变形和轴向热变形趋势几乎一致，预测结果整体热变形略微大于轴向热变形。由此可见电主轴热变形主要表现在轴向热伸长。

4.2. 实验验证

电主轴热态特性实验平台是由硬件部分和软件部分共同搭建的，硬件部分主要是由单片机、模拟量采集器、基恩士电涡流传感器、温度传感器DS18B20、RS485通讯接口等组成。软件部分使用C#编写上位机软件，对本次实验的电主轴温度、热位移数据进行同时采集显示及存储处理。位移传感器为非接触式测量，测量行程为2 mm，调整检测棒末端与位移传感器的距离为1 mm，具体的实验平台如图10所示。

图8. t = 3000 s时电主轴轴向热变形

图9. 电主轴轴向与整体热变形

实验开始前，先编写G代码使其先从原点位置运行到距离传感器1 mm中心位置，停留20 s，测量初始位移，然后将Z轴沿正方向向上移动4 mm，使位移传感器超出测量设置的极限即2 mm，不触发采集程序。然后以10,000 rpm的转速运行14,400 s，采集间隔设置为1 s，当转动时间达到100 s时，Z轴向下移动4 mm，回到初始测量位置，停留5 s，触发采集程序，开始同时采集位移和温度数据，如此往复100次。测得的实验结果和仿真值的结果如图11所示：

图10. 热误差实验测量平台

图11. 电主轴轴向热误差预测值与实验测量值

5. 结论

根据多物理场耦合原理确定了电主轴流–热–固耦合求解的方式及流体分析中的湍流模型，并通过流–固交界面将冷却液的温度、压力数据传递给定子冷却套作为电主轴流–热–固耦合分析中的精确边界条件。对电主轴的角接触球轴承使用热流密度施加方式建立了电主轴的流–热–固多物理场耦合有限元模型。并对电主轴从冷态到热稳态整个时间段中的瞬态温度场及热变形进行预测，该模型充分考虑了电主轴单元的边界条件。通过ANSYS Workbench软件实现了所建模型的仿真，其模型预测高速10,000转热伸长从冷态到稳态呈现前4000秒快速增大并趋于稳态，到4000秒后逐渐呈现稳态，其热误差在51 μm左右。通过实验验证模型的预测残差最大在5.2 μm以内，验证了电主轴流–热–固耦合有限元模型具有较高的预测精度，为机床电主轴热误差建模和补偿提供指导意义。

文章引用

孙玉山,陈光胜,王豪硕. 基于流热固多物理场耦合的电主轴热误差建模研究
Research on Thermal Error Modeling of Electric Spindle Based on Fluid Thermal Solid Multiphysics Coupling[J]. 建模与仿真, 2024, 13(02): 1236-1246. https://doi.org/10.12677/MOS.2024.132116

参考文献