多边界寻优的电主轴热特性数字孪生

Digital Twin for Thermal Characteristics of Motorized Spindle Based on Multi Boundary Optimization

摘要: 针对电主轴热特性有限元仿真误差大，热边界条件无法精确辨识问题，提出数字孪生驱动的电主轴热特性在线监测方法。搭建基于C++和Qt的电主轴热特性数字孪生系统，建立基于最小二乘优化的接触热阻与对流换热系数目标函数，根据关键测温点的实测与仿真温度实时优化热边界条件并输出至ANSYS Parametric Design Language (APDL)，通过后台调用ANSYS完成电主轴热特性分析及数据处理，实现物理实体与虚拟模型实时映射的电主轴热特性数字孪生。实验结果表明，该方法有效提高电主轴热特性分析精度，温度场预测精度达98%，热变形预测精度达96%，为电主轴热优化设计及热误差控制提供依据。

Abstract: In order to improve the finite element simulation accuracy of thermal characteristics of motorized spindle and the identification accuracy of thermal boundary conditions, an online monitoring method of thermal characteristics of motorized spindle driven by digital twin is proposed in this study. A digital twin system for thermal characteristics of motorized spindle is built using C++ and Qt. The objective functions of thermal contact resistance and convective heat transfer coefficient are established based on least square optimization method. The thermal boundary conditions are real- timely optimized using the measured and simulated temperatures of thermal key points and output to ANSYS Parametric Design Language (APDL). The thermal characteristic analysis and data processing of the motorized spindle are completed by calling ANSYS in the background. The digital twin for thermal characteristics of the motorized spindle is realized through real-time mapping of physical entities and virtual models. The experimental results show that the simulation accuracy of thermal characteristics of the motorized spindle is effectively improved using the proposed digital twin-driven online monitoring method, the prediction accuracy of temperature field reaches 98%, and the prediction accuracy of thermal deformation reaches 96%, which provides a basis for thermal optimization design and thermal error control of the motorized spindle.

1. 引言

作为高速机床核心部件的电主轴，其高速运转产生的热量无法及时导出，形成热误差，影响零件加工精度。因此，对电主轴热特性进行精确分析并控制是解决电主轴热变形的关键技术。目前，电主轴热特性分析方法主要是有限元仿真法，通过建立有限元模型，施加热边界条件，应用有限元分析软件获得电主轴温度场及热变形。Uhlmann [1] 等考虑生热量以及零部件之间的接触热阻和对流传热搭建电主轴三维有限元模型，对其热特性进行瞬态模拟。Su [2] 等结合过盈配合、热效应在部件中的耦合作用，建立综合优化的电主轴模型，研究不同变量下电主轴温升的变化。有限元仿真法可以定性模拟特定条件下电主轴热特性，而电主轴热特性就像“黑匣子”，有诸多边界条件无法精确获得，限制了热特性分析的准确性。

内热源生热量、对流换热系数、结合面接触热阻是热特性分析的三个关键参数，其辨识精度决定热特性分析的准确性。为获得更加贴近实际的边界条件，康惠民 [3] 等分析电主轴轴承滚动体的受力情况，根据Palmgren经验公式得到轴承摩擦生热量，通过传热分析建立轴承温度场模型，实验验证轴承温升与滚动体运动状态的关系。陈小安 [4] 等考虑实际功率下各种损耗，确定电主轴生热量等热边界条件，对电主轴的温度场进行分析。张丽秀 [5] 等通过试验测得电主轴总体能量损耗，计算电主轴定/转子及轴承的生热量，通过试验数据优化电主轴的热边界条件，进而得到较为精确的电主轴温升模型。孟令聪 [6] 考虑冷却液实时传热，建立电主轴热–流–固耦合模型，同时利用遗传算法对电主轴的生热量与换热系数进行修正，将修正后的热边界条件重新加载，得到仿真精度较高的温度场模型。

