考虑力热耦合效应的高性能电主轴热动态特性确定方法与流程

本发明涉及机械领域，尤其涉及一种考虑力热耦合效应的高性能电主轴热动态特性确定方法。

背景技术：

随着时代的发展，高性能电主轴已成为现代高速加工机床的核心部件，为机床加工精度提供了极大的保障。而电主轴在运转过程中不可避免的产生大量的热量，使电主轴系统产生不均匀热变形，同时还会导致系统支承刚度软化，从而影响零件的加工质量，严重时还会导致系统失效而无法工作，这严重阻碍了高速加工技术的持续发展。电主轴热动态特性主要包括稳态温度场的分布和系统支承刚度的大小，有效的预测电主轴稳态温度场和准确的计算系统支承刚度的大小对高速加工技术的发展尤为重要。

目前，计算电主轴稳态温度场的方法主要有数值计算法和有限元仿真计算法，数值计算法是通过建立电主轴的热网络模型并求出特定条件下的电主轴热源发热量，然后经过解方程求得电主轴稳态温度场，这种方法求解简单，计算效率高，但需要考虑热源发热量的时变性。有限元仿真计算法是通过建立系统的有限元模型，然后将求得的热源发热量、换热条件等作为边界条件加载到有限元模型上，这种方法可准确获取电主轴每个部分的温度大小，但求解效率往往不高，且忽略了边界条件的时变性。而计算电主轴支承刚度的方法则主要是解析法，对机械式电主轴来说，目前已有的考虑热效应的轴承刚度计算方法则是在轴承运动学方程中加入轴承热位移，目前计算主轴-轴承系统的支承刚度则主要是采用此方法，而对结构复杂、轴承数量多的高性能电主轴来说，这种方法还有待改进。

而且在现有技术中，电主轴热态特性和动态特性往往作为两个独立的研究目标，但实际上，电主轴在高速运转过程中，将不可避免的出现力热耦合效应，即电主轴运转时产生的热会改变电主轴内部力的大小，而力的分布情况又会改变热源发热量的大小，只有综合考虑才能更为准确的得到电主轴的热动态特性。

技术实现要素：

本发明的目的在于提供一种考虑力热耦合效应的高性能电主轴热动态特性确定方法，本发明充分考虑了电主轴高速运转过程中的力热耦合效应，通过数值计算方法得到电主轴稳态温度场分布和系统支承刚度的大小，更为准确的得到电主轴的热动态特性。

为解决上述问题，本发明的技术方案是：

一种考虑力热耦合效应的高性能电主轴热动态特性确定方法，其特征在于，包括如下步骤：

步骤一、建立轴承力热耦合模型，分析考虑力热耦合效应的轴承内部几何相容关系并进行轴承受力分析和轴承产热分析；

步骤二、建立基于热阻的高性能电主轴热传递网络模型，并分析不同类型热阻的计算表达式；

步骤三、建立电机损耗模型，分析电机产热；

步骤四、分析不同安装方式下的轴承套圈热位移变化情况；

步骤五、拟定考虑力热耦合效应的高性能电主轴热动态特性的分析计算流程，并选定算法进行求解；

步骤六、输出计算结果，得到电主轴热动态参数。

进一步的改进，所述电主轴热动态参数包括选定节点温度和系统各轴承支承刚度随主轴转速的变化规律。

进一步的改进，所述步骤一包括如下步骤：建立轴承力热耦合模型，进行轴承几何分析、受力分析和产热分析：

设轴承外圈曲率中心固定不动，轴承内圈随转轴一起旋转，当轴承受到载荷f＝{fx,fy,fz,mx,my}时，轴承内外圈产生相对位移δ＝{δx,δy,δz,θx,θy}，其中，，fx、fy、fz、mx、my分别为轴承受到的x、y、z方向的外载荷和x、y方向的外力矩，j为轴承球体序数，z为轴承球体个数；考虑到轴承产热和电机产热，电主轴系统将会产生热变形，从而轴承内外圈还会产生热位移u＝{ua,ur}，ua，ur分别表示轴承套圈轴向相对热位移和径向相对热位移，此外还有轴承内圈的离心位移δc，得到角位置ψj处的轴承内外圈轴向和径向相对距离：

a1j＝bdbsinα0+δa

