一种高效的电主轴系统复频率计算方法

毛文贵，刘桂萍，韩 旭，郭维祺，林禄生

（1．湖南大学 机械与运载工程学院，长沙 410082；2．湖南工程学院 机械工程学院，湘潭 411101）

随着当代高精密加工的飞速发展，对数控机床加工性能的要求也越来越高，而作为数控机床主要功能部件的电主轴单元是其性能提高的关键。电主轴以液体滑动轴承为支撑，将机床主轴和高速电动机功能从结构上融为一体，实现了变频电机和机床主轴之间的“零传动”。但高速电主轴技术还待改进，结构设计、工艺、甚至控制问题都未得到有效解决［1－3］。液体滑动轴承支承的高速电主轴可以避免滚动轴承较大的振动与噪声、磁力轴承过于复杂的控制系统和气体轴承较低的承载能力。但液体滑动轴承油膜的动特性对整个电主轴系统的动特性有很大影响，因此当前高速电主轴技术研究的一大热点就是对滑动轴承油膜动特性的研究。滑动轴承油膜动特性通常可用8个系数来描述。Lund［4］就采用4个刚度系数和4个阻尼系数建立了挠性转子系统模型，揭示了油膜涡动失稳的物理本质。对电主轴各向异性轴承－转子系统进行动力特性分析时，往往要同时考虑转子在径向两平面内的耦合运动以及轴承的阻尼，此时，系统内各个振动量之间发生了相位差，振动幅值随时间增加或减小。因此，振动量为一复数，特征值为复频率，它反映了轴承－转子系统的稳定性，通常作为评价转子性能的重要依据。Riccati传递矩阵法［5］因可消除奇点的干扰而在转子动力学的振动计算中被广泛采用。Murphy等［6］将复变量引入传递矩阵法，计算出系统的复频率。陆颂元［7］对传递矩阵－多项式方法算法作了改进，只用一次递推就求得特征多项式，大大节省了计算时间。故目前数值法求解复频率一般采用传递矩阵一多项式算法［8］。但电主轴不同于一般主轴，电主轴的转子电流和气隙磁场相互作用产生电磁转磁，使转子沿旋转磁场方向转动，所以在对电主轴进行动力学分析时不能忽略转子部分受到的电磁转矩。电主轴是电机转子与主轴过盈热装一起的轴，在进行动力学分析时应将其作为一个整体，电机的自重不能忽略。本文对Riccati传递矩阵法进行改进，获得应用于电主轴各向异性轴承－转子系统的传递矩阵，通过传递矩阵法直接获得特征多项式方程隐式表达，即得到转子系统的频率方程式，由于这是个具有复数自变量S的复数方程式，因而在计算时需对衰减指数λ和阻尼圆频率Ω两个变量作“平面域”扫描，对此问题至今没有简便的计算方法。目前可以采用的是抛物线法（Muller法）求复方程根［9］，但Muller法中设定的起始试算频率不同，求得各个根的先后次序并没有一定的规律。求解所有的根对于大自由度系统而言工作量太大而不实用，而工程中我们主要关心产生共振的临界转速，一般只需知道前几阶。本文由幅角原理［10］判定系统共振区域中根的个数，采用双种群进化策略［11］求解电主轴系统特征多项式的根，可快速有效地获得电主轴系统的复频率。

1 Riccati传递矩阵法的改进

1.1 电主轴系统的模块化处理

电主轴系统被离散为N个集总圆盘，其间用N－1个无质量的弹性轴段相连，由于考虑了油膜力及其他各种阻尼，支承就是各向异性的。所以在X－Z，Y－Z平面（Z为轴向）同时弯曲，则状态向量的维数是8，为提高传递矩阵的数值稳定性，把状态向量中的元素分为两部分，一部分为对应于初始截面状态向量中具有零值的元素组成，另一部分为其余互补元素组成。考虑转子的材料与主轴不同，在计算用于传递参数的转动惯量等参数时会造成一定的麻烦和误差，转子的重力等效为作用于转子轴段上的外力。同时电主轴转子部分受到的电磁转矩不能忽略，等效成一个集中转矩施加于转子轴段。因此，电主轴的传递矩阵法须考虑外力的加入，为使传递矩阵具有通用性并方便建立整体矩阵，本文将状态向量的维数扩为10。状态向量变为

式中：Mx，My为弯矩，Qx，Qy为剪力，X，Y为位移，θx，θy为转角。向量中各元素均为复变量。如y＝Yest，其中Y＝Yc＋i Ys称为该振动量的复振幅；Yc，Ys分别称为它的余弦分量和正弦分量，S＝λ＋iΩ称为复频率，λ为衰减指数，Ω为阻尼圆频率。

对于图1中集总圆盘i和轴段i左右两端进行受力分析得：