现有研究为揭示电主轴热特性提供依据。但是，由于热边界条件受材料、磨损、接触方式等多种因素影响，通过经验公式计算的热边界条件与实际热边界有一定的误差。导致热边界条件难以准确辨识的主要因素是界面及结合面参数难以定量获得，因此，探索一种多边界寻优的热边界条件辨识方法，可大幅提高热边界条件的辨识精度。在此基础上，将数字孪生技术应用于电主轴热特性分析，通过边界监测及映射搭建热特性数字孪生系统，在虚拟空间中定量分析电主轴热特性，为电主轴热优化设计及热误差控制提供依据。

2. 电主轴热边界优化

图1为热边界优化流程图，将电主轴关键测温点的仿真温度及实测温度输入目标函数，判断目标函数值是否超出阈值，如果超出阈值，则根据优化函数对电主轴热边界条件进行优化，如果小于阈值，则输出电主轴温度场及热变形。电主轴多边界优化的关键是目标函数的建立及关键测温点的选取，根据传热理论，热边界条件与热流密度及边界温度有关，因此，选取组成结合面的界面为关键测温点，通过嵌入微型传感器实时监测边界温度，根据测得的边界温度可由傅里叶定律 [7] 计算流过结合面的热流密度。

$q=\lambda \cdot \text{grad}T$ (1)

式中：q为流过结合面的热流密度，W/m²；λ为导热系数，W/m·K；gradT为温度梯度。根据计算所得流过结合面热流密度及关键测温点温度可对接触热阻及对流换热系数进行优化。

Figure 1. Flowchart of thermal boundary optimization

图1. 热边界优化流程图

2.1. 接触热阻优化

接触热阻的定义为两物体接触界面的温差与通过两物体名义接触面处的热流密度之比 [8] [9] ，

$R_c=\frac{\Delta T}{q}=\frac{T_1-T_2}{q}$ (2)

式中：R_c为接触热阻，m²·K/W；ΔT为两界面温差，℃；q为流过结合面的热流密度，W/m²；T₁、T₂分别为接触面两侧的温度，℃。

由式(2)可知，要计算结合面接触热阻，需要测得结合面两侧的界面温度，由于电主轴的主要结合面包括轴承–轴承套、轴承套–冷却套以及定子–冷却套等，因此，上述结合面的界面即为关键测温点，通过嵌入微型传感器可测量上述关键测温点的温度。由于结合面接触热阻随电主轴温度的升高而变化，而且要保证仿真结果与实际一致，因此，需要通过仿真与实测温度联合对接触热阻进行实时优化，以获得更贴近实际的电主轴热特性。

接触热阻优化的目标应满足以下两点：一是结合面温差的仿真值与实测值之差最小，二是结合面界面温度的实测值与仿真值之差最小。为获得更加精确的接触热阻，本文基于最小二乘法对其进行优化，采用结合面实测温差ΔT_m与仿真温差ΔT_s差值，以及界面实测温度T_m与仿真温度T_s差值的平方和 ${\displaystyle {\sum }_{i=1}^n{r}^2}$ ( $i=1,2,\cdots ,n$ )来衡量误差的大小，构建目标函数：

$F\left(T_s,T_m\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{r}^2}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\left(\Delta T_{si}-\Delta T_{mi}\right)+\left(T_{mi}-T_{si}\right)\right)}^2}$ (3)

式(3)中，根据结合面的数量取n = 5。根据接触热阻与界面温差的关系 $\Delta T_{si}=R_{si}\cdot q_{si}$ 可得接触热阻目标函数：

$R_{\min }={\displaystyle \underset{i=1}{\overset{n}{\sum }}{r}^2}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\left(R_{si}\cdot q_{si}-\Delta T_{mi}\right)+\left(T_{mi}-T_{si}\right)\right)}^2}$ (4)

对式(4)求R_si的最小值：

$\frac{\partial R}{\partial R_{si}}=2{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(\left(R_{si}\cdot q_{si}-\Delta T_{mi}\right)+\left(T_{mi}-T_{si}\right)\right)\cdot q_{si}}=0$ (5)

