CN113051686B - 倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法及热误差模型建模方法 - Google Patents

Abstract

本发明公开了一种倾斜工作条件下的主轴系统热‑结构耦合特性模型创建方法，包括如下步骤：1)创建轴倾斜工作条件下主轴系统的热‑结构耦合特性仿真模型，并将倾斜工作条件下的主轴系统的热源和热边界条件应用于所述热‑结构耦合特性仿真模型；2)根据所述热‑结构耦合特性仿真模型对主轴系统进行瞬态分析；3)比较相邻两个迭代子步骤中的主轴的温度，若满足收敛条件，则进入步骤4)；若不满足收敛条件，则进入步骤5)；4)结束分析，保存轴系统中所有温度节点的温度；5)更新主轴系统的热源和热边界条件；6)循环步骤2)和步骤3)，直至满足收敛条件。本发明还公开了一种倾斜工作条件下的主轴系统热误差建模方法。

Description

倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法 及热误差模型建模方法

技术领域

本发明属于机械误差分析技术领域，具体的为一种倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法及热误差模型建模方法。

背景技术

轴系统包括轴和分别位于轴两端的轴承、轴承箱和轴承套圈等。高速主轴系统的热误差与轴承的发热密切相关，轴承的倾斜会严重影响工作条件，因而需要考虑轴承倾斜状态下对主轴系统热误差的影响。

发明内容

有鉴于此，本发明的目的在于提供一种倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法及热误差模型建模方法，考虑倾斜工作条件对轴承的发热以及轴向和径向刚度的影响，进而得到主轴系统热-结构耦合特性模型，通过该主轴系统热-结构耦合特性模型建立了温度场，进而获得了主轴系统各部件的热变形，从而构建主轴系统热误差模型。

为达到上述目的，本发明提供如下技术方案：

本发明首先提出了一种倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，包括如下步骤：

1)创建主轴系统的热-结构耦合特性仿真模型，并将倾斜工作条件下的主轴系统的热源和热边界条件应用于所述热-结构耦合特性仿真模型；

倾斜工作条件下的主轴系统的热源包括：轴承的发热、电动机的发热；

倾斜工作条件下的主轴系统的热边界条件包括：强制对流换热系数、自由对流换热系数、轴承内圈与轴颈之间的接触热阻、轴承外圈与轴承座之间的接触热阻以及轴承滚珠与环之间的接触热阻；

倾斜工作条件下的主轴系统力学边界条件包括：轴承的径向接触刚度、轴承的轴向接触刚度。

2)根据所述热-结构耦合特性仿真模型对主轴系统进行瞬态分析；

3)比较相邻两个迭代子步骤中的主轴的温度，若满足收敛条件，则进入步骤4)；若不满足收敛条件，则进入步骤5)；所述收敛条件为：

T_i-T_i-1＜ΔT

其中，T_i表示第i次迭代子步骤中得到的主轴温度；ΔT表示收敛温差；

4)结束分析，保存轴系统中所有温度节点的温度；

5)更新主轴系统的热源和热边界条件；

6)循环步骤2)和步骤3)，直至满足收敛条件。

进一步，轴承的发热Q_f为：

Q_f＝Q_l+Q_v+Q_s+Q_gy

其中，Q_l、Q_v、Q_s和Q_gy分别为外部负载、润滑剂的粘滞摩擦、球的旋转运动和陀螺运动产生的热量；

Q_l＝M_l·ω

F_β＝max(0.9F_a/tanα-0.1F_r,F_r)

其中，M_l表示由外部载荷引起的摩擦扭矩；ω表示轴承的转动角速度；f₁表示与轴承结构和载荷油管的系数，F_s表示等效静载荷；C_s表示基本额定静载荷；d_m表示轴承的节圆直径；F_β表示等效载荷；F_a表示轴承的轴向载荷；F_r表示轴承的径向载荷；α表示实际接触角；

Q_v＝M_v·ω

其中，M_v表示由润滑剂的粘滞摩擦引起的摩擦扭矩；n表示内圈的转速；f_o表示与轴承类型和润滑条件有关的常数；v_o表示在工作温度下润滑剂的运动粘度；

Q_s＝Z·M_si·ω_si+Z·M_so·ω_so

其中，M_s表示由滚珠的旋转运动引起的旋转摩擦扭矩；η_oil表示滚珠与槽之间的摩擦系数；Q_i/o表示施加在内沟道或外沟道上的负载；a_i/o表示接触椭圆的半长轴；ε_i/o表示第二类的完整椭圆积分；M_si表示旋转运动引起的滚珠与内沟道之间的旋转摩擦扭矩；ω_si表示滚珠相对于内沟道的转动角速度；M_so表示旋转运动引起的滚珠与外沟道之间的旋转摩擦扭矩；ω_so滚珠相对于外沟道的转动角速度；Z表示滚珠的数量；

Q_gy＝M_gy·ω

其中，M_gy表示由陀螺力矩引起的摩擦转矩；β表示滚珠的运动姿势；ω_c表示方位角

处的公转角速度；ω_b表示每个滚珠的转动角速度；ρ表示滚动体密度；D_b表示滚珠的直径。

进一步，轴承的轴向接触刚度K_a为：

其中，K_aj表示第j^th个滚珠的轴向接触刚度；K_aij表示第j^th个滚珠/内沟道界面之间的轴向接触刚度；K_aoj为第j^th个滚珠/外沟道界面之间的轴向接触刚度；Z表示滚珠的数量；

轴承的径向接触刚度K_r为：

其中，K_rj表示第j^th个滚珠的径向接触刚度；K_rij表示第j^th个滚珠/内沟道界面之间的径向接触刚度；K_roj为第j^th个滚珠/外沟道界面之间的径向接触刚度；

为第j^th个滚珠的方位角；

且：

其中，Γ表示第一类和第二类的完整椭圆积分；E'表示等效弹性模量；κ表示椭圆偏心率参数；∑ρ表示在正交平面上的接触点处的接触体的曲率之和；ε表示接触变形；α_ij表示内接触角；α_oj表示外接触角；Q_ij和Q_oj表示内沟道和外沟道分别施加到j^th滚珠的法向接触载荷；且：

其中，λ_ij和λ_oj为系数；D_b表示滚珠直径；M_gj表示陀螺力矩；F_cj表示j^th滚珠的离心力；且：

其中，m表示滚珠质量；d_m表示轴承的节圆直径；ω_c和ω_b表示滚珠的公转和自转角速度；

α_o表示初始接触角；ω_cj表示第j个滚动体公转角速度；ω表示轴承内圈的转动角速度；

且：

其中，X_1j表示内滚道沟曲率中心到球中心最终位置的轴向距离；X_2j表示内滚道沟曲率中心到球中心最终位置的径向距离；A_1j和A_2j分别表示内沟道、外沟道曲率中心之间的轴向和径向距离；f_i表示内沟道曲率系数；δ_ij表示接触载荷作用下的正常接触变形；f_o表示外沟道曲率系数；δ_oj表示接触载荷作用下的正常接触变形；其中：

(A_1j-X_1j)²+(A_2j-X_2j)²-[(f_o-0.5)D_b+δ_oj]²＝0

其中，A_o表示内外沟道之间的初始距离；R_o表示外沟道曲率中心的走线半径；

表示方位角；α^o表示静载荷作用下内圈与外圈的接触角；δ_a和δ_r分别为内圈和外圈之间的相对轴向位移和径向位移；θ表示轴承的倾斜角，且

其中ξ和d分别表示间隔环的平行度误差和外径。

进一步，所述步骤5)中，轴承发热量的更新方法为：

更新润滑剂粘度v：

v＝32×e^{-0.0242×(T-40)}

其中，T为从热-结构耦合特性仿真模型仿真中提取的温度；

更新主轴的轴向变形，得到更新后的主轴轴向变形δ_a为：

δ_a＝δ_preload±Δδ_y

其中，δ_a为更新后的主轴的轴向变形；δ_preload表示初始预紧力作用下的主轴的轴向变形；Δδ_y表示通过热-结构耦合特性仿真模型获得主轴的轴向热伸长；

根据轴向变形与预紧力之间的相关性，采用牛顿-拉夫森迭代算法计算实际接触角α：

δ_a＝BD_b sin(α-α_o)/cosα

其中，α_o为初始接触角；B为曲率半径和；D_b表示滚珠的直径；

通过将实际接触角α获取实际预紧力P_a：

其中，K_n＝K·D_b ^0.5/B^1.5表示轴向变形常数，K表示接触刚度；Z表示滚珠的数量；

接触载荷与接触变形的关系为：

对该等式两端求导，即可将接触面积的接触刚度表示为

其中，Γ表示第一类和第二类的完整椭圆积分；E'表示等效弹性模量；κ表示椭圆偏心率参数；∑ρ表示在正交平面上的接触点处的接触体的曲率之和；ε表示接触变性；Q表示接触载荷；