a2j＝bdbcosα0+δr+δc

式中，a1j为轴承内外圈轴向相对距离，a2j为轴承内外圈径向相对距离，b＝fi+fo-1，fi、fo为内、外滚道曲率半径系数，db为球体直径，α0为轴承初始接触角，δa为轴承内外圈轴向总相对位移，δr为轴承内外圈径向总相对位移，δc为轴承内圈离心位移；

建立o-xyz主轴坐标系，o1-x1y1z1为轴承一坐标系，o2-x2y2z2为轴承二坐标系；o点建在转轴端面中心，o1表示轴承一的几何中心，o2表示轴承二的几何中心；

δa＝δz+riθxcosψj-riθysinψj+ua

对轴承二：

δa＝-δz-riθxcosψj+riθysinψj+ua

式中，δz为轴承内外圈轴向相对位移，ri为内圈曲率中心圆半径，θx为轴承内外圈x向相对角位移，θy为轴承内外圈y向相对角位移，ψj为轴承滚动体位置角，ua为轴承内外圈轴向相对热位移；

轴承内外圈径向总相对位移：

式中，δx、δy分别为轴承内外圈x向和y向相对位移，θr为最大变形处的滚动体与位置角为0的滚动体之间的夹角，且ur为轴承内外圈径向相对热位移。

轴承内圈的离心位移δc：

式中，ρi表示轴承内圈密度，ωi表示轴承内圈转速，ei表示轴承内圈弹性模量，μi表示轴承内圈泊松比，d表示轴承内圈直径，di表示内圈沟道直径；

得到角接触球轴承沟道接触的变形几何相容方程为：

(a1j-x1j)²+(a2j-x2j)²-[(fi-0.5)db+δij]²＝0

x1j²+x2j²-[(fo-0.5)db+δoj]²＝0

式中，x1j、x2j分别为球体中心与外沟道曲率中心的轴向距离和径向距离，δij、δoj分别为内、外圈与球体接触的弹性变形趋近量，db为球体直径；

在通过球体中心和轴承轴线的平面上，角位置ψj处的球体所受载荷，设λij和λoj分别表示内、外圈沟道控制系数，对于高速情况，设球体所受陀螺力矩完全被球体与外圈沟道的接触摩擦力抵消，这时，取λij＝0,λoj＝2；根据hertz点接触变形关系及球体受力平衡关系可得：

式中，kij、koj分别为内、外圈滚道与球体接触的载荷变形系数，αij、αoj分别为轴承内、外圈与球体接触角，mgj、fcj分别为第j个球体受到的陀螺力矩和离心力；

为了求出稳定状态时轴承的各向位移值，建立轴承所受载荷的平衡方程，假设轴承受到联合载荷f＝{fx,fy,fz,mx,my}的作用在主轴坐标系下，轴承载荷与位移的关系可用下式表达：

定压预紧时：

对轴承一：

对轴承二：

式中，fx、fy、fz、mx、my分别为轴承受到的x、y、z方向的外载荷和x、y方向的外力矩，j为轴承球体序数，z为轴承球体个数；

定位预紧时，轴承载荷平衡方程只需去掉定压预紧载荷平衡方程中的z向方程即可。

电主轴工作过程中，由于轴承滚动体与内外圈的摩擦作用，轴承将产生热量，下面给出轴承产热量的计算方法：

对于位置角ψj处的球体有：

式中，mij、moj分别为球体与内、外滚道的摩擦力矩，f0、f1是与轴承类型和润滑、负载有关的系数，η0为润滑油的运动粘度，ωcj为第j个球体的公转角速度，qij、qoj分别为球体与内、外滚道的接触应力，qimax、qomax分别为球体与内、外滚道接触的最大应力；

对于外滚道控制，仅需考虑球体与内滚道的自旋摩擦力矩，即：

式中，msij为球体与内滚道的自旋摩擦力矩，μsi为球体与内滚道的摩擦系数，aij为球体与内滚道的赫兹接触椭圆长半轴，σij为球体与内滚道赫兹接触椭圆的第二类完全积分；

故内、外滚道接触区的摩擦发热量分别为：

hij＝ωcj.mij+ωbj.msij

hoj＝ωcj.moj

式中，hij、hoj分别为球体与内、外滚道接触区的摩擦发热量，ωbj为球体的自旋角速度。

