许 琦，姚红良，刘子良，闻邦椿

(东北大学 机械工程与自动化学院，沈阳 110819)

由于当前转子系统的密封间隙和轴承间隙越来越小，碰摩故障成为了一种常见的故障，不仅影响转子系统的运行稳定性，而且会引起灾难性的后果。及早的诊断转子碰摩状态信息，对处理碰摩故障有着重要的意义[1]。

目前，学者们对转子碰摩故障问题做了许多研究。在碰摩故障位置诊断方面，目前主要有两种方法：一是基于信号的诊断方法，另一种是模型求解方法。基于信号的诊断方法是采集转子本身发出的振动信号，经过各种信号处理方法，如傅里叶变换、小波变换等将采集的信号分解进行碰摩诊断的方法。如Hall等[2]利用声发射监测技术诊断了转子-定子碰摩故障，并通过同时检测相邻轴承相位延迟估算了碰摩位置；同时Hall等[3]利用Kolmogorov-Smirnov统计特性处理声发射信号诊断不同类型的轴密封碰摩；文献[4]以声发射传感器在不同轴承支承处采集的数据样本为基础，通过小波变换将信号分解成包含特定倍频带的一系列时域信号，采用互关联法(method of cross-correlation)诊断碰摩位置，这种小波变换互关联法提高了互关联法的计算精度，采样频率越高，诊断的精度越高，并通过实验验证了该方法的有效性；文献[5]采用经验模态分解法(EMD)将碰摩故障转子的复杂振动信号分解，得出碰摩信号，从而诊断出转子的局部碰摩故障。

目前基于信号的诊断方法能够快速直接地判断故障的存在，但诊断精度不高，只能估算碰摩故障位置，或者需要较高的条件来提高诊断精度，难以满足需要。模型求解方法是建立碰摩故障转子系统的动力学模型，通过对模型的分析求解，可以得出精确的碰摩位置、碰摩力、严重程度等故障信息，便于对碰摩故障的掌握和解决。张娅等[6]采用有限元动力学模型，同时考虑圆盘轴向和径向碰摩分析了双盘转子系统非线性动力学行为，结合系统响应特征图可以识别系统的碰摩故障；Chu等[7]采用最小二乘法获得多圆盘转子系统的动态刚度，通过动态刚度在转轴碰摩位置和其它位置的变化率不同诊断出了多圆盘转子系统的碰摩位置，同时发现阻尼系数也有类似的特性；Cong等[8]建立了冲击能量模型(IEM)评估转子系统的碰摩故障的存在和严重程度，并通过波、谱和轴心轨迹等实验验证了该模型的有效性；姚红良等[9]计算了转子系统动力学方程的剩余振动量，采用最小二乘法、利用挠曲线诊断出了周期故障力单点故障发生的位置，并能够确定故障力的大小，诊断出系统故障的严重程度，碰摩故障的故障力即为周期故障力，故此方法适用于转子碰摩故障；同时姚红良等[10]根据碰摩故障转子系统中任意两节点之间高次谐波分量之比等于无故障转子系统频率响应矩阵的相关元素之比的关系，提出了诊断转子系统单碰摩故障的方法，并通过数值模拟和实验验证了方法是非常有效的。

目前模型求解方法对于碰摩信息的诊断精度比较高，但大多只是对针对转子系统单碰摩故障，对于双碰摩甚至是多碰摩故障的诊断研究还是比较少的。本文以单碰摩故障诊断理论以及谐波分量的特征为基础，推导出诊断转子系统双碰摩故障的方法，通过碰摩故障前后转子三个测点的振动响应数据就可以分别确定双碰摩故障转子系统的两个碰摩位置，同时根据系统的动力学模型即可求出故障力的数值，并通过数值模拟仿真实验验证了方法的有效性。

1 转子系统的动力学模型

转子系统由弹性轴段单元组成，如图1所示。

图1 转子轴段单元有限元模型

轴段单元的广义坐标为两端节点两个方向的位移和转角，忽略轴段轴向变形，仅考虑弯曲变形和扭转变形，则广义坐标复数表示为

u=[xA-jyAθyA-jθxAxB-jyBθyB-jθxB]T

(1)

式中：x，y和θ分别为节点广义坐标中两个方向的位移和转角。

设转子系统有J个节点，则有J-1个轴段，其动力学方程为：

(2)

式中：M为质量矩阵，C=D+ωG，D为阻尼矩阵，G为陀螺力矩矩阵，ω为转子转速，K为刚度矩阵，矩阵为2J×2J方阵，u为振动响应向量，F为质量不平衡引起的激励矢量。