实现对接触刚度的更新。

进一步，电动机的发热Q为：

Q＝P_in(1-η)

其中，P_in为电机输入功率；η为电机效率。

进一步，自然对流换热系数和强制对流换热系数表示为：

其中，h为对流换热系数；Nu为努塞尔数；L为水力直径；λ为流体的导热系数。

进一步，轴承内圈与轴颈之间的接触热阻表示为：

其中，k₁和k₂分别为两种接触固体材料的导热系数；k₃为间隙中介质的导热系数；

为无因次接触面积；A_n表示名义接触面积；H_s为接触热导率；且：

其中，h_cp表示塑性变形引起的接触热导率；h_ep表示弹性变形引起的接触热导率；h_cep弹塑性变形引起的接触热导率；n(a)表示面积介于a和a+da之间的粗糙体的数量；a为接触面积；a_c1为接触点开始发生塑性变形时的临界接触面积；a_c2为完全塑性变形时的临界塑性变形；

其中，R_s为表示接触半径；D为分形维数；f(a)为模板函数；a_L为最大接触面积；

所述步骤5)中，轴承内圈与轴颈之间的接触热阻的更新方法为：

为了保持平衡，由接触粗糙物的变形产生的总接触载荷F应该等于外部载荷F_ext，即：

F＝F_ext＝P·A_a

其中，A_a为标称表面积；P为接触压力；且接触压力P是由轴承内圈与轴颈之间的过盈配合I产生，并表示为：

轴颈与轴承内圈之间的过盈配合随主轴系统的热变形而变化，并且偏离初始过盈量I；动态过盈配合表示为：

I'＝I+ΔI

其中，ΔI标致轴颈和轴承内圈的径向热变形之间的差值，ΔI可通过热-结构耦合特性仿真模型的每个迭代子步骤获取，从而求取得到更新的最大接触面积a_L，以更新轴承内圈与轴颈之间的接触热阻。

进一步，轴承外圈与轴承座之间的接触热阻R为：

其中，δ_ring表示外圈厚度；λ_air表示空气的热导率；λ_ring表示轴承外圈的热导率；R_s为表示接触半径；δ_gap表示间隙，且：

δ_gap＝δ_gap0-(T_ring-T_housing)·α·r_housing

其中，δ_gap0表示初始间隙；T_ring表示外圈的平均温度；T_housing表示轴承座的平均温度；α表示热膨胀系数；r_housing表示接触区域的半径。

进一步，轴承滚珠与环之间的接触热阻为：

其中，λ_b表示滚珠的热导率；λ_ring表示环的热导率；ψ表示与接触表面的尺寸有关的几何参数；a为赫兹接触区域的长半轴。

本发明还提出了一种倾斜工作条件下的主轴系统热误差建模方法，主轴系统的热误差模型为：

δ(t)＝δ₀+(δ_t-δ₀)(1-e^-t/τ)

其中，δ₀表示主轴系统的初始状态热误差，δ_t表示主轴系统在热平衡状态下的热误差；τ表示热误差平衡过程的时间常数；

采用如权利要求1-9任一项所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，在迭代子步骤中计算主轴系统的热误差，且当

时，主轴系统达到热平衡状态，此时的时间t_i为热平衡时间，则

此时的热误差E_i为主轴系统的热平衡误差δ_t；其中，ΔE_max＝E_max-E₀，E_max和E₀分别表示主轴系统的最大热误差和初始热误差；f_i表示t_i时刻主轴系统的热误差变化相对于主轴系统的最大热误差变化的百分比；ΔE_i表示t_i时刻的热误差变化，且ΔE_i＝E_i-E₀。

本发明的有益效果在于：

本发明倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法及热误差模型建模方法，需要考虑轴系统倾斜状态下轴承的变形，导致的轴承发热的变化、电动机发热的变化、强制对流换热系数、自由对流换热系数、轴承内圈与轴颈之间的接触热阻、轴承外圈与轴承座之间的接触热阻以及轴承滚珠与环之间的接触热阻的变化，以及主轴轴承轴向接触刚度、轴承的径向接触刚度的变化；将倾斜角引入轴承的准静态力学模型中，以考虑倾斜工作条件对轴承的发热以及轴向和径向刚度的影响，然后计算了倾斜状态下散热边界条件，包括对流系数和轴承接头的热接触电阻(TCR)，将上述热源和热边界条件应用于热-结构耦合特性仿真模型，从而用于对主轴系统进行瞬态分析，在满足收敛条件时，主轴系统达到热平衡状态，从而可以确定主轴系统热误差模型的热误差和时间常数。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为倾斜工作条件下的主轴系统的示意图；

图2为倾斜工作条件下内沟道、外沟道曲率中心和滚珠中心的相对位置示意图；

图3为倾斜工作条件下倾斜角度的示意图；

图4为施加在滚珠上的力的示意图；

图5为倾斜工作条件下轴承部件的受力示意图；

图6为倾斜工作条件对热量产生的影响的曲线图；

图7为倾斜工作条件下每个接触椭圆处的热量分布曲线图；

图8为倾斜工作条件下接触负载分布曲线图；

图9为轴承的刚度示意图；

图10为不同倾斜工作条件下轴承的轴承刚度曲线图；

图11为摩擦转矩与温度的关系曲线图；

图12为预紧力降低对热量产生和总摩擦扭矩的影响曲线图；

图13为倾斜工作条件下每个接触区域的径向刚度分布曲线图；

图14为倾斜工作条件下接触角分布曲线图；

图15为轴承接头的总热接触电阻示意图；

图16为粗糙的表面形态示意图；

图17为粗糙的表面形态图；

图18为双对数坐标中的功率谱函数曲线拟合图；

图19为两个圆柱面之间的接触示意图；

图20为TCR和表面接触系数的曲线图；图20(a)为TCR与干扰之间的关系曲线图；图20(b)为表面接触系数曲线图；

图21为本实施例主轴系统的实验装置的示意图；

图22为倾斜工作条件下主轴系统温度场；

图23为倾斜工作条件下前轴承组件的温度场；

图24为倾斜工作条件下轴承部件的温度曲线图；

图25(a)为倾斜工作条件下后轴承外圈的温度曲线图；图25(b)为倾斜工作条件下前轴承外圈的温度曲线图；

图26为倾斜工作条件下热变形的仿真图；

图27为倾斜工作条件下轴承的热应力仿真图；

图28为倾斜工作条件下主轴系统分别在XYZ方向上的热误差的曲线图；

图29为倾斜工作条件下热误差平衡状态准则曲线图；

图30为倾斜工作条件下主轴系统的稳态热误差；

图31为倾斜工作条件下主轴系统的转速曲线图；

图32为倾斜工作条件下主轴系统的热误差补偿结果曲线图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。

本实施例倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，包括如下步骤：

1)创建主轴系统的热-结构耦合特性仿真模型，并将倾斜工作条件下的主轴系统的热源和热边界条件应用于所述热-结构耦合特性仿真模型；

倾斜工作条件下的主轴系统的热源包括：轴承的发热、电动机的发热；

倾斜工作条件下的主轴系统的热边界条件包括：强制对流换热系数、自由对流换热系数、轴承内圈与轴颈之间的接触热阻、轴承外圈与轴承座之间的接触热阻以及轴承滚珠与环之间的接触热阻；

倾斜状态条件下力学边界条件包括：轴承的轴向接触刚度、轴承的径向接触刚度。

2)根据所述热-结构耦合特性仿真模型对主轴系统进行瞬态分析；

3)比较相邻两个迭代子步骤中的主轴的温度，若满足收敛条件，则进入步骤4)；若不满足收敛条件，则进入步骤5)；所述收敛条件为：

T_i-T_i-1＜ΔT

其中，T_i表示第i次迭代子步骤中得到的主轴温度；ΔT表示收敛温差，本实施例的ΔT＝0.1℃；

4)结束分析，保存轴系统中所有温度节点的温度；

5)更新主轴系统的热源和热边界条件；

6)循环步骤2)和步骤3)，直至满足收敛条件。

具体的，对于高速主轴系统，冷却系统对于散发电机和轴承的热量很重要，有关冷却系统的信息如下：油冷却器是DAIKIN制造的AKZ328-C-D210，使用的冷却介质为R410，R410的入口温度为15℃,循环冷却油的流量为24/28.8L/min,电源电压为三相电源，AC200V50/60Hz,制冷量为2.8/3.2KW,工作温度范围为

冷却套的矩形凹槽的长度和宽度分别为18.5mm和8.25mm。材料特性对于主轴系统的热-结构耦合特性分析很重要，因为这些参数直接确定温度升高和热变形。在主轴系统的闭环迭代热分析中，有必要根据表1设置主轴系统主要零件的初始材料属性。

表1主要零件的材料特性

1、主轴系统的热源

主轴系统中的主要热源是前后轴承和内置电动机。

1.1、轴承的发热和刚度