由式(5)可得接触热阻优化函数：

$R_{si}=\left(T_{si}-T_{mi}+\Delta T_{mi}\right)/q_{si}$ (6)

式(6)中，流过结合面的热流密度可由式(1)计算，结合面的界面温度分别通过仿真和实验获得，因此，应用式(6)可实现电主轴结合面接触热阻的优化。

2.2. 对流换热系数优化

对流换热系数定义为当流体与固体表面之间的温度差为1 K时，一平方米壁面面积在每秒所能传递的热量，可由牛顿冷却方程 [10] [11] 计算：

$h_c=\frac{q}{\Delta T}=\frac{q}{T_1-T_2}$ (7)

由式(7)可知，要计算对流换热系数，需要测得冷却套界面温度、流体温度及热流密度，根据1.1节所述传感器布置，冷却套温度可由嵌入的温度传感器测得，流体温度可通过流体出口温度获得，热流密度可由式(1)根据关键测温点温度计算。由于对流换热系数与冷却介质的温度、流速等参数有关，不同条件下的对流换热系数会有所不同，为获得更贴近实际的电主轴热特性，需要通过仿真与实测温度联合对对流换热系数进行实时优化。

同接触热阻优化，对流换热系数优化的目标应满足以下两点：一是固流温差的仿真值与实测值之差最小，二是冷却套界面温度的实测值与仿真值之差最小。为获得更加精确的对流换热系数，本文基于最小二乘法对其进行优化，采用固流实测温差ΔT_m与仿真温差ΔT_s差值，以及冷却套界面实测温度T_m与仿真温度T_s差值的平方和 ${\displaystyle {\sum }_{i=1}^n{r}^2}$ ( $i=1,2,\cdots ,n$ )来衡量误差的大小，构建目标函数：

$F\left(T_s,T_m\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{r}^2}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\left(\Delta T_{mi}-\Delta T_{si}\right)+\left(T_{mi}-T_{si}\right)\right)}^2}$ (8)

式(8)中，根据对流界面的数量取n = 3。根据对流换热系数与固流温差的关系 $\Delta T_{si}=q_{si}\cdot /h_{si}$ 可得对流换热系数目标函数：

$h\left(T_s,T_m\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{r}^2}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\left(\Delta T_{mi}-q_{si}\cdot /h_{si}\right)+\left(T_{mi}-T_{si}\right)\right)}^2}$ (9)

对式(9)求h_si的最小值：

$\frac{\partial h}{\partial h_{si}}=2{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(\left(\Delta T_{mi}-q_{si}\cdot /h_{si}\right)+\left(T_{mi}-T_{si}\right)\right)\cdot q_{si}/h_{si}^2}=0$ (10)

由式(10)可得对流换热系数优化函数：

$h_{si}=q_{si}/\left(T_{mi}-T_{si}+\Delta T_{mi}\right)$ (11)

式(11)中，流过冷却套的热流密度可由式(1)计算，冷却套的界面温度分别通过仿真和实验获得，因此，应用式(11)可实现电主轴对流换热系数的优化。

3. 热特性数字孪生系统

本文采用C++和Qt Creator开发电主轴热特性数字孪生系统。首先使用SolidWorks对某型电主轴进行有限元建模；然后利用C++对热边界条件进行优化并生成APDL命令流，通过后台调用ANSYS 对电主轴进行热特性分析；使用Qt Creator开发热特性数字孪生系统；采用IoT数据采集系统，通过TCP/IP通信协议搭建电主轴物理空间与数字孪生系统的数据交互。

3.1. 数字孪生系统架构

如图2所示，电主轴热特性数字孪生系统构架由物理维度、孪生维度、虚拟维度三部分组成，三者相互关联。物理维度将采集到的数据传输至孪生维度；孪生维度根据实测温度及仿真温度通过热边界优化函数实时优化电主轴热边界条件，优化后的热边界条件输出至APDL并后台调用ANSYS进行电主轴热特性分析；虚拟维度实时展示电主轴的温度场及热变形。