进一步的改进，所述步骤二：建立电主轴热传递网络模型，分析不同类型热阻的计算表达式步骤如下：

电主轴在运转过程中，存在着两种主要热传递方式：热传导、热对流；

a.对于热传导热阻：

热传导热阻分为轴向、径向热阻；根据fourier定理，圆筒和圆柱的轴向热传导热阻r1为：

式中，l为圆筒或圆柱轴向长度，k为圆筒或圆柱热导率，sa为圆筒或圆柱轴向截面积；

圆筒径向热传导热阻r2为：

式中，d2、d1分别为圆筒外径和内径；

圆柱径向热传导热阻r3为：

球体热传导热阻r4为：

式中，kb为球体热导率；

油脂润滑时内、外滚道的润滑脂热传导阻抗r5和r6分别为：

式中，kl为润滑脂热导率，do为轴承外圈沟道直径，bi、bo分别为轴承内、外圈宽度；

b.对流换热热阻

根据fourier定理，圆筒和圆柱的轴向、径向对流换热热阻r7为：

式中，h为对流换热系数，s为圆筒或圆柱与外界对流换热面积；

换热系数的计算方法如下：

轴承座与外界环境发生自然对流换热和辐射换热，其复合换热系数h1为：

式中，nu1为轴承座表面与空气对流换热的努赛尔数，λe为环境介质热导率，de为换热表面特征直径。

其中，g为重力加速度，βe为环境介质热膨胀系数，th为轴承座表面温度，te为周围环境温度，pr为环境介质普朗特数，ηe为环境介质的运动粘度；

冷却套中的流体与周围壁面之间的换热系数h2为：

式中，nu2为流体与壁面对流换热的努赛尔数，λf为流体热导率，df为通道特征尺寸；

其中，nu2＝0.0225re1^0.8pr^0.3，re1为雷诺数，且uf为流体的运动速度，ηf为流体的运动粘度；

流体运动速度qf为流体体积流量，af为通道截面积；

对于冷却水套，通道特征尺寸cf为通道截面周长；对于定转子气隙，df为定转子气隙长度；

电主轴系统旋转部分与周围环境流体对流换热系数h3为：

式中，deq为旋转部分当量直径，nu3为流体与旋转壁面对流换热的努赛尔数，λe为环境介质热导率

其中，re2为雷诺数，且ω为旋转部分转速；

结合面的接触热阻r8表示为：

式中，hc为结合面接触热导系数；a为结合面名义接触面积；

考虑微凸峰、界面流体介质的热传导，忽略气隙间的辐射换热，接触热导系数hc表示为：

式中，lg为结合面空隙厚度，k1、k2分别为两个配合零件的热导率，kf为界面间隙介质的热导率，a^*为结合面无量纲实际接触面积，且ac为实际接触面积。

建立电主轴热网络模型：

以轴承、电机的径向几何中心平面为分割线，沿径向布置一系列节点，其中结合面两端各布置一个节点，同时保证各个径向节点沿轴向方向对齐，建立系统热平衡方程。

进一步的改进，所述步骤三中：建立电机损耗模型，分析电机产热步骤如下：

工作过程中，定子损耗包括定子绕组上的铜耗和铁耗，忽略转子铁耗，考虑转子涡流损耗；

定子绕组上的铜耗pcu：

pcu＝mi²r

式中，m为定子绕组相数，i为定子绕组相电流，r为定子绕组阻值；

定子铁心损耗pfe由磁滞损耗ph、涡流损耗pe以及附加损耗pexc三者之和确定，即：

pfe＝ph+pe+pexc

其中，

式中，ms为定子质量，kh为磁滞损耗系数，ke为涡流损耗系数，kexc为附加损耗系数，f为电机工作频率，bm为最大磁通密度，α为无量纲磁滞损耗计算参数；

转子涡流损耗pr：

pr＝ρl∫∫j²ds

式中，ρ为永磁体的质量密度，l为永磁体导路长度，j为永磁体的感应涡流电流密度，s为永磁体导路横截面积。

进一步的改进，所述步骤四中分析不同安装方式下的轴承套圈热位移变化情况如下：

根据传热与热变形基础理论，等截面杆一维变温位移微分方程式为：