轴承发热的计算方法是准静态力学分析。由于倾斜的工作条件，在接触过程中在轴承外圈和间隔圈之间存在倾斜间隙。此外，在径向载荷F_r和轴向载荷F_a的共同作用下，由于外圈/轴承箱的间隙配合，将产生轴承内圈和外圈之间的相对轴向位移δ_a、径向位移δ_r和角位移θ以及轴承的内部游隙，如图1所示。

内沟道和外沟道的曲率中心之间的初始距离为A_o，如图2所示。球/内圈的接触角不等于球/外圈的接触角，轴承滚珠的中心当轴向和径向载荷施加到轴承上时，也会发生变化。内沟道的固定曲率中心与球中心的最终位置之间的距离Δ_ij表示为：

Δ_ij＝(f_i-0.5)D_b+δ_ij (4-a)

其中，f_i表示内沟道曲率系数；D_b表示滚珠直径；δ_ij表示接触载荷作用下的正常接触变形。

外沟道的固定曲率中心与滚珠中心的最终位置之间的距离Δ_oj表示为：

Δ_oj＝(f_o-0.5)D_b+δ_oj (4-b)

其中，f_o表示外沟道曲率系数；δ_oj表示接触载荷作用下的正常接触变形。

内圈和外圈之间的相对轴向位移、径向位移分别为δ_a、δ_r，轴倾斜工作条件下引起的内圈和外圈之间的相对角度位移θ。即角位移θ是由倾斜的工作条件引起的。然后，内沟道、外沟道曲率中心之间的轴向和径向距离表示为：

其中，A_1j表示内沟道、外沟道曲率中心之间的轴向距离；A_2j表示内沟道、外沟道曲率中心之间的径向距离；A_o表示内外沟道之间的初始距离；R_o表示外沟道曲率中心的走线半径；

表示方位角；α^o表示静载荷作用下内圈与外圈的接触角。

在研究中，倾斜工作条件的影响等效于间隔环的平行度误差，以及根据图3，倾斜角

其中ξ和d分别表示间隔环的平行度误差和外径。

由于倾斜工作条件改变了内沟道、外沟道曲率中心之间的轴向和径向距离A_1j和A_2j，导致轴承的几何相容方程发生了变化，根据勾股定理和图3所示的几何关系，将(4-c)所示倾斜工作条件下的内沟道、外沟道曲率中心之间的轴向和径向距离A_1j和A_2j引入，获得以下关系：

其中，X_1j表示内滚道沟曲率中心到球中心最终位置的轴向距离；X_2j表示内滚道沟曲率中心到球中心最终位置的径向距离；

而且，内接触角和外接触角满足以下关系：

其中，α_ij表示内接触角；α_oj表示外接触角。

如等式(4-c)所示，倾斜工况对A_1j和A_2j有影响，然后影响几何变形相容方程。最后，影响内沟道和外沟道的变形δ_ij和δ_oj。滚珠受离心力F_cj、陀螺力矩M_gj以及法向接触载荷Q_ij和Q_oj的受到倾斜工况的影响，导致与不倾斜条件下的结果有差异。如图4所示。对于每个滚珠，应平衡在水平和垂直方向上施加的载荷。倾斜工作条件下的几何变形相容方程与正常情况下不同。

其中，Q_ij和Q_oj表示内沟道和外沟道分别施加到j^th滚珠的法向接触载荷；λ_ij和λ_oj为系数，且在低速范围内，λ_ij＝λ_oj＝1，在高速范围内，λ_ij＝0andλ_oj＝2；M_gj表示陀螺力矩；

j^th滚珠的离心力F_cj和陀螺力矩M_gj表示为：

其中，m表示滚珠质量；d_m表示轴承的节圆直径；ω_c和ω_b表示滚珠的公转和自转角速度；

α_o表示初始接触角；ω_cj表示第j个滚动体公转角速度；ω表示轴承内圈的转动角速度；

通过结合式(4)、等式(5)、式(6)和等式(7)，如果确定X_1j,X_2j,δ_ij和δ_oj的值，则δ_a和δ_r可以求解。为了确定δ_a和δ_r，应建立整个轴承的平衡方程。施加在轴承组件上的力如图5所示，获得轴承外圈的力平衡方程。倾斜工作条件下的轴承外圈平衡方程与正常情况下不同。

其中，

表示j^th滚珠的方位角；F_a和F_r分别表示轴承的轴向和径向载荷；Z表示滚珠的数量。

然后，可得内、外轴承套圈的倾斜力矩，倾斜工作条件下的轴承外圈力矩平衡方程与正常情况下不同。

其中，M表示轴承的弯矩；R_i表示内滚道沟曲率中心轨迹半径；r_i滚珠班级；

表示第j个球所在的位置角。

可以看出，倾斜角度导致内圈和外圈之间的相对倾斜，从而导致载荷分布不均匀。此外，还会产生一个额外的力矩，然后该力矩通过轴承内圈传递到主轴系统。

1.1.1、轴承发热

轴承倾斜工作条件导致用于计算轴承发热的动力学和运动学参数发生了变化，从而影响轴承生热量。轴承的发热是由外力、润滑剂的粘滞摩擦、滚珠的旋转和陀螺运动四种摩擦转矩引起的，因此轴承的发热表示为：

Q_f＝Q_l+Q_v+Q_s+Q_gy (10)

其中，Q_l、Q_v、Q_s和Q_gy分别为外部负载、润滑剂的粘滞摩擦、球的旋转运动和陀螺运动产生的热量。

其中，由外部载荷引起的摩擦扭矩及由外部负载引起的热量表达为：

F_β＝max(0.9F_a/tanα-0.1F_r,F_r) (12)

Q_l＝M_l·ω (13)

其中，M_l表示由外部载荷引起的摩擦扭矩；ω表示轴承的转动角速度；f₁表示与轴承结构和载荷油管的系数，F_s表示等效静载荷；C_s表示基本额定静载荷；d_m表示轴承的节圆直径；F_β表示等效载荷；F_a表示轴承的轴向载荷；F_r表示轴承的径向载荷；α表示实际接触角；

由于润滑剂的粘滞摩擦而产生的热量及由粘滞摩擦产生的热量表示为：

Q_v＝M_v·ω (15)

其中，M_v表示由润滑剂的粘滞摩擦引起的摩擦扭矩；n表示内圈的转速；f_o表示与轴承类型和润滑条件有关的常数；v_o表示在工作温度下润滑剂的运动粘度；

由滚珠的运动引起的旋转摩擦扭矩及由所有球的旋转运动引起的总热量表示为：

Q_s＝Z·M_si·ω_si+Z·M_so·ω_so (17)

其中，M_s表示由滚珠的旋转运动引起的旋转摩擦扭矩；η_oil表示滚珠与槽之间的摩擦系数；Q_i/o表示施加在内沟道或外沟道上的负载；a_i/o表示接触椭圆的半长轴；ε_i/o表示第二类的完整椭圆积分；M_si表示旋转运动引起的滚珠与内沟道之间的旋转摩擦扭矩；ω_si表示滚珠相对于内沟道的转动角速度；M_so表示旋转运动引起的滚珠与外沟道之间的旋转摩擦扭矩；ω_so滚珠相对于外沟道的转动角速度；Z表示滚珠的数量；

陀螺力矩引起的摩擦转矩及由陀螺力矩产生的热量表示为：

Q_gy＝M_gy·ω (19)

其中，M_gy表示由陀螺力矩引起的摩擦转矩；β表示滚珠的运动姿势；ω_c表示方位角

处的公转角速度；ω_b表示每个滚珠的转动角速度；ρ表示滚动体材料密度；D_b表示滚珠的直径。

在本实施例的主轴系统中，使用了一对型号为FAG HS7014-C-T-P4S的角接触轴承。表2中列出了轴承的参数。为了研究倾斜工作条件对轴承动态特性和发热的影响，轴承的具体工作条件设置如下，径向和轴向载荷分别为300N和1000N。如图6所示，通过准静态力学模型计算轴承在不同倾斜角度下的发热量，可以看出该热量随着倾斜工作条件的增加而增加。当倾斜误差为40μm时，与没有倾斜角的情况相比，热量产生增加了29.23％。而且，热量的产生随着转速的增加而显着增加，从而导致轴承的温度明显升高。

表2 FAG HS7014-C-T-P4S轴承参数

图7显示了在滚珠/轴承圈的每个接触椭圆处产生的热量。倾斜工作条件会导致滚珠/内圈和滚珠/外圈接触区域的热量产生不均匀，这与图8中所示的接触负载分布一致。即，当接触椭圆时，在每个接触椭圆中的热量都会增加。热量和接触负荷都随着倾斜角而增加。当倾斜角为0°时，在每个方位角处的滚珠的发热和接触载荷是均匀的。但是，当存在倾斜角时，在每个方位角处的滚珠的发热和接触载荷是不均匀的。当倾斜角度为0时，每个方位角度的滚珠的接触载荷都是均匀的，随着倾斜角的增加，接触载荷分布将越来越不均匀。而且，倾斜角越大，最大接触载荷越大。当倾斜角度为40μm时，最大接触负载比没有倾斜角度时增加了99.13％。结论是，倾斜角导致接触负载显着增加。此外，接触负载取决于方位角，并且接触负载在两个方位角