Figure 2. Digital twin system for thermal characteristics of motorized spindle

图2. 电主轴热特性数字孪生系统构架

3.2. 物理维度构建

物理维度是热特性数字孪生系统的基础，包括设备层、采集层、传输层。设备层由电主轴、控制系统、冷却系统组成；采集层由温度传感器和位移传感器组成；传输层由IoT数据采集系统实时采集电主轴工作过程中的动态数据并上传至孪生维度的数据层。热特性数字孪生系统物理维度如图3所示，其中温度传感器嵌入到电主轴内部，以采集结合面的界面温度，位移传感器布置在芯轴的前端，以采集主轴热变形，传感器采集的数据由IoT数据采集系统传输至上位机PC，上位机根据采集到的数据对电主轴热边界条件进行优化，并后台调用ANSYS进行热特性仿真。

Figure 3. Physical dimension of digital twin for thermal characteristic

图3. 热特性数字孪生物理维度

3.3. 孪生维度研发

孪生维度由C++语言编译实现，通过QODBC在数据库中存储电主轴关键测温点的实测与仿真温度，并通过输入输出流传递到优化函数进行热边界优化；初始热边界条件及电主轴运行参数由QT开发的可视化窗口输入；电主轴热特性仿真由C++后台调用ANSYS通过APDL命令流文件实现。

如图4所示孪生维度的孪生流程，首先将物理维度中的主轴系统映射到孪生维度中，完成孪生层中电主轴热特性数字孪生模型的构建。逻辑处理是数字孪生系统正常运转的支撑，逻辑处理主要包括数据库的连接、多边界条件优化、ANSYS后台分析，数据库通过I/O命令流实现与各个模块之间的通信，多边界条件优化基于优化函数实现接触热阻与对流换热系数的修正，通过批处理加载至ANSYS后台进行热-结构耦合分析。分析后的结果由QODBC完成数据处理与存储，孪生维度的核心是孪生数据交互，孪生数据主要包括温度数据、位移数据、热边界数据，最终在虚拟维度实现孪生数据可视化。

Figure 4. Flow chart of twin dimension

图4. 孪生维度孪生流程图

物理维度与孪生维度的数据交互是电主轴热特性数字孪生系统的桥梁，数字孪生系统通过互联网协议(TCP/IP)实现与IoT数据采集系统之间的连接。两者之间使用Scoket通信传递信息，通过connectToHost()函数创建服务器于客户端之间的连接 [12] ；最后通过read()和write()函数来读写数据。数据传输原理如图5所示。

Figure 5. Schematic diagram of data transmission

图5. 数据传输原理图

3.4. 虚拟维度设计

虚拟维度由初始设置层、显示层、存储层组成。初始设置层将热边界条件初始值发送到APDL供ANSYS调用；显示层将实体电主轴工作过程中的温度与热变形映射到虚拟空间实时显示；存储层将电主轴温升、热变形、误差等数据存储到计算机磁盘。热特性数字孪生系统界面如图6所示，包括温度场云图显示区、热变形云图显示区、关键测温点温度及误差显示区、温度曲线显示区、连接设置及仿真设置等。通过虚拟维度可实时监测电主轴的热特性，为电主轴优化设计及热变形控制提供依据。

Figure 6. Digital twin interface for thermal characteristics

图6. 热特性数字孪生界面

热特性数字孪生系统主要对电主轴工作过程中的温度数据与热变形数据进行孪生，孪生空间通过数据采集系统将主轴热特性参数传输到数据库中进行热边界映射及修正，结合虚拟维度中的初始热特性参数生成APDL命令流文件，热特性数字孪生系统后台调用命令流文件完成热特性仿真，并将电主轴温度数据、热变形数据、温度云图、热变形云图实时传递到数字孪生系统界面中。

4. 实验验证

4.1. 实验设置