式中：I是单位矩阵，O是零矩阵，上标L、R分别表示集总圆盘和轴段的左端或右端。

式中：m为各节点的质量；kxx，kxy，kyx，kyy，Cxx，Cxy，Cyx，Cyy为支承轴承的刚度和阻尼系数；Jd，JP为直径和极转动惯量；ω为自转角速度；V为考虑轴段的剪切影响，V＝6EJ／aGAL2；a为截面系数（实心圆轴为 0．886，空心圆轴为2／3）；G为剪切弹性模量；E为弹性模量；A为截面面积；L为轴段长；J为轴段的截面惯性矩。

1.2 特征多项式

引入Riccati变换并经推导，展开式（2）可得：

参考 Riccati推导［5］和陆颂元的多项式法［7］，引入表征同一截面从位移向量到力向量的转换关系矩阵B，在每一截面上有：

式（4）式代入（3）式整理得：

引入

其中：Ti－1＝Ui－1＋Ri－1

整理可得状态向量之间的递推关系式：

引入边界条件，如具有N个节点的转子，两端自由，则左节点的左边及右节点的右边两截面力向量为零，可得：＝0＝0，由BL1＝0

可得＝0且C1＝I，D1＝W1，由＝0可得：

从i＝2起反复使用式（7），递推到 i＝N得到DN（S）特征多项式。采用了4×4阶矩阵代替了10×10阶矩阵，从而机时进一步减少，具有计算速度快，占用内存少。解此特征多项式便得复频率S＝λ＋iΩ，其中实部为一系统阻尼，反映转子轴承系统的稳定性。对数衰减率为δ＝－2πλ／Ω。对数衰减率大于零则系统稳定。特征多项式的最高次数为整个系统的自由度数。对于节点数为N的转子，有N×4对共扼复根。

2 高效电主轴系统复频率计算方法

当轴承是各向异性轴承时，S是复数，则特征多项式是复函数，当转子节点数较多时，易出现数值溢出问题，虽然张才稳提出特征多项式系数放大gi倍的方法使系数处于机器数系所能表达的范围内［12］。g的取值没有界定，g的寻找需要大量的时间，而且增加了程序的复杂性。现实中并不实用。工程中我们只对系统的某些低阶特征值和特征向量感兴趣，并不需要求解所有的根，对于临界转速，主要是了解与工作转速存在共振可能的那几阶临界转速。因此计算系统在某一范围内的部分特征值和特征向量便是一种非常实用而又经济的方法。Muller法求解复函数，对初值选取较敏感。而且复数的求解不同于实数，由于设定的起始试算频率不同，求得各个根的先后次序并没有一定的规律，也不能确保在共振区域中所有的根都已经被找到。因此能了解在求解的范围里的根的个数才能确保所有的根已经被找到，而不会漏根。以前实际计算中，采用先不计算阻尼和支承的各向异性，获得剩余量的曲线，初步了解此系统的固有频率的分布，然后求解复频率［9］，这样的过程过于复杂和耗时。本文基于复函数方程根的分布理论－幅角原理［13］判定共振区域中根的个数，在相应的区域提出基于双种群进化策略求解特征多项式DN（S），获得电主轴系统的复频率。其流程如图2所示。

图2 高效电主轴系统复频率计算方法流程Fig．2 A high efficient method for calculating complex frequency of the motorized spindle

2.1 复函数方程根的分布理论

定理 幅角原理［13］

若f（z）在正向简单闭曲线C上每点都非零解析，在C的内部是亚纯函数，那么f（z）在C内零点的个数等于1／（2π）乘以当 z沿 C的正向绕行一周 f（z）的幅角的改变量。

推论 设f（z）＝u＋iv在简单闭曲线C上和在C内解析，且在C上不等于零，点z0＝x0＋iy0沿C的正向绕行一周，设向量（u，v）作正方向旋转的次数为N1，作负方向旋转的次数为N2，那么在封闭曲线C内f（z）＝0的根的个数N＝N1－N2。

通过幅角原理对所定义的区间进行根的个数的确定。在本文中共振区域内根的个数为电主轴的特征多项式的值在原点右边的圈数。

2.2 双种群进化策略寻根

进化策略是一种全局优化搜索算法，更新进化很快，特别适合函数优化［14］。在确定了编码方案，适应度的计算及相应进化算子后，以“优胜劣汰”的方式，不断逼近最优解。算法中主要用变异算子作为进化算子，变异算子主要有高斯变异和柯西变异，较大的高斯变异算子可以保证种群具有较好的探索能力，较小的高斯变异算子可以保证种群在局部具有良好的搜索能力。柯西分布因有一条很长的尾巴，容易产生一个远离原点的随机数，比高斯变异产生的随机数具有更宽的分布范围，这样有可能很快跳出局部极小的区域。双种群进化策略使得进化分别在两个不同的种群间并行进行，两个种群中采用不同的变异算子，种群1的目标变量与决策变量均采用高斯变异，种群2的目标变量采用柯西变异、决策变量采用高斯变异。种群1和种群2均采用（μ，λ）选择策略，得出优良的个体作为进化的结果。双种群方法具有计算精度高、自适应强而在求复函数方程中得到较好的应用。

