吴 磊，于慎波，于言明，窦汝桐，翟凤晨

(1.沈阳工业大学机械工程学院, 沈阳 110870；2.中国邮政集团公司辽宁省机要通信局, 沈阳 110000)

0 引言

永磁同步电主轴作为高速机床的核心部件，将电机和主轴组成一个整体，实现了机床主轴系统的高效传动[1]。随着永磁材料的不断创新与发展，使得永磁同步电主轴在航空航天、机械加工和精密制造等领域得到了越来普遍的应用[2]。永磁同步电主轴转子细长，且转速较高，和其它转子相比，其转子系统接近临界转速的机率更大。另外，由于外部环境、负载的变化、加工误差等原因会导致转子发生动态偏心[3]，进而产生不平衡磁拉力。在正常运行过程中，电主轴转子固有频率和转子不平衡力导致的强迫振动频率相耦合，当转速接近该阶模态下的临界转速时，主轴将会发生多重共振，导致噪声过大，严重时，甚至导致转子不稳定或损坏[4]。为了避免转子在运转中产生共振，影响电主轴正常运行，需要对转子偏心引起的不平衡磁拉力以及转子临界转速和稳定性进行相关的动力学分析[5-6]。

文献[7-9]在转子系统的不平衡响应分析、临界转速和稳定性等领域做了许多的研究，为旋转机械的发展奠定了理论基础，但主要集中在水、汽轮机，对永磁同步电主轴的研究涉及很少。目前求解临界转速常采用传统Prohl传递矩阵法、Riccati传递矩阵法和有限元法等[10]，其中Prohl传递矩阵法应用在各项同性转子系统，Riccati应用在各项异性转子系统，但由于复杂的外部环境以及算法本身的缺陷，计算时会出现根的溢出和丢根的现象，导致临界转速计算不精确。有限元法计算转子系统临界转速时，对于复杂转子由于建模复杂，占用计算机内存较大，计算速度较慢，在一定程度上限制了其应用范围。因此，如何提高解析法计算转子临界转速的精确度以及进行机电耦合的分析，成为了转子动力学研究热点。

为了解决上述问题，本文中提出了一种用加时间因t的方法构建转子传递矩阵，运用了双重QR法[11]求解传递矩阵特征值，进而可求得转子系统的临界转速，既提高了计算速度，又解决了根的上下溢出问题。同时运用编程计算了转子运行中由于动态偏心而导致的不平衡磁拉力，使其与机械系统进行耦合，分析其对动力学特性的影响，对研究转子系统动力学特性有着重要的参考价值。

1 永磁同步电主轴临界转速计算

1.1 转子离散建模与传递矩阵的建立

根据文献[12]中节点划分原则，应用集总质量法对所求转子进行离散处理[8]。本文中，永磁同步电主轴转子如图1所示，将其简化为25个轴段，26个节点(节点7和19为轴承支撑点)并在对应节点附加上相应的质量和转动惯量。

图1 电主轴转子总成模型

对于离散后的第i个单元，将i和i+1分别定义为其左右两端的截面编号，根据转子受力特点，考虑轴承支承为各项异性，将截面的x方向和y方向弯矩、剪力、位移和转角组合成该截面的状态向量，如式(1)所示：

(1)

式中，Mx、My、θx、θy、x、y、Qx、Qy分别为水平和垂直方向的弯矩、转角、挠度和剪力。

引入以复变量表示的振动量w=τ+σi，其中，w为复涡动频率；实部τ为衰减指数；虚部σ为阻尼圆周率。根据圆盘的受力分析，可得出加时间因子t的刚性薄圆盘的传递矩阵，如式(2)所示：

(2)

根据轴段左右截面状态向量关系，可以得出等截面弹性轴的传递矩阵，如式(3)所示：

(3)

将状态向量{Z}i分别写成力向量{Z}g和位移向量{Z}h，根据边界条件以及转子系统各截面状态向量的关系，可得出方程(4)：

(4)

其中，[Q] =[B]K[D]K[B]K-1[D]K-1…[B]1[D]1。

由方程(4)可得：

[Q12]{Z}h,1=0

(5)