式中，u为轴向位移，x为轴向长度，t为温度，β表示表示杆件热膨胀系数；

当杆件无内热源稳定导热时，可得：

式中，c1、c2、c3、c4为积分常数；

不同安装方式下的轴承套圈热位移计算方法：

对于正向安装，轴承套圈轴向相对热位移ua：

ua＝δs-δh-ub(sinαi+sinαo)

对于反向安装，轴承套圈轴向相对热位移ua：

ua＝δh-δs-ub(sinαi+sinαo)

式中，δs为转轴与轴承配合处的轴向热位移，δh为轴承座与轴承配合处的轴向热位移，ub为轴承球体热膨胀位移，且ub＝βbδtbdb，βb为球体热导率，δtb为球体温度变化量，αi和αo分别为球体与轴承内外圈的接触角；

正反装轴承的轴承套圈径向相对热位移ur的计算方法为：

式中，βi为轴承内圈热导率，δti为轴承内圈温度变化量，βs为转轴热导率，δts为转轴与轴承配合处的温度变化量，μs为转轴泊松比，βh为轴承座热导率，δth为轴承座与轴承配合处的温度变化量，μh为轴承泊松比，do为外圈沟道直径。

进一步的改进，所述步骤五，拟定分析计算流程，选定算法进行求解包括如下步骤：

其中，采用newton-raphson迭代法求解轴承拟静力学模型和轴承力热耦合模型，采用基尔霍夫电流定律求解电主轴热网络模型，即在任一瞬时，流向某一节点的热流量之和恒等于由该节点流出的热流量之和，对当前热网络模型来说，设系统共有n个温度节点，且第i个节点有mi个热交换途径，应用基尔霍夫电流定律，则有：

式中，i,j＝1,2,…,n，ti表示第i个节点的温度，rij表示第i个节点与第j个节点之间的热阻，qi表示第i个节点自身产生的热流量，视其为流入热流量，并取流入热流量方向为正，流出热流量方向为负；mi表示热交换途径；将所有节点的热流平衡方程组装起来便构成了电主轴系统热流平衡方程，即：

{q}+[g]{t}＝0

上式即为电主轴热网络数学模型表示形式，其中{q}为系统生热矩阵，由轴承产热和电机产热计算求得，[g]为系统热阻矩阵的倒数，即热导矩阵，确定了{q}和[g]，便可得到系统温度分布矩阵{t}。

输入已知参数后，首先通过求解轴承拟静力学模型，得到轴承接触区的摩擦生热量，再根据电机损耗模型，计算得到电机产热量；将计算得到的各部分热阻值和轴承产热值、电机产热值代入电主轴热网络数学模型中，获得电主轴温度分布，由此可计算得到轴承与转轴和轴承座配合处的热位移，将得到的热位移值代入轴承力热耦合模型中并求解，若是收敛，则输出计算结果，否则便根据各个轴承的热位移情况修正各个轴承的内部几何关系，之后再对轴承拟静力学模型进行求解，如此反复，直到各模型达到设定的收敛精度后，方可输出最终结果。

进一步的改进，所述已知参数包括电主轴系统各零件材料、结构数据、轴承预紧载荷、主轴转速。

本发明的优点和有效权益在于：

(1)明确了考虑力热耦合效应的角接触球轴承的内部几何相容关系，建立了不同安装方式下的轴承载荷平衡方程；(2)考虑了电主轴内部各热源附近的接触热阻，建立了基于热阻的高性能电主轴热传递网络模型，模型概念清晰易懂；(3)建立了电机损耗模型，分析了电机产热，相较于经验公式计算法更加准确；(4)分析了考虑系统热效应的轴承套圈热位移变化情况，给出了相关计算公式；(5)充分考虑了电主轴运转过程中的力热耦合效应，拟定了准确清晰的计算流程，将这种建模方法和耦合求解思想引入电主轴这种复杂装配体的热动态性能预测领域中。

附图说明：

图1为本发明中建立的角接触球轴承拟静力学模型示意图

图2为本发明中建立的角接触球轴承力热耦合模型示意图

图3为本发明中建立的角接触球轴承安装形式图

图4为本发明中建立的位置角ψj处的球载荷示意图

图5为本发明中建立的简化后的电主轴二维结构图

图6为本发明中建立的电主轴热网络节点图