处达到最大值，而接触负载在方位角

处达到最小值。主要原因是在两个极限方位角

处的滚珠承受着主要的径向载荷。同时，球在其他两个极限方位角

处承受的径向载荷负担大大降低。因此，根据等式(17)、(18)和(19)，热量的产生取决于方位角，热量在两个极限方位角

处达到最大值，而在其他两个方位角

处达到最小值。在交变接触载荷的作用下，轴承更有可能加剧轴承元件的疲劳损伤。因此，必须减小外圈的倾角，以减少接触载荷和每个滚动体的发热，并提高轴承的使用寿命。

1.1.2、轴向和径向刚度

轴向和径向刚度非常重要，因为它们决定了主轴系统的热变形，而间轴承倾斜工作条件也会影响轴承的刚度。接触载荷和接触变形之间的关系为：

对该等式两端求导，即可将接触面积的接触刚度K表示为：

其中，δ表示接触变形；Γ表示第一类和第二类的完整椭圆积分；E'表示等效弹性模量；κ表示椭圆偏心率参数；∑ρ表示在正交平面上的接触点处的接触体的曲率之和；ε表示接触变形；Q表示接触载荷，具体的，可以表示内沟道和外沟道分别施加到j^th滚珠的法向接触载荷Q_ij和Q_oj；

当已知接触载荷和接触角时，将内沟道和外沟道在径向和轴向的法向接触刚度的分量表示为：

根据弹簧的串联和并联原理，如图9所示，得出j^th滚珠/内沟道界面的径向刚度K_rij和轴向接触刚度K_aij：

则轴承的总径向刚度K_r和轴向刚度K_a表示为：

其中，K_aj表示第j^th个滚珠的轴向接触刚度；K_aij表示第j^th个滚珠/内沟道界面之间的轴向接触刚度；K_aoj为第j^th个滚珠/外沟道界面之间的轴向接触刚度；Z表示滚珠的数量；K_rj表示第j^th个滚珠的径向接触刚度；K_rij表示第j^th个滚珠/内沟道界面之间的径向接触刚度；K_roj为第j^th个滚珠/外沟道界面之间的径向接触刚度；

为第j^th个滚珠的方位角。

从以上分析可以看出，球的每个接触椭圆的径向刚度和整个轴承的径向刚度是通过改变方位角

来计算的。轴承的径向刚度K_r和轴向刚度K_a是在不同倾斜度下获得的，如图10所示，倾斜工作条件导致径向刚度的不均匀分布，并且径向刚度在与倾斜方向平行的方向上呈现出增加的趋势。然而，径向刚度在垂直于倾斜方向的方向上减小。而且，轴承的轴向刚度随倾斜角而减小。因此，通过减小倾斜角可以显着提高轴向刚度。即倾斜状态下轴向刚度减小，降低了轴系统抵抗外力的能力，同时，温度升高后也更容易发生轴向变形。

1.1.3、轴承发热的更新

然后根据粘度温度效应和轴承的热诱导预紧力更新轴承的热量生成，以实现热量和刚度的精确建模。温度升高导致润滑剂粘度降低，并且由于粘度-温度效应而直接影响球轴承的发热。

v＝32×e^{-0.0242×(T-40)} (25)

其中，T为从热-结构耦合特性仿真模型仿真中提取的温度。

如图11所示，将润滑剂粘度变化引入模型以更新轴承的发热。粘度随着温度升高而降低，总摩擦扭矩逐渐升高，然后在温度为70℃时开始降低。更新润滑剂的粘度，从而能够更新轴承发热。

主要的预紧力变化归因于轴的热伸长。更新后的轴向变形表示为：

δ_a＝δ_preload±Δδ_y (26)

其中，δ_a为更新后的主轴的轴向变形；δ_preload表示初始预紧力作用下的主轴的轴向变形；Δδ_y表示通过热-结构耦合特性仿真模型获得主轴的轴向热伸长；

根据轴向变形与预紧力之间的相关性，采用牛顿-拉夫森迭代算法计算实际接触角α：

δ_a＝BD_bsin(α-α_o)/cosα (27)

其中，α_o为初始接触角；B为曲率半径和；D_b表示滚珠的直径；

通过将实际接触角α获取实际预紧力P_a：

其中，K_n＝K·D_b ^0.5/B^1.5表示轴向变形常数，K表示接触刚度；Z表示滚珠的数量；

更新后的预紧力被引入轴承热模型，然后更新轴承发热，如图12所示。初始预紧力为600N，当预紧力减小400N时，总摩擦扭矩仅减小约30N.mm。因此，热引发的预紧力对轴承发热的影响远小于润滑剂粘度对轴承发热的影响。

如图13所示，获得了滚珠/环的每个接触区域的径向刚度分布。根据等式(24)，倾斜角度导致每个球上的接触载荷分布不均匀，因此轴承刚度的分布不均匀。此外，径向刚度在两个方位角

处达到最大值，而接触载荷在方位角

处达到最小值。原因是两个极端方位角

处的滚珠承受了主要的径向载荷，而另外两个极端方位角

处的滚珠所承受的径向载荷大大减小了，从而导致在两个极端的方位角

和其他两个极端的方位角

的最小接触载荷。然后，径向刚度在两个方位角

处达到最大值，而接触载荷在方位角

处达到最小值。

如图14所示，也可以得到外沟道与滚珠之间的接触角分布。接触角在两个方位角

处达到最大值，而在其他两个方位角

达到最小值。即，接触角也会交替变化，验证了考虑倾斜工作条件影响的准静态力学分析模型的有效性。此外，倾斜误差越大，在不同位置角下的接触角越大。

1.2、电机发热

电机产生的热量是重要的热源，效率分析方法用于计算热量的产生。

Q＝P_in(1-η) (28)

其中，P_in为电机输入功率；η为电机效率。

转子和定子之间的发热分布由同步频率f_sync和滑差比f_slip确定。

其中，Q_r和Q_s分别表示转子和定子的产生的热量。

2、对流换热

对流系数表示为：

其中，h为对流换热系数；Nu为努塞尔数；L为水力直径；λ为流体的导热系数。

固定表面的对流状态属于自由对流，标准判据方程表示为：

其中，C、m和n均为与表面取向有关的常数；Gr和Pr分别为Grashof数和Prandtl数；g表示重力加速度；β₂表示流体的膨胀系数；v表示运动粘度；Δt表示主轴系统各组件与空气之间的温差，且Δt＝T-T_a，T表示主轴部件的温度，T_a表示环境温度；L表示水力半径。

运动表面的对流状态属于强制对流，Nu表示为：

Nu＝0.133Re^2/3Pr^1/3 (32)

其中，Re表示雷诺数，且Re＝ω×d/v,ω和d分别角速度和直径。

然后，计算自由对流系数和强迫对流系数，如表3所示。

表3对流系数

3、轴承界面的接触热阻

轴承界面的TCR直接决定主轴系统的温度场，因为轴承界面的TCR限制了从轴承流到轴颈和轴承座的热量。因此，在热-结构耦合特性仿真建模中，应考虑轴承界面的TCR。

3.1、内圈/轴颈接口的TCR

轴承套圈、轴颈和轴承箱的粗糙表面由许多不同尺度的粗糙物组成，当外部载荷施加到接触件上时，接触界面的粗糙物会发生弹性、弹塑性和塑性变形。接触粗糙体的变形方式与临界弹性变形和塑性变形密切相关。

3.1.1、分形粗糙表面表征

(1)临界变形

当最大赫兹接触压力达到0.6H时，在接触粗糙体上开始发生塑性流动，然后将临界塑性变形表示为：

其中，H表示较软材料的硬度；k＝0.6，E表示弹性模量；G表示分形尺度常数；a表示接触面积；D表示分形维数；ω_c1表示临界塑性变形。

当ω＝ω_c1时，接触点开始发生塑性变形，然后相应的临界接触面积a_c1为：

当ω＝ω_c2时，最小接触压力达到H，并且纯塑性变形发生在接触粗糙体上。

然后将完全塑性变形的相应临界弹性接触面积a_c2表示为：

(2)单对接触粗糙体的TCR

根据图15获得单对接触粗糙体的TCR。

其中，k_s表示等效热导率，并满足2/k_s＝1/k₁+1/k₂，k₁和k₂分别表示两种接触材料的热导率；b和c分别接触面积半径和热流通道，它们之间的关系表示为：

然后，将单对接触粗糙体的TCR表示为：