由于{Z}h,1≠0时，上式有非零解的条件为：

|Q12|=0

(6)

由于式(6)中的元素为多项式，运用代数余子式法将其展开，频率方程实际上是关于w的8K阶一元高次方程，即：

e0+e1w+e2w2+···+e8K-1w8K-1+w8K=0

(7)

求方程(7)的根可以转化为求其对应友矩阵全部特征值的问题。

电主轴转子为细长阶梯轴，为了准确模拟实际转子工况，需要离散足够的节点数，本文中离散节点数K=26，方程最高阶会达到208，求解过程中会出现根的溢出，求解难度特别大。文中，为了缩小求解过程中大数与小数之间的差异，防止根的溢出，提出加时间因子t的方法建立转子系统传递矩阵。计算过程中其他参数与时间因子相乘，转速与时间因子相除，运用双重步QR法编程求解出复频率的实部与虚部，计算完成后将虚部乘以时间因子得到临界转速。

1.2 解析法计算永磁同步电主轴临界转速

本文中以额定功率9.5 kW永磁同步电主轴转子为例，取时间因子t等于15 000构建传递矩阵，取3000 r/min为转速步长，运用编程的方法求解出复频率并画出坎贝尔图。同频线与各阶涡动速度曲线的交点即是转子系统临界转速，如图2所示。从图中可以看出，随着转子自转速度逐步升高，一阶正向和反向涡动速度相近且变化很小；二阶反向涡动速度逐步减小，正向涡动速度逐步增大。

图2 解析法计算转子坎贝尔图

1.3 有限元法计算永磁同步电主轴临界转速

利用ANSYS计算临界转速时，运用实体模型建模复杂，且计算结果精度达不到要求。因此本文运用有限元法时，以MASS21、BEAM188和COMBI214分别模拟质量点、弹性轴和轴承单元。同时质量点单元、轴承单元分别施加旋转约束和全约束，并充分考虑了陀螺效应的影响，通过编写ANSYS APDL命令流的方式建立二维有限元模型如图3所示。

图3 ANSYS二维质量点单元模型

本文用QR阻尼法进行模态提取，利用PLCAMP命令画出轴承-转子系统的CAMP-BELL图，如图4所示。

图4 ANSYS有限元计算转子坎贝尔图

1.4 解析法和有限元法计算临界转速结果对照

各阶临界转速为同频线与正向涡动速度的交点。转子临界转速计算结果对比如表1所示。由表1可知，解析法与有限元计算结果对照，误差均在3.2%以内，验证了本文中解析法计算临界转速的有效性。

表1 临界转速结果对比 (r/min)

2 稳定性分析

稳定性裕度体现在系统自身的修复能力，也就是能够回到平稳状态的能力，度量方法主要有对数衰减率、系统阻尼和系统抗定常干扰界限值法，其中对数衰减率是以振幅衰减的快慢程度为指标，对系统受干扰后恢复到平衡状态过程的度量[8]。

本文以对数衰减率为依据判定轴系的稳定性裕度。由本文第一部分解析法得到复频率w=τ+σi，则第i阶特征值对应的对数衰减率：

(8)

γi越大表明系统稳定程度越高，当γi<0时，表明系统处于失稳状态。

计算得到的轴系对数衰减率曲线如图5所示。图中表明，一阶反向、二阶正向和三阶反向对数衰减率呈逐渐下降趋势，其中三阶反向对数衰减率最小且有加速下降的趋势，最易引起系统失稳，设计时应尽量避开三阶反向涡动临界转速。

图5 永磁同步电主轴转子对数衰减率

3 不平衡磁拉力

3.1 不平衡磁拉力的计算

通常气隙磁场分布均匀，转子不承受不平衡磁拉力，但由于转子偏心导致径向气隙磁场畸变，会产生不平衡磁拉力，使得转子受力不均匀。传统公式计算不平衡磁拉力时，将模型简化为理想状态，计算结果不够精确。因此本文基于Ansoft有限元软件对某型号8极36槽永磁同步电主轴转子动态偏心引起的不平衡磁拉力进行分析。电主轴的基本参数如表2所示。