为验证多边界优化的电主轴热特性数字孪生，以某型电主轴为对象，应用上述开发的热特性数字孪生系统，对电主轴进行热特性数字孪生实验。孪生系统物理维度的传感器布置如图7所示，在电主轴前后轴承、定子、冷却套等结合面位置嵌入微型温度传感器，使用导热硅胶将温度传感器与被测位置紧密贴合 [13] ，实现1.2节所述电主轴结合面界面温度监测；为验证电主轴热变形的孪生精度，在电主轴轴端布置三个电涡流位移传感器，分别测量电主轴芯轴XYZ三个方向的热变形。实验环境温度为20℃，冷却系统进水口水温为20℃，冷却水流量为0.2 m³/min，电主轴转速为12,000 r/min。

Figure 7. Motorized spindle sensor layout

图7. 电主轴传感器布置

4.2. 热边界条件计算与优化

4.2.1. 初始边界条件计算

轴承的摩擦生热和电机的损耗发热是电主轴生热的主要热源，电主轴与外部环境、零件之间的热量交换是一个非常复杂的过程，主要包括热传导、热对流、热辐射三种形式。在进行模拟仿真前需要计算初始接触热阻、初始对流换热系数、初始生热量等热边界条件的数值，电主轴传热计算几何模型如图8所示。为了能够准确测量电主轴各零件的温度数据，将温度传感器合理的布置在前后轴承、前后轴承座、定子、冷却水套等结合面位置，具体热电偶测温位置如图8所示。

热特性有限元分析需要施加内热源生热量、结合面接触热阻、对流换热系数等边界条件。轴承生热量可由Palmgren经验公式计算：

$Q_f=1.047\times {10}^{-4}nM$ (12)

式中：Q_f为轴承生热量，W；n为主轴转子转速，r/min；M为总摩擦力矩，N·mm，包括粘性摩擦力矩M_n和载荷摩擦力矩M_z：