根的个体由目标变量和标准差两部分组成，每个个体有2个分量，分别代表根的实部和虚部。种群1和2的解代入特征多项式中，获得其函数值，函数值越接近零，则根越接近最优解。规定适应度为最大值，取适应度函数为：

并终止条件为一个很接近1的数ε。当适应度值大于ε时终止，从而求解出一个复频率。

2．3 综合除法降次

用双种群进化策略算法求出特征多项式一个根后，考虑转子特征多项式根的共轭性，用综合除法进行降次，方程（8）变为如下：

式中：Sr为已求根；为已求根的共轭根；F为已求出的根的个数。

获得下次计算的特征多项式，可知求一个根后特征多项式可降低二次。

3 计算实例

3.1 电主轴系统

根据上述方法，对图3所示模型转子进行了计算，电主轴系统最高工作转速为1 000 rad／s，电机功率为35．3 kW，主轴段密度 ρ为7 650 kg／m3。弹性模量 E为210 GPa。电主轴上各轴长 Lc、轴内径 dc，、轴外径Dc和圆盘宽度Lp，圆盘外径Dp参数如表1，本模型被分为9个节点，液体滑动轴承支承在第3，5节点上，两支承为相同的圆柱轴承，其参数为：长径比L／D＝0．5，间隙比ψ＝2．0‰，动力粘度 η＝0．018 4 N．s／m2。本文采用熊万里提出的基于动网格模型的计算方法［15］获得液体动静压轴承7个偏心率下的8个刚度阻尼系数。轴承的刚度和阻尼无量纲系数随偏心率的变化关系如图4。利用改进的Riccati传递矩阵法获得主轴段参数和节点总圆盘参数，得到包含复频率的特征多项式的隐式表达。

图3 电主轴各向异性轴承—转子系统模型Fig．3 Anisotropic bearings-rotor systems of motorized spindle

表1 电主轴系统参数Tab．1 Parameters of motorized spindle system

图4 轴承刚度阻尼无量纲系数随偏心率的变化关系Fig．4 Dimensionless stiffness and damping coefficients vs．Journal eccentricity

3.2 高效电主轴系统复频率计算方法

因本模型最高工作转速n为1 000 rad／s，对于在第一临界转速nc1以下工作的转子，应使工作转速小于0．75 nc1，对于在第一临界转速nc1以上工作的转子，应使1．4 nc1＜n＜0．7nc2。由此确定系统关心的临界转速不超过 2 000 rad／s，此系统关心的共振区域定为［0，2 000］。对工作频率所可能激振的共振区域采用幅角原理进行根的个数的确定，如图5所示，可知在［0，2 000］中电主轴系统特征多项式有5个根。用双种群进化策略求解复频率。双种群进化策略中μ＝20，λ＝140，终止条件ε＝0．999，5个复频率的计算结果如表2所示。为了验证结果，采用传统的Muller法求得本模型的特征值。本模型为有9个节点的轴系，共有36个复频率，Muller法求复频率的无序性使得必须求出全部的36个复频率，然后按复频率虚部值从小到大排列，取出前6个复频率，如表2所示。对比表2中两种方法求解的5个复频率，可以看出，共振区间上两种方法的结果误差很小，说明本文所提方法求解的结果是准确的。第五个复频率的虚部小于2 000 rad／s，而第六个复频率的虚部大于2 000 rad／s，说明幅角原理对共振区间里根的个数的判断是对的。从图6显示的适应度值随迭代次数变化的5条曲线图可知迭代次数不超过45代，说明双种群进化策略收敛速度很快。同时采用本文方法只须求5个复频率就能获得工程中所关心的临界转速，而采用传统的Muller法要求36个复频率（本模型求出全部解需要9分钟）。因此采用本文方法可以大大节省工作量。

图5 基于辐角原理获得共振区内复频率的个数Fig．5 Number of complex frequency among resonance region based on argument principle

表2 复频率计算结果及时间对比Tab．2 Comparison of complex frequency and time consumed

图6 适应度值随迭代次数变化的5条曲线图Fig．6 five fitness value and iterative number of times change tendency

4 结 论

（1）本文通过对Riccati传递矩阵法进行改进，考虑电机外力的加入，将状态向量的维数扩为10。采用Riccati分块矩阵递推，对所求的解不赋以具体数值而将其作为未知数保存，一次递推求出以此为变量的特征多项式方程。