表2 永磁同步电主轴部分参数

电主轴有限元模型建立完成并求解后可提取气隙磁密，然后由式(9)求出不平衡磁拉力[13-15]。

(9)

其中，Fx、Fy分别是不平衡磁拉力在X和Y方向上的分量；l是转子硅钢片的轴向长度；rp是不平衡磁拉力积分半径。

转子偏心后，径向畸变气隙磁场可看作由多个同心机磁场叠加而成，且切向气隙磁密对径向不平衡磁拉力影响极小，可忽略不计。因此可以将气隙积分圆周离散成N个相等的子区域，求得各区域径向畸变气隙磁场产生的不平衡磁拉力然后进行叠加，所以式(9)可以表示成式(10)的形式。

(10)

式中，δiω、δiρ分别为第i个区域的起始和结束角度。

则不平衡磁拉力的合力为：

(11)

为了防止转子和定子出现扫膛现象，转子偏心距不应大于平均气隙的10%[16]。本文中电主轴气隙总长为1.6 mm，因此偏心距不应超过0.16 mm。根据此原则，以无偏心和动态偏心0.16 mm为例，将模型积分圆周划分为1000个区域，分别提取其径向气隙磁密，如图6、图7所示。

图6 无偏心径向气隙磁密

图7 偏心距0.16 mm径向气隙磁密

应用式(10)，利用编程求出不平衡磁拉力在X和Y轴的分力，如图8所示。图中可看出，当无偏心时，X和Y方向不平衡电磁力均在0 N附近上下波动，且浮动很小；当偏心距为0.16 mm时，X和Y方向上均产生了随时间周期性变化的幅值为230 N的不平衡磁拉力。

图8 不平衡磁拉力在X和Y轴的分力

应用上述方法，可求出偏心距分别为0.032 mm、0.064 mm、0.096 mm和0.128 mm时在X和Y方向的不平衡磁拉力，由式(10)、式(11)可求得不同偏心距下不平衡磁拉力在Y方向分量Fy和在X方向分量Fx随圆周的变化以及不平衡磁拉力的合力，如图9、图10所示。

图9 Fy随Fx在圆周内的变化

图10 不平衡磁拉力合力随时间变化

由图9可看出，相同偏心距下，不同时刻的X和Y方向的不平衡磁拉力在圆周上组成了一个同心圆，且圆的半径即为该偏心状态下的不平衡磁拉力，这是由于动态偏心时最小气隙随着时间在变化。由图9和图10可得到不同偏心距下不平衡磁拉力的数值如表3所示，且不平衡磁拉力随着偏心距的增大呈线性增加。

表3 不同偏心距下的不平衡磁拉力

3.2 机电耦合

当电主轴转子在高速旋转时，不平衡磁拉力被引入机械系统引发共振[17]，严重情况下会造成定转子接触，刮蹭永磁体。因此在动力学分析时，考虑不平衡磁拉力的影响十分必要。将不同动态偏心下产生的不平衡磁拉力与转子机械系统进行耦合后，临界转速的变化如图11所示。

图11 不平衡磁拉力对临界转速的影响

通过分析图11可以看出，考虑不平衡磁拉力时的临界转速低于无不平衡磁拉力时的临界转速。随着不平衡磁拉力的增加，一阶临界转速由23 657 r/min降低到23 339 r/min，降幅为1.4%，二阶临界转速几乎无明显变化。

4 结论

(1)运用加时间因子t的方法构建转子传递矩阵，同时运用双重步QR法计算转子系统临界转速，既提高了计算精度又解决了根的上下溢出问题。同时将解析法与有限元法计算的各阶临界转速进行对比，误差均小于3.2%，验证了解析法的有效性。

(2)通过分析对数衰减率曲线，发现三阶反向涡动衰减率曲线最低，且下降速度有逐渐加快的趋势。说明三阶反向涡动稳定性裕度最小，最易引起转子系统失稳。

(3)分析了动态偏心工况下的不平衡磁拉力，结果表明，不平衡磁拉力随着偏心距的增大呈线性增加。随着不平衡磁拉力增加，一阶临界转速将下降，二阶临界转速几乎不受影响。