设转子系统碰摩故障发生在节点L1和L2处，则故障转子系统的动力学方程为：

jPy1)ej(ωt-ξ1)-T2L2-1(Px2+jPy2)ej(ωt-ξ2)

(3)

由式(2)，式(3)，得：

jPy1)ej(ωt-ξ1)-T2L2-1(Px2+jPy2)ej(ωt-ξ2)

(4)

式中：Δu=ur-u。

式(4)即为故障转子系统剩余振动量动力学方程[11]，从式(4)中可以看出，碰摩力(式(4)的右侧)可以被视为外加等效力，质量、阻尼和刚度矩阵保持不变。

2 等效系统谐波分量的特征

将Δu展开成各阶谐波分量和的形式，即：

Δu=Δu0+Δu1ejωt+Δu2e2jωt+Δu3e3jωt+…

(5)

式中：Δui(i=0，1,2,…)为碰摩故障转子系统剩余振动响应向量的第i阶谐波分量。

并将等效力也展开成各阶谐波分量和的形式，即：

(6)

式中：Ai、Bi(i=0，1,2,…)分别为碰摩力的第i阶谐波分量。

根据谐波平衡法，有：

Δu0=K-1(T2L1-1A0+T2L2-1B0)=

K-1T2L1-1A0+K-1T2L2-1B0=

E2L1-1(0)A0+E2L2-1(0)B0

Δu1=[-Mω2+jωC+K]-1(T2L1-1A1+T2L2-1B1)=

E2L1-1(jω)A1+E2L2-1(jω)B1

…

Δul=[-Ml2ω2+jωlC+K]-1(T2L1-1Al+T2L2-1Bl)=

E2L1-1(jlω)Al+E2L2-1(jlω)Bl

(7)

式中：E(jlω)=[K+jωlC-ω2l2M]-1为转子系统的频率响应函数矩阵，l为谐波分量的阶数，E2L1-1(jlω)和E2L2-1(jlω)分别为矩阵E(jlω)的第2L1-1和2L2-1列。

设碰摩故障转子节点i剩余振动响应向量第l阶谐波分量为Δul(2i-1)ejlωt，节点k剩余振动响应向量第l阶谐波分量为Δul(2k-1)ejlωt，由式(7)可得：

(8)

式中：E2i-1,2L1-1(jlω)、E2k-1,2L1-1(jlω)和E2i-1,2L2-1(jlω)、E2k-1,2L2-1(jlω)分别为列向量E2L1-1(jlω)和E2L2-1(jlω)的第2i-1和2k-1行元素，Δul(2i-1)和Δul(2k-1)分别为列向量Δui的第2i-1和2k-1行元素，Al和Bl由下式求得。

(9)

设转子系统的两个任意节点n1和n2，由式(8)，得：

(10)

转子系统各个节点在第l阶的谐波分量Δul为：

(11)

取系统的响应中的前N阶谐波分量，则故障转子系统的故障力分别为：

(12)

式中：H(jlω)=[-Ml2ω2+jlωC+K]，H(2L1-1)(jlω)和H(2L2-1)(jlω)分别为矩阵H(jlω)的第2L1-1和2L2-1行。

3 计算步骤

本方法的具体计算步骤如下：

(1) 分别测得转子系统无故障时的任意三个节点i、，k和h的响应为ui=xi+jyi、uk=xk+jyk和uh=xh+jyh；

(2) 分别测得转子系统有故障时三个节点i,k和h的响应为uri=xri+jyri、urk=xrk+jyrk和urh=xrh+jyrh；

(3) 求得系统这三个节点剩余振动量的响应为Δui=uri-ui、 Δuk=urk-uk、 Δuh=urh-uh，并展开成谐波分量和的形式；

(4) 求得转子系统的频率响应函数矩阵E(jlω)=[K+jωlC-ω2l2M]-1和Al、Bl的值；

(5) 建立下式并计算δn1n2：

δn1n2=abs

(13)

(6) 取δL1L2=δn1n2，L1和L2即分别为故障转子系统的碰摩位置，根据式(12)计算故障力。

其中各个符号的含义参见上节。

4 数值模拟实验

采用图2所示的模型检验此方法的可行性。

图2 转子系统有限元模型

具体参数为：轴径d=10 mm，长度l=420 mm，两个圆盘直径和厚度分别为d1=d2=10 mm和h=20 mm。将系统划分为24个节点23个轴段。设钢的弹性模量为2.1×1011Pa，支承刚度为7×106N/m，系统的偏心在节点9处，大小为mr=5×10-6N·m。