其中，r表示单个粗糙接触区域a的接触半径；A_r ^*表示无因次接触面积，且A_r ^*＝A_r/A_n,A_r和A_n分别表示总实际接触面积和标称接触面积。然后，获得单对接触粗糙体的TCC，如下所示：

接触界面的大多数缝隙充满了普通的介质，包括空气/润滑剂，并且介质的热导率比金属的导热率小得多。因此，不能忽略媒体对TCR的影响。如图15所示，体热阻R_s与介质热阻R_f并联。然后，轴承接头的总TCR表示为：

介质热阻R_f与体热阻R_s之比表示为：

其中，

k₃表示空气/润滑剂的热导率。

总热阻R_c表示为：

其中，k₁和k₂分别为两种接触固体材料的导热系数；k₃为间隙中介质的导热系数；

为无因次接触面积；A_n表示名义接触面积；H_s为总的接触热导率；且：

(3)多重接触粗糙体的TCR

由单对接触粗糙体的弹性、弹塑性和塑性变形引起的TCR表示为：

其中，h_ce表示弹性接触热导率；h_cep表示弹塑性接触热导率；h_cp表示塑性接触热导率；k_s表示等效热导率；f(a)表示模板函数；

在给定的平均表面分离距离下，当分形维数D不等于1.5时，总TCC为：

其中，n(a)表示面积介于a和a+da之间的粗糙体的数量；a为接触面积；a_c1为接触点开始发生塑性变形时的临界接触面积；a_c2为完全塑性变形时的临界塑性变形；

若a_L＜a_c2，则：

若a_c2＜a_L＜a_c1，则：

若a_L＞a_c1，则：

热阻R_f表示为：

最后，总热阻R_c表示为：

从等式(38)和(39)中可以看出，TCR是以下参数的函数：固体的导热系数k₁和k₂，间隙中介质的导热系数k₃，分形维数D，分形粗糙度参数G，最大接触面积a_L，无因次接触面积

临界弹性变形a_c1和临界塑性变形a_c2。k₁,k₂,k₃的热导率是已知的。临界弹性变形a_c1和临界塑性变形a_c2通过公式(33-b)和(34-b)计算。粗糙表面的形貌由分形几何特征来表征，并且分形维数D和分形粗糙度参数G将通过结构函数方法来识别。建立多尺度接触力学模型，以计算最大接触面积a_L和无量纲接触面积

3.2、分形粗糙表面表征

(1)分形几何理论

表面形态的变化是具有固有分形特性的非平稳随机过程，即在各处连续，但不可区分。分形几何中的3DWeierstrass-Mandelbrot函数是最常用的函数，用于表征具有无序、非平稳、随机和自亲和行为的粗糙表面形态。