（2）针对复频率的平面域扫描求解，传统的Muller法求解根的无序性。依据复函数方程根的分布理论，由幅角原理判定工程中所关心的共振区域内根的个数。

（3）本文中采用双种群进化策略，进化分别在两个不同变异算子的种群中并行进行。避免了大规模范围搜索时易产生的局部极值问题，是一种全局优化搜索方法。计算精度高、自适应强。

（4）以某高速主轴系统为例，采用本文方法对其进行了复频率求解，并与Muller法进行对比，表明本方法计算结果准确，收敛速度很快，只须求几个复频率就能获得工程中所关心的临界转速，工作量大大减小。这对做出可靠的稳定性分析是很重要的。

［1］熊万里，阳雪兵，吕浪，等．液体动静压电主轴关键技术综述［J］．机械工程学报，2009，45（9）：2－18．XIONG Wan-li，YANG Xue-bing，LULang，et al．Review on key technology of hydrodynamic and hydrostatic highfrequency motor spindles［J］． Journal of Mechanical Engineering，2009，45（9）：2－18．

［2］孟杰，陈小安，合烨．高速电主轴电动机－主轴系统的机电耦合动力学建模［J］．机械工程学报，2007，43（12）：160－165．MENG Jie， CHEN Xiao-an， HE Ye． Electromechanical coupled dynamical modeling of high speed motorized spindle’s motor-spindle subsystems［J］． Journal of Mechanical Engineering，2007，43（12）：160－165．

［3］陈小安，单文桃，周进明，等．高速电主轴驱动控制技术研究综述［J］．振动与冲击，2013，32（8）：39－47．CHEN Xiao-an，SHAN Wen-tao，ZHOU Jin-ming，et al．A survey of drive and control technology for high-speed motorized spindle systems［J］．Journal of vibration and shock，2013，32（8）：39－47．

［4］Lund JW．The stability of an elastic rotor in journal bearing with flexible damped supports［J］．ASME Journal of Applied Mechanics，1965，32（4）：911－920．

［5］王正．Riccati传递矩阵法的奇点及其消除方法［J］．振动与冲击，1987，（2）：74－78．WANG Zheng．Singular point and elimination method of Riccati transfer matrix method［J］．Journal of Vibration and Shock，1987，（2）：74－78．

［6］Murphy B T，Vance JM．An improved method for calculating critical speeds and rotor dynamic stability of turbo machinery［C］．Trans of ASME for Power，1983，（5）：141－145．

［7］陆颂元．转子系统复频率的传递矩阵－多项式计算法［J］．应用力学学报，1986，31（1）：73－82．LU Song-yuan．A method for calculating damped critical speeds and stability of rotor-bearing systems［J］．Journal of Applied Mechanics，1986，31（1）：73－82．

［8］Zu JW，Ji Z．An improved transfer matrix method for steadystate analysis of nonlinear rotor-bearing systems ［J］．Transactions of the ASME，APRIL，2002（124）：302－310．

［9］闻邦椿，顾家柳，夏松波，等．高等转子动力学［M］．北京：机械工业出版社，2000．

［10］Beardon A F．Complex analysis［J］．The Argument Principle in Analysis and Topology．Chichester，U．K．：Wiley，1979．

［11］夏慧明，梁华，周永权．用双种群进化策略算法求解复函数方程的根［J］．计算机工程与应用，2008，44（7）：78－81．XIA Hum-ming，LIANG Hua，ZHOU Yong-quan．Novel bigroup evolution strategy algorithm for solving complex functional equation ［J］． Computer Engineering and Applications，2008，44（7）：78－81．

［12］张才稳，瞿红春，陆颂元．转子系统复频率算法中数值溢出问题的研究［J］．振动、测试与诊断，1998，18（3）：206－210．ZHANG Cai-wen，QU Hong-chun，LU Song-yuan．Research on the problem of overflow in calculating the damped critical speeds of rotor-bearing systems［J］．Journal of Vibration，Measurement＆Diagnosis，1998，18（3）：206－210．

［13］Rusu C，Astola J．Note on argument principle［J］．Ieee Signal Processing Letters，2004，11（10）：817－820．

［14］Lee Wei-ping，Chien Wan-jou．A novel differential evolution with co-evolution［J］．Journal of Computers，2011，6（3）：594－602．

［15］熊万里，侯志泉，吕浪，等．基于动网格模型的液体动静压轴承刚度阻尼计算方法［J］．机械工程学报，2012，48（23）：118－126．XIONG Wan-li，HOU Zhi-quan，LULang，et al．Method for calculating stiffness and damping coefficients of hybrid bearings based on dynamic mesh model［J］．Journal of Mechanical Engineering，2012，48（23）：118－126．