图7为本发明中建立的电主轴热网络模型图

图8为本发明中建立的电主轴热动态特性的分析计算流程图

图9a为轴承球体温度

图9b为电机温度

图10a轴承轴向支承刚度

图10b轴承径向支承刚度

图10c为轴承角向支承刚度。

具体实施方式：

下面结合附图对本发明的一种实施例做进一步的详细阐述，包括以下步骤：

步骤一：建立轴承力热耦合模型，进行轴承几何分析、受力分析和产热分析

轴承拟静力学模型如图1所示，在其基础上建立的轴承力热耦合模型如图2所示，假设轴承外圈曲率中心固定不动，轴承内圈随转轴一起旋转，当轴承受到载荷f＝{fx,fy,fz,mx,my}时，轴承内外圈产生相对位移δ＝{δx,δy,δz,θx,θy}，考虑到轴承产热和电机产热，电主轴系统将会产生热变形，从而轴承内外圈还会产生一定的热位移u＝{ua,ur}，此外还有轴承内圈的离心位移δc。由图2得到角位置ψj处的轴承内外圈轴向和径向相对距离：

a1j＝bdbsinα0+δa

a2j＝bdbcosα0+δr+δc

式中，a1j为轴承内外圈轴向相对距离，a2j为轴承内外圈径向相对距离，b＝fi+fo-1，fi、fo为内、外滚道曲率半径系数，db为球体直径，α0为轴承初始接触角，δa为轴承内外圈轴向总相对位移，δr为轴承内外圈径向总相对位移，δc为轴承内圈离心位移。

如图3所示，o-xyz为主轴坐标系，o1-x1y1z1为轴承一坐标系，o2-x2y2z2为轴承二坐标系。

对轴承一：

δa＝δz+riθxcosψj-riθysinψj+ua

对轴承二：

δa＝-δz-riθxcosψj+riθysinψj+ua

式中，δz为轴承内外圈轴向相对位移，ri为内圈曲率中心圆半径，θx为轴承内外圈x向相对角位移，θy为轴承内外圈y向相对角位移，ψj为轴承滚动体位置角，ua为轴承内外圈轴向相对热位移。

轴承内外圈径向总相对位移：

式中，δx、δy分别为轴承内外圈x向和y向相对位移，θr为最大变形处的滚动体与位置角为0的滚动体之间的夹角，且ur为轴承内外圈径向相对热位移。

轴承内圈的离心位移δc：

式中，ρi表示轴承内圈密度，ωi表示轴承内圈转速，ei表示轴承内圈弹性模量，μi表示轴承内圈泊松比，d表示轴承内圈直径，di表示内圈沟道直径。

由图2可得到角接触球轴承沟道接触的变形几何相容方程为：

(a1j-x1j)²+(a2j-x2j)²-[(fi-0.5)db+δij]²＝0

x1j²+x2j²-[(fo-0.5)db+δoj]²＝0

式中，x1j、x2j分别为球体中心与外沟道曲率中心的轴向距离和径向距离，δij、δoj分别为内、外圈与球体接触的弹性变形趋近量，db为球体直径。

在通过球体中心和轴承轴线的平面上，角位置ψj处的球体所受载荷如图4所示，其中λij和λoj分别表示内、外圈沟道控制系数，对于高速情况，可假设球体所受陀螺力矩完全被球体与外圈沟道的接触摩擦力抵消，这时，取λij＝0,λoj＝2。根据hertz点接触变形关系及球体受力平衡关系可得：

式中，kij、koj分别为内、外圈滚道与球体接触的载荷变形系数，αij、αoj分别为轴承内、外圈与球体接触角，mgj、fcj分别为第j个球体受到的陀螺力矩和离心力。

为了求出稳定状态时轴承的各向位移值，还需要建立轴承所受载荷的平衡方程，假设轴承受到联合载荷f＝{fx,fy,fz,mx,my}的作用，根据图3所示，在主轴坐标系下，轴承载荷与位移的关系可用下式表达：

定压预紧时：

对轴承一：

对轴承二：

式中，fx、fy、fz、mx、my分别为轴承受到的x、y、z方向的外载荷和x、y方向的外力矩，j为轴承球体序数，z为轴承球体个数。

定位预紧时，轴承载荷平衡方程只需去掉定压预紧载荷平衡方程中的z向方程即可。