其中，x和表示每个粗糙的坐标；L表示样本长度；D_s表示3D形态的分形维数(2＜D_s＜3)；G表示分形粗糙度；γ表示缩放参数，它控制表面轮廓中的频率密度；φ_m,n表示随机阶段；n表示频率指数，且n_max＝int[loLg(_sL)/γl代表n的上限,L_s为截止长度。根据表面平整度和频率分布密度的考虑，γ＝1.5.如在单变量W-M函数中一样，将φ_m,n相视为均匀随机分布在[0,2π]中的独立随机变量。然后，使用随机数生成器将随机相位φ_m,n的值均匀分布在间隔[0,2π]中。3D W-M函数用于表征D_s＝2.67和G＝9.46×10^-6m时的粗糙表面形态，如图16所示。

(2)分形参数辨识

一般而言，二维表面形态用于识别分形维数D(D_s＝D+1)和缩放常数G。轴承套圈的表面形态是通过Olympus激光显微镜测量的，以识别分形参数。2D W-M函数的功率谱函数表示为：

等式(41-a)在两边都取对数，得到：

lgS(ω)＝2(D-1)lgG-lg2lnγ+(5-2D)lgω (41-b)

假设k＝2D-5以及b＝2(D-1)lgG-lg2lnγ，则等式(41-b)被改写为：

lgS(ω)＝klgω+b (41-c)

从式(41-c)可以看出，功率谱函数在双对数坐标中接近线性函数。然后，从直线的斜率获得2D分形维数D，并通过与纵轴的交点获得缩放常数G，即：

其中，γ表示尺度参数；ω表示频率。

如图17所示，通过Olympus激光显微镜OLS8000测量轴承套的粗糙表面形貌，并获得3D粗糙表面形貌。将数据作为输入来识别分形参数和功率谱函数。计算出X方向上的轮廓高度，如图18所示。X方向上的轮廓高度的功率谱函数几乎是一条直线，然后使用最小二乘拟合获得双对数中的红线坐标。最后，分形维数和比例常数由斜率和与纵轴的交点确定。

(3)表面接触系数

当忽略较小范围的细节并采用W-M函数的第一项时，在l的范围内，W-M函数如下：

其中，l表示接触椭圆半轴；

粗糙体的顶部曲率半径表示为：

两个接触粗糙表面的接触会产生大量不同大小的微接触点，这些微接触点共存于接触界面处。然后将这些接触点的数量和接触面积表示为：

其中，a_L表示接触点的最大接触面积。

然后，面积介于a和a+da之间的粗糙体的数量表示为

其中，N表示发生接触的微凸体数目。

轴承内圈和轴承箱之间的接触是两个粗糙曲面之间的接触。但是，等式(44)仅适用于两个粗糙表面之间接触的情况。因此，方程(44)不能用于计算两个一般粗糙曲面的接触界面上的粗糙分布数。粗糙的曲面形式会改变粗糙度分布，如图19所示。此外，接触参数(包括实际接触面积和总接触载荷)会受到粗糙度分布的影响。因此，两个粗糙曲面之间的TCR与两个粗糙曲面之间的TCR不同。当表面不是弯曲的，将构造传统的TCR模型。但是，内圈和轴颈的表面是圆柱表面，然后它们之间的接触就是圆柱表面之间的内部接触。表面形状的变化导致分形接触面积的变化，如图19所示。可以看出，当表面形状为圆柱面时，总接触面积会减小，这意味着传统的TCR模型不适合计算圆柱表面之间的接触参数。相反，它们仅适用于计算两个没有曲率的粗糙表面之间的接触的接触参数。为了计算上述接触参数和两个普通粗糙曲面之间的TCR，应修改粗糙度分布，然后应建立多尺度接触力学模型，以考虑由于表面形式的变化而考虑到的变化粗糙度分布。将表面接触系数引入等式(44-b)以修改粗糙分布。

其中R_s表示接触半径。

得到接触区域在a和a+da之间的粗糙物的数量n'(a)为

(4)接触力学模型

总TCR中包括最大接触面积a_L和无因次接触面积A_r ^*在内的接触参数，这些接触参数由接触粗糙物的接触状态决定。要计算接触参数，应建立接触力学模型。

(41)整个界面上的总触点负载

当变形ω足够小时，根据赫兹的弹性接触理论，弹性变形发生在接触粗糙体上。当a＜a_pc时，所有接触粗糙面都在塑性变形下，然后在内圈/轴颈界面上的接触载荷表示为：

当a_pc＜a＜a_ec时，内圈/轴颈接口上的总接触载荷表示为：

其中，

当a＞a_ec时，内圈/轴颈接口上的总接触载荷表示为：

若D＝1.5，则内圈/轴颈接口上的总接触载荷表示为：

(42)整个界面上的实际接触面积

当参数a_L,a_ec,和a_pc已知时，计算总接触面积A_r。a_ec表示临界弹性接触面积，a_pc表示临界塑性接触面积。

若a_L＜a_c2，则在粗糙处会发生完全塑性变形，总的实际接触面积表示为：

若a_c2＜a_L＜a_c1，粗糙处会发生弹塑性变形，总的实际接触面积表示为：

若a_L＞a_c1，粗糙处会发生弹性变形，总的实际接触面积表示为：

至此，确定了总热阻R_c中除最大接触面积a_L以外的参数。为了确定最大接触面积a_L，应考虑接触界面的力平衡。为了保持平衡，由接触粗糙物的变形产生的总接触载荷F应该等于外部载荷F_ext。即，总接触载荷F等于接触压力P与标称表面积A_a的乘积。

F＝F_ext＝P·A_a (48)

接触压力P是由轴承内圈与轴颈之间的过盈配合I产生的，表示为：

其中，D₁,D₂和D₃分别表示内圈的外径、轴颈的内径以及接触直径；E₁表示轴承内圈材料弹性模量；E₂表示轴径材料弹性模量；v₁表示轴承内圈材料泊松比；v₂表示轴径材料泊松比。

因此，可以通过等式(46)、(48)和(49)确定接触点的最大接触面积a_L，然后根据等式(39-b)计算总热阻R_c。

如图20(a)所示，获得了轴承接头的TCR与过盈配合之间的关系。随着过盈配合增大，TCR降低。轴承内圈与轴颈之间的过盈配合为1μm，则当考虑表面接触系数λ时，轴承接头的TCR约为4200m^2℃/W，而当表面接触系数时，轴承接头的TCR约为5000m^2℃/W。如图20(b)所示，未考虑表面接触系数λ的TCR大于考虑表面接触系数λ的TCR，这是因为表面接触系数λ小于1。接触应力不是无限大的，当接触半径小且趋于零时具有一定的值。接触面积随接触半径而增加，接触应力随接触半径而减小。而且，当接触半径逐渐增加到无穷大时，表面接触系数λ趋于为1。

(5)更新内圈/轴颈接口的TCR

轴颈与轴承内圈之间的过盈配合随主轴系统的热变形而变化，并且偏离初始过盈量I；动态过盈配合表示为：

I'＝I+ΔI (50)

其中，ΔI标致轴颈和轴承内圈的径向热变形之间的差值，ΔI可通过热-结构耦合特性仿真模型的每个迭代子步骤获取，从而求取得到更新的最大接触面积a_L，以更新轴承内圈与轴颈之间的接触热阻。

3.2、外圈/轴承座之间的TCR

一般而言，轴承外圈与轴承箱间隙配合。同样，表面接触系数如等式(45-a)所示，考虑圆柱表面对TCR的影响。间隙配合的TCR，也即轴承外圈与轴承座之间的接触热阻R为：

其中，δ_ring表示外圈厚度；λ_air表示空气的热导率；λ_ring表示轴承外圈的热导率；R_s为表示接触半径；δ_gap表示间隙，且：δ_gap＝δ_gap0-(T_ring-T_housing)·α·r_housing，δ_gap0表示初始间隙；T_ring表示外圈的平均温度；T_housing表示轴承座的平均温度；α表示热膨胀系数；r_housing表示接触区域的半径。径向热膨胀发生在轴承外圈上，然后减小了游隙。

3.3、滚珠/环之间的TCR

滚珠与轴承圈之间的TCR由赫兹接触面的形状和尺寸决定。赫兹接触表面与轴承的几何形状和轴承的内部接触应力分布有关。滚珠/环界面的TCR表示为：

其中，λ₁表示导热系数；u表示积分变量；a和b表示赫兹接触区域的长半轴和短半轴，且a和b表示为：

其中，k表示第一类文积分；Σ表示第二类椭圆积分；Q表示接触载荷；

通过使用第一种椭圆积分，方程(52)被改写为：

ψ表示与接触表面的尺寸有关的几何参数，为：

其中，

每个球和环之间的TCR为：

其中，λ_b和λ_ring表示球和环的热导率。

修改后的预载荷被代入轴承的非线性方程组，以修改滚珠和环之间的法向接触载荷Q。正常接触负载Q被引入到等式(53)中，修改赫兹接触区域的长半轴和短半轴。最后，每个球和环之间的TCR可以通过等式(56)进行计算。根据以上分析，计算出外圈/轴承座界面，球/圈槽界面和其他触点界面的TCR，如表4所示。

表4接触点的初始TCR

下面对本实施例倾斜工作条件下的主轴系统热误差建模方法的具体实施方式进行详细说明。

主轴系统对单个热源的响应为：

其中，C₁和C₂表示要确定的系数；T(x,t)表示轴温；d和T_∞表示轴的直径和环境温度；h和λ分别为对流系数和热导率；x和t表示轴向位置和运行时间；ρ和c表示密度和比热容；

表示积分变量。

在运行过程中，主轴系统中的轴被加热，对于带有固定支撑装置的主轴系统，从固定端到支撑/自由端的热膨胀不断增加。因此，获得的热膨胀为：

其中，α表示热膨胀系数；L表示轴长；C₃,C₄,C₅和C₆表示要确定的系数。

用等式(2)直接计算热误差很困难，原因如下：一方面，由于无法确定系数C₁、C₂、C₃、C₄、C₅和C₆，因此很难用等式(2)预测热误差。另一方面，广义的热误差模型具有重要意义，等式(2)为误差建模提供了有用的思路，等式(2)中所示的定积分是运行时间的指数函数。根据热误差的行为，热误差模型本质上是指数的，其中运行时间t、热误差平衡过程的时间常数τ、初始热误差δ₀与稳态下的热误差δ_t作为输入，从而得到主轴系统的热误差模型为：

δ(t)＝δ₀+(δ_t-δ₀)(1-e^-t/τ) (3)

其中，δ₀和δ_t分别为初始状态和稳态下的热误差；τ为热误差平衡过程的时间常数。

初始状态下的热误差δ₀为0，热误差建模的核心是确定热误差平衡过程的时间常数τ和稳态下的热误差δ_t，这是通过热-结构耦合特性仿真建模获得的。具体的，通过热-结构耦合特性仿真建模获得时间常数τ和稳态下的热误差δ_t的方法如下：

采用如上所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，在迭代子步骤中计算主轴系统的热误差，且当

时，主轴系统达到热平衡状态，此时的时间t_i为热平衡时间，则

此时的热误差E_i为主轴系统的热平衡误差δ_t；其中，ΔE_max＝E_max-E₀，E_max和E₀分别表示主轴系统的最大热误差和初始热误差；f_i表示t_i时刻主轴系统的热误差变化相对于主轴系统的最大热误差变化的百分比；ΔE_i表示t_i时刻的热误差变化，且ΔE_i＝E_i-E₀。

本实施例的主轴系统热误差建模方法的具体方式如下。

4、热误差建模

4.1、实验装置和测量方法

主轴系统严格按照安装规范进行安装，以确保滴定工作条件效果的可靠性。测量设备和功能如下：开发同步采集系统，并且该同步采集系统是基于NI-SCXI-1600构建的。它用于通过PT100精密磁性温度传感器收集轴承外圈的温度，并通过涡流位移传感器来测量主轴系统的热变形，如图21所示。使用PT100温度传感器T1和T2收集前后轴承的温度。涡流位移传感器S1，S2和S3用于分别测量Y，X和Z方向的热位移。测量系统的架构为NI-SCXI-1600，温度和位移变送器用于实现温度和位移信号的滤波、放大以及A/D转换。

4.2、稳态热误差和时间常数

(1)温度场建模

主轴系统的热-结构耦合特性在6000r/min，8000r/min，10000r/min，12000r/min和15000r/min的条件下进行仿真。图22示出了当转速为15000r/min并且间隔环的平行度误差为30μm时，轴倾斜工作条件下的的温度场。温度场分布不均匀，转子温度高于其他部件。整个主轴系统的最高温度约为76.375℃，并且由于TCR的作用，轴承内圈与轴颈之间以及轴承外圈与轴承座之间的温度差异非常明显。因此，对于热-结构耦合特性的建模和分析，TCR对温度场的影响已不容忽视。表面接触系数的引入对于预测轴承接头的TCR是有效的。循环冷却系统可以消散定子产生的大部分热量，因此定子的温度不高。转子与轴接触，然后由转子产生的大量热量传导到轴中，然后热量沿着轴传导。内圈的温度高于轴承外圈的温度，因为内圈的散热比外圈的散热差。

图23显示了轴倾斜工作条件下的前轴承各部件的温度。外圈和内圈的最高温度分别为41.371℃和48.921℃。外圈和内圈的最高温度之差约为7.5℃。球的最高温度为47.268℃。外圈的最高温度发生在外圈与球的接触区域。内圈的最高温度发生在内圈与间隔圈接触的末端。内圈在轴向上的温度差为7.9℃。内圈的最高温度发生在球/内圈的接触区域。

图24所示，获得轴倾斜工作条件下的轴承部件的温度。轴承部件的温度随着转速的升高而逐渐升高。在15000r/min时，最高轴承温度达到64.4℃，由于内圈和球的散热条件相对较差，内圈和球的温度差为

内圈的热量通过与轴的热传导而消散。内圈和外圈之间的温度差不会随转速显着变化，大约为7℃。轴承外圈的最高温度为54.4℃。

通过热特性分析模型预测前后轴承的温度。此外，将轴倾斜工作条件下的轴承外圈的测量数据与热特性分析模型的预测结果进行比较，如图25所示。由于轴承的热量增加，因此预测温度和测量温度均随着间隔圈的平行度误差而增加。带有间隔环的平行度误差，如图6所示。当达到热平衡状态时，轴倾斜工作条件下的测得的温度围绕热平衡温度波动，并且该波动是由循环冷却系统的工作原理引起的。循环冷却系统在出口温度高于临界温度时开始工作，然后在出口温度低于临界温度时停止工作。测得的温度和预测的温度都显示出趋势，即温度从一开始就缓慢升高，并在早期急剧升高，然后在热量产生与向大气散热的平衡过程中逐渐饱和到最终温度。更重要的是，在不同的平行度误差下，预测结果与测得的温度非常匹配。因此，验证了轴承准静态力学模型的有效性。结论是，平行度误差越大，轴承产生的热量就越大，轴承套圈的温升也就越高。以上分析初步验证了热-结构耦合特性建模方法的有效性。

(2)稳态热误差和时间常数

将温度作为主轴系统变形分析的输入，并通过热-结构耦合特性仿真模型获得倾斜工作条件下的高速主轴系统的热变形，如图26所示。，轴倾斜工作条件下主轴系统不均匀，主轴系统的最大热变形为71.825μm，发生在轴的最右端。后轴承受到径向和轴向支撑，并且轴在后轴承方向上的膨胀受到限制。前轴承受到径向支撑，轴在前轴承方向上的膨胀不受限制。因此，最大的热变形发生在轴的最右端。

轴倾斜工作条件下轴承的热应力分布是通过热-结构耦合特性分析获得的，如图27所示。在滚珠/环界面接触区域产生的应力通常很高。可以看出，滚珠的热应力以15000r/min的速度达到约367.7MPa。滚珠与凹槽之间有两个高应力区域，凹槽的最大应力约为647.3MPa，它发生在凹槽的边缘。

主轴系统倾斜工作条件下以特定的转速运行。采集系统在热平衡过程中每隔一秒钟收集一次前后轴承的温度。进行了四组实验。然后，测量主轴系统在X，Y和Z方向上的热变形，并将测得的数据与预测结果进行比较，如图28所示。在整个热平衡过程中测量数据，验证了热-结构耦合特性建模方法。模拟的热误差与实测数据吻合，验证了热-结构耦合特性建模方法。X方向，Y方向和Z方向的热误差很明显，并且它们会随着时间的推移而增加，并在达到热平衡状态时最终变得稳定。在整个热平衡过程中，在X，Y和Z方向上测得的热误差与预测结果非常吻合。此外，平行度误差越大，X，Y和Z方向的热误差就越大。如图10(b)所示，倾斜工作条件下轴承的轴向刚度随倾角而减小，这意味着当倾角增大时，轴向热膨胀更容易发生。如图10(a)所示，倾斜工作条件下在平行于倾斜方向的方向上的径向刚度随着倾斜角度而增加，而在垂直于倾斜方向的方向上的径向刚度随着倾斜角度而减小。而且，如图6所示，轴承的发热随着倾斜角度而增加，这导致轴的更大的热膨胀。因此，轴的轴向热伸长率随着倾斜角度而增加。对于X和Z方向的热变形，由于径向刚度的变化，也可以获得类似的趋势。

通过使用温度平衡过程的时间常数τ的定义作为参考来定义轴倾斜工作条件下热误差平衡过程的时间常数τ。即，主轴系统已经达到热平衡，并且当f_i＝95％时热平衡时间为t_i，如图29所示。然后获得以下关系。

然后，获得轴倾斜工作条件下的稳态热误差δ_t，如图30所示。可以看出，轴倾斜角度越大，稳态热误差就越大，热平衡越短时间。稳态下的热误差δ_t随倾斜角ξ线性变化。因此，使用直线拟合确定热误差与倾斜角ξ之间的关系，如等式(58)所示。热误差作为倾斜角ξ的线性函数的拟合极大地促进了热误差补偿。本实施例进行了热-结构耦合特性仿真建模，以模拟在不同转速轴倾斜工作条件下的热误差，然后预先拟合了不同方向上的稳态热误差方程。在热误差补偿中，调用稳态热误差方程来计算稳态下的热误差δ_t。

x方向的热误差拟合为：

δ_x＝-0.352ξ-12.44 (58-a)

y方向的热误差拟合为：

δ_y＝-0.4675ξ-3.334 (58-b)

轴向的热误差拟合为：

δ_z＝-0.3375ξ-15.220 (58-c)

通过轴倾斜工作条件下热-结构耦合特性模型确定热误差平衡过程的时间常数。根据以上分析，当热误差变化达到最大热误差的95％时，判断主轴单元已达到热平衡。根据图28中所示的建模结果，获得每个速度下的热误差平衡过程的时间常数。热平衡时间是时间常数_τ的三倍。然后获得用于不同并行度误差的时间常数_τ，计算出不同倾角下的时间常数。可以看出，倾斜角越大，热误差越大，并且热误差平衡过程的时间常数越短。然后表5列出了在不同倾斜角下的热平衡过程的时间常数。各速度下的热平衡的时间常数相差不大，主轴系统的热误差平衡过程的时间常数为各速度下的热平衡时间常数的平均值。即，主轴系统在X方向，Y方向和Z方向上的热误差平衡的时间常数分别约为τ_x＝283.2s,τ_y＝210.6s,andτ_z＝270.3s。

表5热误差平衡过程的热平衡时间

根据以上分析，建立了倾斜工作条件下主轴系统的热误差模型。根据不同方向上的稳态热误差和热平衡的时间常数，可以得出热误差方程：

δ(x)＝(-0.352ξ-12.44)·(1-e^-t/282.2) (59-a)

δ(y)＝(-0.4675ξ-3.334)·(1-e^-t/260.1) (59-b)

δ(z)＝(-0.3375ξ-15.220)·(1-e^-t/270.3) (59-c)

5、热误差补偿与验证

为了充分验证倾斜工作条件下热误差模型的有效性，将转速分布设置为图31。转速逐步增加，然后下降并趋于稳定。将等式(59)中建立的热误差补偿模型嵌入补偿器中，然后进行误差补偿。

本实施例的热误差模型和MLRA模型用于补偿热误差。图32显示了三种热误差模型的补偿结果。MLRA模型的残差波动远大于本实施例的热误差模型，本实施例的热误差模型减少了三个方向的热误差，其补偿效果优于MLRA模型，验证了热误差建模和补偿方法。

如表6所示，评估了倾斜工作条件下以上两种补偿模型的补偿效果。对于本实施例的热误差模型，X方向上的绝对平均残基和最大残基分别约为0.506μm和1.871μm。对于MLRA模型，X方向的平均残留量和最大残留量分别约为2.395μm和8.678μm。本实施例的热误差模型的补偿效果远优于MLRA模型。然后，验证了当前的热误差建模方法，该方法足够准确，可以在后续的加工过程中进行热误差补偿。本实施例的热误差模型的出色补偿效果归因于误差产生机理的揭示和精确的热-结构耦合特性仿真建模。误差产生机制使其能够将热误差表示为指数函数，其中稳态时的热误差为δ_t，热误差平衡过程的时间常数为_τ，运行时间为其自变量。通过热-结构耦合特性仿真建模，可以准确确定热误差平衡过程的时间常数_τ和稳态下的热误差δ_t。

表6补偿效果评估。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (9)

1.一种倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：包括如下步骤：

1)创建倾斜状态下的主轴系统的热-结构耦合特性仿真模型，并将倾斜工作条件下的主轴系统的热源、热边界条件和力学边界条件应用于所述热-结构耦合特性仿真模型；