图3 ω=300 rad/s时的诊断结果

假设碰摩分别发生在节点12和节点15处，接触刚度均为2×105N/m，碰摩间隙分别为e1=4×10-5m和e2=3×10-5m，碰摩类型为局部碰摩，碰摩角度分别为0°≤α1≤45°和0°≤α2≤45°。

设三个测点分别为3，7和22。当系统转速为ω=300 rad/s时，采用本文提出的方法得到的计算结果如图3(a)所示，碰摩位置在节点12和15处。从图中可以看出，计算结果与精确值一致。同时也可以用等高线图表示此计算结果，如图3(b)所示，等高线图更直观地看出碰摩位置在节点12和15处；当系统转速为ω=600 rad/s时，计算结果如图4所示，从图中可以看出，诊断的碰摩位置为节点12和15处，验证了本文提出的诊断方法的正确性和可行性。

图4 ω=600 rad/s时的诊断结果

5 结 论

(1) 以单碰摩故障诊断理论以及谐波分量的特征为基础，推导出诊断转子系统双碰摩故障的方法，用数值仿真验证了方法的有效性；

(2) 通过碰摩故障前后转子三个测点的振动响应数据和转速就可以分别确定双碰摩故障转子系统的两个碰摩位置，方法简单有效；

(3) 不仅限于碰摩故障的诊断，对不平衡故障等问题也具有良好的适用性，即故障力可转化为外加等效力、原系统的质量、刚度、阻尼和陀螺力矩矩阵不变；

(4) 基于谐波平衡理论中的谐波分量特性，适用于故障力为周期力，即激励和响应能够展开成谐波分量和形式的故障诊断。

[1] 闻邦椿，顾家柳，夏松波,等.高等转子动力学[M].北京:机械工业出版社，2000.

[2] Hall L D，Mba D.Diagnosis of continuous rotor stator rubbing in large scaleturbine units using acoustic emissions[J].Ultrasonics，2004，41(9):765-773.

[3] Hall L D，Mba D.Acoustic emissions diagnosis of rotor-stator rubs using the KS statistic[J].Mechanical Systems and Signal Processing，2004，18(4):849-868.

[4] Wang Q，Chu F.Experimental determination of the rubbing location by means of acoustic emission and wavelet transform[J].Journal of Sound and Vibration，2001，248(1):91-103.

[5] Cheng Jun-sheng，Yu De-jie，Tang Jia-shi，et al.Local rub-impact fault diagnosis of the rotor systems based on EMD[J].Mechanism and Machine Theory，2009，44(4):784-791.

[6] 张娅，王维民，姚剑飞.双盘转子系统轴向-径向碰摩非线性动力学特性分析[J].振动与冲击，2012，31(12):141-145.

ZHANG Ya，WANG Wei-min，YAO Jian-fei.Nonlinear dynamic behavior of a double-disk isotropic rotor system with axial and radial rub-impacts[J].Journal of Vibration and Shock，2012，31(12):141-145.

[7] Chu F，Lu W.Determination of the rubbing location in a multi-disk rotor system by means of dynamic stiffness identification[J].Journal of Sound and Vibration，2011，248(2):235-246.

[8] Cong Fei-yun，Chen Jin，Dong Guang-ming，et al.Experimental validation of impact energy model for the rub impact assessment in a rotor system[J].Mechanical Systems and Signal Processing，2011，25(7):2549-2558.

[9] 姚红良，李鹤，李小彭，等.旋转机械局部故障力的模型诊断及瞬时故障力识别[J].机械工程学报，2007，43(1):120-124.

YAO Hong-liang，LI He，LI Xiao-peng，et al.Diagnosis of local fault and identification of transient fault force in rotating machinery[J].Chinese Journal of Mechanical Engineering，2007，43(1):120-124.

[10] 姚红良，韩清凯，李凌轩，等.基于谐波分量的转子系统碰摩故障定量诊断方法[J].机械工程学报，2012，48(5):43-48.

YAO Hong-liang，HAN Qing-kai，LI Ling-xuan，et al.Method for detecting rubbing fault in rotor system based on harmonic components[J].Chinese Journal of Mechanical Engineering，2012，48(5):43-48.

[11] Han Qing-kai，Yao Hong-liang，Wen Bangchun.A model-based identification approach for local fault of rotating machinery[J].Advances in Vibration Engineering，2008，7(4):365-376.