电主轴工作过程中，由于轴承滚动体与内外圈的摩擦作用，轴承将产生一定的热量，下面给出轴承产热量的计算方法。

对于位置角ψj处的球体有：

式中，mij、moj分别为球体与内、外滚道的摩擦力矩，f0、f1是与轴承类型和润滑、负载有关的系数，η0为润滑油的运动粘度，ωcj为第j个球体的公转角速度，qij、qoj分别为球体与内、外滚道的接触应力，qimax、qomax分别为球体与内、外滚道接触的最大应力。

对于外滚道控制，仅需考虑球体与内滚道的自旋摩擦力矩，即：

式中，msij为球体与内滚道的自旋摩擦力矩，μsi为球体与内滚道的摩擦系数，aij为球体与内滚道的赫兹接触椭圆长半轴，σij为球体与内滚道赫兹接触椭圆的第二类完全积分。

故内、外滚道接触区的摩擦发热量分别为：

hij＝ωcj.mij+ωbj.msij

hoj＝ωcj.moj

式中，hij、hoj分别为球体与内、外滚道接触区的摩擦发热量，ωbj为球体的自旋角速度。

步骤二：建立电主轴热传递网络模型，分析不同类型热阻的计算表达式

电主轴在运转过程中，存在着三种基本的热传递方式：热传导、热对流、热辐射。通常在电主轴系统中装有冷却系统，因此系统内部的热辐射对整个电主轴温度场的分布影响很小，可忽略不计。若将电主轴系统用有限个节点通过传热热阻连接起来，则构成了一种基于热阻的电主轴热网络模型。按系统换热方式的不同，可将传热热阻分为热传导热阻和对流换热热阻以及零部件结合面之间的接触热阻。针对电主轴系统主要零部件的结构，将电主轴轴承座或空心转轴简化为若干个圆筒，实心转轴则简化为若干个圆柱，轴承内外圈、电机定转子简化为圆筒，然后以传热学理论来计算各种类型的传热热阻，本实例中简化后的电主轴二维结构图如图5所示，下面介绍各种类型热阻的计算方法。

1.热传导热阻

热传导热阻可分为轴向、径向热阻。根据fourier定理，圆筒和圆柱的轴向热传导热阻r1为：

式中，l为圆筒或圆柱轴向长度，k为圆筒或圆柱热导率，sa为圆筒或圆柱轴向截面积。

圆筒径向热传导热阻r2为：

式中，d2、d1分别为圆筒外径和内径。

圆柱径向热传导热阻r3为：

球体热传导热阻r4为：

式中，kb为球体热导率。

油脂润滑时内、外滚道的润滑脂热传导阻抗r5和r6分别为：

式中，kl为润滑脂热导率，do为轴承外圈沟道直径，bi、bo分别为轴承内、外圈宽度。

2.对流换热热阻

对流换热主要表现在：轴承座、转轴端部与外界环境，定子、转子端部与周围环境，定子与冷却水，定子内表面、转子外表面与冷却空气。根据fourier定理，圆筒和圆柱的轴向、径向对流换热热阻r7为：

式中，h为对流换热系数，s为圆筒或圆柱与外界对流换热面积。

下面介绍系统各处换热系数的计算方法。

轴承座与外界环境发生自然对流换热和辐射换热，其复合换热系数h1为：

式中，nu1为轴承座表面与空气对流换热的努赛尔数，λe为环境介质热导率，de为换热表面特征直径。

其中，g为重力加速度，βe为环境介质热膨胀系数，th为轴承座表面温度，te为周围环境温度，pr为环境介质普朗特数，ηe为环境介质的运动粘度。

冷却套中的流体与周围壁面之间及定转子气隙冷却空气与定转子壁面之间的换热系数h2为：

式中，nu2为流体与壁面对流换热的努赛尔数，λf为流体热导率，df为通道特征尺寸。

其中，nu2＝0.0225re1^0.8pr^0.3，re1为雷诺数，且uf为流体的运动速度，ηf为流体的运动粘度。

流体运动速度qf为流体体积流量，af为通道截面积。

对于冷却水套，通道特征尺寸cf为通道截面周长；对于定转子气隙，df为定转子气隙长度。

电主轴系统旋转部分与周围环境流体对流换热系数h3为：