倾斜工作条件下的主轴系统的热源包括：轴承的发热、电动机的发热；

倾斜工作条件下的主轴系统的热边界条件包括：强制对流换热系数、自由对流换热系数、轴承内圈与轴颈之间的接触热阻、轴承外圈与轴承座之间的接触热阻以及轴承滚珠与环之间的接触热阻；

倾斜工作条件下的主轴系统的力学边界条件包括：倾斜状态下的主轴轴承的轴向接触刚度和径向接触刚度；

2)根据所述倾斜状态下的热-结构耦合特性仿真模型对主轴系统进行瞬态分析；

3)比较相邻两个迭代子步骤中的主轴的温度，若满足收敛条件，则进入步骤4)；若不满足收敛条件，则进入步骤5)；所述收敛条件为：

T_i-T_i-1＜ΔT

其中，T_i表示第i次迭代子步骤中得到的主轴温度；ΔT表示收敛温差；

4)结束分析，保存倾斜状态下的轴系统中所有温度节点的温度；

5)更新倾斜状态下的主轴系统热源和热边界条件；

6)循环步骤2)和步骤3)，直至满足收敛条件；

轴承发热量的更新方法为：

更新润滑剂粘度v：

v＝32×e^{-0.0242×(T-40)}

其中，T为从倾斜状态系统热-结构耦合特性仿真模型仿真中提取的温度；

更新主轴的轴向变形，得到更新后的主轴轴向变形δ_a为：

δ_a＝δ_preload±Δδ_y

其中，δ_a为更新后的主轴的轴向变形；δ_preload表示初始预紧力作用下的主轴的轴向变形；Δδ_y表示通过热-结构耦合特性仿真模型获得倾斜状态主轴的轴向热伸长；

根据轴向变形与预紧力之间的相关性，采用牛顿-拉夫森迭代算法计算实际接触角α：

δ_a＝BD_bsin(α-α_o)/cosα

其中，α_o为初始接触角；B为曲率半径和；D_b表示滚珠的直径；

通过将实际接触角α获取实际预紧力P_a：

其中，K_n＝K·D_b ^0.5/B^1.5表示轴向变形常数，K表示接触刚度；Z表示滚珠的数量；

接触载荷与接触变形的关系为：