$\{\begin{array}{l}M=M_n+M_z\\ M_n={10}^{-7}{f}_0{\left(vn\right)}^{\frac{2}{3}}{D}_m^3\left(vn\ge 2000\right)\\ M_n=160\times {10}^{-7}{f}_0{D}_m^3\left(vn<2000\right)\\ M_z=f_1{P}_1{D}_m\end{array}$ (13)

式中：M为总摩擦力矩，N·mm；M_n 为与转速有关的粘性摩擦力矩，N·mm；f₀为与润滑方式有关的系数；v为润滑剂运动粘度，mm²/s；D_m为轴承平均直径，mm；M_z为与载荷有关的载荷摩擦力矩，N·mm；f₁为与轴承类型、载荷有关的系数，P₁为轴承所受载荷，N。

Figure 8. Geometric model for calculating heat transfer of motorized spindle

图8. 电主轴传热计算几何模型

电机生热量包括机械损耗、电损耗、磁损耗，可由下式计算 [14] ：

$Q_m=\left(1-\eta \right)P{P}_r$ (14)

式中：η为电机的机械效率；P是电机的比例负载；P_r为电机额定功率，W。因此，可以根据电主轴的机械效率与比例负载之间的关系来估算内置电机的最大发热量。

接触热阻初始值R可由半经验公式得出：

$\frac{r}{Rk}=105.6{\left(\frac{P}{E}\right)}^{0.68}$ (15)

式中：R为初始接触热阻，m²/(K·W)；r为粗糙度均方根，Ra；k为调谐平均导热系数，W/(m·K)；P为结合面所受压力，MPa；E为材料弹性模量，GPa。

4.2.2. 边界条件优化

由于上述热边界条件随电主轴组件热变形及工作状态而变化，为获得更贴近实际的电主轴热特性，需要对上述热边界条件进行实时优化，接触热阻的优化和对流换热系数的优化由式(6)和(11)通过关键测温点的实测及仿真温度进行优化。内热源生热量可通过经验公式进行优化 [15] ：

$Q^{\prime }=\frac{t_m-t_a}{t_f-t_a}Q$ (16)

式中：Q'为修正后的生热量，W；t_m为关键测温点实测温度，℃；t_f为有限元仿真温度，℃；t_a为环境温度，℃。

表1为电主轴部分热边界条件的初始值及优化值，其中接触热阻1为前轴承与轴承套、接触热阻2为前轴承套与冷却套、接触热阻3为后轴承与轴承套、接触热阻4为后轴承套与冷却套、接触热阻5为电机定子与冷却套。

Table 1. Initial and optimized values of thermal boundary conditions

表1. 部分热边界条件的初始值及优化值

4.3. 实验结果分析

将3.2节计算所得电主轴初始热边界条件输入热特性数字孪生系统，开启电主轴及热特性数字孪生系统进行电主轴热特性数字孪生实验。图9为电主轴运行2400 s时的孪生温度场云图，由图9可知，电主轴最大温升在转子处，为14.72℃，其次为定子和前后轴承；最小温升在主轴前端及对应前后轴承的冷却套，为5.63℃。由于转子及芯轴散热困难，导致其温升最大，因此，在电主轴热优化设计时应考虑改善转子、芯轴及轴承内圈的散热。另外，图9还显示，电主轴的温度场非常不均匀，这样会引起电主轴组件的不均匀热变形，因此，还需考虑冷却套的优化设计，以减小定子处的温升，使电主轴温度场更加均匀。

Figure 9. Temperature field of themotorized spindle

图9. 电主轴温度云图

图10为电主轴运行2400 s时的热变形云图。由图10可知，电主轴最大热变形发生在主轴轴端，为11.86 μm，由于主轴装夹刀具进行加工，该热变形会严重影响其加工精度，因此，可通过变形约束或误差补偿的方式减小主轴热变形对加工精度的影响。

Figure 10. Thermal deformation of the motorized spindle

图10. 电主轴热变形云图

图11为电主轴关键测温点实测温度与数字孪生仿真温度的对比图；由图11可知，电主轴前轴承的温升比后轴承快，主要是由于前轴承类型、载荷、润滑等导致前轴承生热量大。冷却套处由于冷却水的对流换热温度明显降低。图11还显示，电主轴组件在800 s前温升速度大，这主要是由于电主轴冷态时，组件与热源的温差大，导致初始热传导速度快，随着主轴组件温度升高，其温差逐渐减小，800 s后各组件的温升速度减小。由各测温点实测温度与仿真温度对比可知，数字孪生系统的跟随性好，温度仿真精度能达98%。

Figure 11. Temperature rise diagram of key position of motorized spindle

图11. 电主轴关键位置温升图

图12为电主轴轴端径向与轴向热变形实测值与仿真值对比图；由图12可知，电主轴轴端热变形在2400 s后趋于稳定，轴端轴向最大热变形为12.36 μm，径向最大热变形为11.24 μm。热变形仿真精度达96%，说明热特性数字孪生系统能较精确的反应电主轴热特性。

Figure 12. Radial and axial thermal deformation of the shaft end

图12. 轴端径向与轴向热变形

5. 结论

本文提出一种基于最小二乘优化边界的电主轴热特性数字孪生方法，建立电主轴热边界条件目标函数，设计开发基于C++的热特性数字孪生系统，实现电主轴热特性数字孪生。实验结果表明，电主轴温度场孪生值精度误差控制在2%之内，热变形孪生精度误差控制在4%之内。为电主轴热特性精确分析及控制提供依据。本文得出主要结论如下：

1) 温度是电主轴热特性的主要特征参数之一，采用关键测温点的实测及仿真温度可实现热特性有限元仿真所施加边界条件的实时优化，进而提高电主轴热特性分析精度。

2) 基于热边界的实时监测、优化及映射，可实现电主轴物理实体热特性在虚拟实体的同步映射，为电主轴热特性实时监测及热优化设计提供依据，同时可为误差补偿提供热变形数据。

3) 由于电主轴内热源生热量不同，传热路径不同，均匀分布流道的冷却套难以实现电主轴各组件冷却量的匹配，导致电主轴温度场不均匀。