式中，deq为旋转部分当量直径，nu3为流体与旋转壁面对流换热的努赛尔数，

其中，re2为雷诺数，且ω为旋转部分转速。

3.接触热阻

宏观上，相互结合的两个机械结构表面，如轴承与转轴的配合面，看似完全结合在一起，但实际上的真实接触只发生在接触面的有限个微凸峰上。当有热量通过结合面时，存在于界面间隙中的流体给传热过程带来了一定的附加阻力，称之为接触热阻，本实例中考虑了轴承内圈与转轴、轴承外圈与轴承座及转子与转轴等热源附近结合面的接触热阻。结合面的接触热阻表示为：

式中，hc为结合面接触热导系数；a为结合面名义接触面积。

考虑微凸峰、界面流体介质的热传导，忽略气隙间的辐射换热，接触热导系数表示为：

式中，lg为结合面空隙厚度，k1、k2分别为两个配合零件的热导率，kf为界面间隙介质的热导率，a^*为结合面无量纲实际接触面积，且ac为实际接触面积。

4.电主轴热网络模型

按照系统节点布置原则：以轴承、电机等径向几何中心平面为分割线，沿径向布置一系列节点，其中结合面两端各布置一个节点，同时保证各个径向节点沿轴向方向对齐，以便建立系统热平衡方程。划分好的电主轴热网络节点如图6所示，建立的电主轴热网络模型如图7所示，图7中的te代表环境温度，tcw代表冷却水温度，tca代表冷却空气温度，其余各温度则为待求的节点温度。

步骤三：建立电机损耗模型，分析电机产热

本实例中的电主轴为永磁同步电机电主轴，工作过程中，定子损耗包括定子绕组上的铜耗和铁耗，转子铁耗很小可忽略不计，考虑转子涡流损耗。

定子绕组上的铜耗pcu：

pcu＝mi²r

式中，m为定子绕组相数，i为定子绕组相电流，r为定子绕组阻值。

根据经典铁心损耗分离计算方法，定子铁心损耗pfe由磁滞损耗ph、涡流损耗pe以及附加损耗pexc三者之和确定，即：

pfe＝ph+pe+pexc

其中，

式中，ms为定子质量，kh为磁滞损耗系数，ke为涡流损耗系数，kexc为附加损耗系数，f为电机工作频率，bm为最大磁通密度，α为无量纲磁滞损耗计算参数。

转子涡流损耗pr：

pr＝ρl∫∫j²ds

式中，ρ为永磁体的质量密度，l为永磁体导路长度，j为永磁体的感应涡流电流密度，s为永磁体导路横截面积。

步骤四：分析不同安装方式下的轴承套圈热位移变化情况

根据传热与热变形基础理论，等截面杆一维变温位移微分方程式为：

式中，u为轴向位移，x为轴向长度，t为温度。

当杆件无内热源稳定导热时，可得：

式中，c1、c2、c3、c4为积分常数，c3和c4可由温度边界条件求得，c1和c2可由位移边界条件求得。

根据图5可知，轴承一、轴承二和轴承六正向安装，轴承三、轴承四和轴承五反向安装，下面介绍不同安装方式下的轴承套圈热位移计算方法。

对于正向安装，轴承套圈轴向相对热位移ua：

ua＝δs-δh-ub(sinαi+sinαo)

对于反向安装，轴承套圈轴向相对热位移ua：

ua＝δh-δs-ub(sinαi+sinαo)

式中，δs为转轴与轴承配合处的轴向热位移，δh为轴承座与轴承配合处的轴向热位移，其中，δs和δh可根据前面介绍的杆件无内热源稳定导热计算公式求得，ub为轴承球体热膨胀位移，且ub＝βbδtbdb，βb为球体热导率，δtb为球体温度变化量，αi和αo分别为球体与轴承内外圈的接触角。

正反装轴承的轴承套圈径向相对热位移uc的计算方法相同，为：

式中，βi为轴承内圈热导率，δti为轴承内圈温度变化量，βs为转轴热导率，δts为转轴与轴承配合处的温度变化量，μs为转轴泊松比，βh为轴承座热导率，δth为轴承座与轴承配合处的温度变化量，μh为轴承泊松比，do为外圈沟道直径。

步骤五：拟定分析计算流程，选定算法进行求解