对该等式两端求导，即可将接触面积的接触刚度表示为

其中，Γ和κ分别表示第一类和第二类的完整椭圆积分；E'表示等效弹性模量；κ表示椭圆偏心率参数；∑ρ表示在正交平面上的接触点处的接触体的曲率之和；ε表示接触变形；Q表示接触载荷；

实现对接触刚度的更新。

2.根据权利要求1所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：轴承的发热Q_f为：

Q_f＝Q_l+Q_v+Q_s+Q_gy

其中，Q_l、Q_v、Q_s和Q_gy分别为外部负载、润滑剂的粘滞摩擦、球的旋转运动和陀螺运动产生的热量；

Q_l＝M_l·ω

Fβ＝max(0.9F_a/tanα-0.1F_r,F_r)

其中，M_l表示由外部载荷引起的摩擦扭矩；ω表示轴承的转动角速度；f₁表示与轴承结构和载荷油管的系数，F_s表示等效静载荷；C_s表示基本额定静载荷；d_m表示轴承的节圆直径；F_β表示等效载荷；F_a表示轴承的轴向载荷；F_r表示轴承的径向载荷；α表示实际接触角；

Q_v＝M_v·ω

其中，M_v表示由润滑剂的粘滞摩擦引起的摩擦扭矩；n表示内圈的转速；f_o表示与轴承类型和润滑条件有关的常数；v_o表示在工作温度下润滑剂的运动粘度；

Q_s＝Z·M_si·ω_si+Z·M_so·ω_so

其中，M_s表示由滚珠的旋转运动引起的旋转摩擦扭矩；η_oil表示滚珠与槽之间的摩擦系数；Q_i/o表示施加在内沟道或外沟道上的负载；a_i/o表示接触椭圆的半长轴；ε_i/o表示第二类的完整椭圆积分；M_si表示旋转运动引起的滚珠与内沟道之间的旋转摩擦扭矩；ω_si表示滚珠相对于内沟道的转动角速度；M_so表示旋转运动引起的滚珠与外沟道之间的旋转摩擦扭矩；ω_so滚珠相对于外沟道的转动角速度；Z表示滚珠的数量；

Q_gy＝M_gy·ω

其中，M_gy表示由陀螺力矩引起的摩擦转矩；β表示滚珠的运动姿势；ω_c表示方位角

处的公转角速度；ω_b表示每个滚珠的转动角速度；ρ表示滚动体密度；D_b表示滚珠的直径。

3.根据权利要求2所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：轴承的轴向接触刚度K_a为：

其中，K_aj表示第j^th个滚珠的轴向接触刚度；K_aij表示第j^th个滚珠/内沟道界面之间的轴向接触刚度；K_aoj为第j^th个滚珠/外沟道界面之间的轴向接触刚度；Z表示滚珠的数量；

倾斜状态下轴承的径向接触刚度K_r为：

其中，K_rj表示第j^th个滚珠的径向接触刚度；K_rij表示第j^th个滚珠/内沟道界面之间的径向接触刚度；K_roj为第j^th个滚珠/外沟道界面之间的径向接触刚度；

为第j^th个滚珠的方位角；

且：

其中，Γ表示第一类和第二类的完整椭圆积分；E'表示等效弹性模量；κ表示椭圆偏心率参数；∑ρ表示在正交平面上的接触点处的接触体的曲率之和；ε表示接触变形；αij表示内接触角；α_oj表示外接触角；Q_ij和Q_oj表示内沟道和外沟道分别施加到j^th滚珠的法向接触载荷；且：

其中，λ_ij和λ_oj为系数；D_b表示滚珠直径；M_gj表示陀螺力矩；F_cj表示j^th滚珠的离心力；且：

其中，m表示滚珠质量；d_m表示轴承的节圆直径；ω_c和ω_b表示滚珠的公转和自转角速度；ω表示轴承内圈的转动角速度；

α_o表示初始接触角；ω_cj表示第j个滚动体公转角速度；

且：

其中，X_1j表示内滚道沟曲率中心到球中心最终位置的轴向距离；X_2j表示内滚道沟曲率中心到球中心最终位置的径向距离；A_1j和A_2j分别表示内沟道、外沟道曲率中心之间的轴向和径向距离；f_i表示内沟道曲率系数；δ_ij表示接触载荷作用下的正常接触变形；f_o表示外沟道曲率系数；δ_oj表示接触载荷作用下的正常接触变形；其中：

(A_1j-X_1j)²+(A_2j-X_2j)²-[(f_o-0.5)D_b+δ_oj]²＝0

其中，A_o表示内外沟道之间的初始距离；R_o表示外沟道曲率中心的走线半径；

表示方位角；α°表示静载荷作用下内圈与外圈的接触角；δ_a和δ_r分别为内圈和外圈之间的相对轴向位移和径向位移；θ表示轴承的倾斜角，且

其中ξ和d分别表示间隔环的平行度误差和外径。

4.根据权利要求1所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：电动机的发热Q为：

Q＝P_in(1-η)

其中，P_in为电机输入功率；η为电机效率。

5.根据权利要求1所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：自然对流换热系数和强制对流换热系数表示为：

其中，h为对流换热系数；Nu为努塞尔数；L为水力直径；λ为流体的导热系数。

6.根据权利要求1所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：轴承内圈与轴颈之间的接触热阻表示为：

其中，k₁和k₂分别为两种接触固体材料的导热系数；k₃为间隙中介质的导热系数；

为无因次接触面积；A_n表示名义接触面积；H_s为接触热导率；且：

其中，h_cp表示塑性变形引起的接触热导率；h_ep表示弹性变形引起的接触热导率；h_cep表示弹塑性接触热导率；n(a)表示面积介于a和a+da之间的粗糙体的数量；a为接触面积；a_c1为接触点开始发生塑性变形时的临界接触面积；a_c2为完全塑性变形时的临界塑性变形；

其中，R_s为表示接触半径；D为分形维数；f(a)为模板函数；a_L为最大接触面积；

所述步骤5)中，轴承内圈与轴颈之间的接触热阻的更新方法为：

为了保持平衡，由接触粗糙物的变形产生的总接触载荷F应该等于外部载荷F_ext，即：

F＝F_ext＝P·A_a

其中，H表示较软材料的硬度；A_a为标称表面积；P为接触压力；且接触压力P是由轴承内圈与轴颈之间的过盈配合I产生，并表示为：

其中，D₁,D₂和D₃分别表示内圈的外径、轴颈的内径以及接触直径；E₁表示轴承内圈材料弹性模量；E₂表示轴径材料弹性模量；v₁表示轴承内圈材料泊松比；v₂表示轴径材料泊松比；

轴颈与轴承内圈之间的过盈配合随主轴系统的热变形而变化，并且偏离初始过盈量I；动态过盈配合表示为：

I'＝I+ΔI

其中，ΔI表示轴颈和轴承内圈的径向热变形之间的差值，ΔI可通过热-结构耦合特性仿真模型的每个迭代子步骤获取，从而求取得到更新的最大接触面积a_L，以更新轴承内圈与轴颈之间的接触热阻。

7.根据权利要求1所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：轴承外圈与轴承座之间的接触热阻R为：

其中，δ_ring表示外圈厚度；λ_air表示空气的热导率；λ_ring表示轴承外圈的热导率；R_s表示接触半径；δ_gap表示间隙，且：

δ_gap＝δ_gap0-(T_ring-T_housing)·α·r_housing

其中，δ_gap0表示初始间隙；T_ring表示外圈的平均温度；T_housing表示轴承座的平均温度；α表示热膨胀系数；r_housing表示接触区域的半径。

8.根据权利要求1所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，其特征在于：轴承滚珠与环之间的接触热阻为：

其中，λ_b表示滚珠的热导率；λ_ring表示环的热导率；ψ表示与接触表面的尺寸有关的几何参数；a为赫兹接触区域的长半轴。

9.一种倾斜工作条件下的主轴系统热误差建模方法，其特征在于：主轴系统的热误差模型为：

δ(t)＝δ₀+(δ_t-δ₀)(1-e^-t/τ)

其中，δ₀表示主轴系统的初始状态热误差，δ_t表示主轴系统在热平衡状态下的热误差；τ表示热误差平衡过程的时间常数；

采用如权利要求1-8任一项所述倾斜工作条件下的主轴系统热-结构耦合特性模型创建方法，在迭代子步骤中计算主轴系统的热误差，且当

时，主轴系统达到热平衡状态，此时的时间t_i为热平衡时间，则

此时的热误差E_i为主轴系统的热平衡误差δ_t；其中，ΔE_max＝E_max-E₀，E_max和E₀分别表示主轴系统的最大热误差和初始热误差；f_i表示t_i时刻主轴系统的热误差变化相对于主轴系统的最大热误差变化的百分比；ΔE_i表示t_i时刻的热误差变化，且ΔE_i＝E_i-E₀。

Images