为准确掌握高性能电主轴热动态参数的变化规律，设计了一个考虑力热耦合效应的高性能电主轴热动态特性的分析计算流程，如图8所示。其中，采用newton-raphson迭代法求解轴承拟静力学模型和轴承力热耦合模型，采用基尔霍夫电流定律(kcl)求解电主轴热网络模型，即在任一瞬时，流向某一节点的热流量之和恒等于由该节点流出的热流量之和，对当前热网络模型来说，设系统共有n个温度节点，且第i个节点有mi个热交换途径，应用kcl，则有：

式中，i,j＝1,2,…,n，ti表示第i个节点的温度，rij表示第i个节点与第j个节点之间的热阻，qi表示第i个节点自身产生的热流量，视其为流入热流量，并取流入热流量方向为正，流出热流量方向为负。将所有节点的热流平衡方程组装起来便构成了电主轴系统热流平衡方程，即：

{q}+[g]{t}＝0

上式即为电主轴热网络数学模型表示形式，其中{q}为系统生热矩阵，由轴承产热和电机产热计算求得，[g]为系统热阻矩阵的倒数，即热导矩阵，确定了{q}和[g]，便可得到系统温度分布矩阵{t}。

整个计算过程如下：输入已知参数后，首先通过求解轴承拟静力学模型，得到轴承接触区的摩擦生热量，再根据电机损耗模型，计算得到电机产热量。将计算得到的各部分热阻值和轴承产热值、电机产热值代入电主轴热网络数学模型中，获得电主轴温度分布，由此可计算得到轴承与转轴和轴承座配合处的热位移，将得到的热位移值代入轴承力热耦合模型中并求解，若是收敛，则输出计算结果，否则便根据各个轴承的热位移情况修正各个轴承的内部几何关系，之后再对轴承拟静力学模型进行求解，如此反复，直到各模型达到设定的收敛精度后，方可输出最终结果。

步骤六：根据计算流程求解并输出计算结果，得到电主轴热动态参数随主轴转速的变化规律。

本发明以一种永磁同步电机电主轴为算例，根据拟定的计算流程，求解得到电主轴各热动态性能参数，包括电主轴节点温度以及轴承支承刚度等随转速的变化规律。本算例中的轴承选用角接触混合陶瓷球轴承，具体参数示于表1，轴承座材料为ht300，转轴材料为40cr，前轴承预紧方式为定位预紧，后轴承预紧方式为定压预紧，预紧载荷均为50n，环境温度为25℃，冷却水入口温度为20℃，螺旋冷却水套中的冷却水体积流量为4l/min，电机定转子气隙间的轴向冷却空气体积流量为1.5l/min，主轴工作转速为1000～8000r/min。

图9a和9b示出了本算例中转速在1000～8000r/min范围内的电主轴轴承球体和电机定转子的温度随转速的变化曲线，从图中可以看到，随转速上升，各轴承球体温度逐渐上升，在转速范围内最大上升约3.2℃，且前轴承球体温度上升速率高于后轴承，而定子温度和转子温度随转速上升逐渐下降，这是因为冷却水套中的冷却水流量不变，定转子发热量也基本不变，而定转子与周边环境的换热系数逐渐增大，所以定转子温度逐渐降低。

图10a-10c示出了本算例中转速在1000～8000r/min范围内的电主轴轴承各向支承刚度变化曲线，从图中可以看到，电主轴前轴承即轴承一、2、3、4的各向支承刚度均随转速升高而逐渐增加，其中轴承一、2的刚度增加趋势较轴承三、4明显，这是因为轴承一、2的球体温度上升速率大于轴承三、4的球体，使得轴承一、2的内外圈与球体的接触载荷增加更快，故刚度增加更快。后轴承即轴承五、6的轴向刚度和角向刚度随转速上升略有下降，径向刚度则略有上升，这是因为定压预紧下，轴承内外圈允许有轴向位移，而内圈接触载荷基本不变，但接触角增大，外圈接触载荷增加，使得外圈法向接触刚度增加，但内圈法向接触刚度基本不变，刚度的串联作用使轴承径向刚度有所增加，但增加幅度较小，而外圈接触角的减小使得轴承轴向刚度和角向刚度均有一定幅度的下降。

表1轴承参数

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。