刘杨, 石拓, 姚红良, 闻邦椿
东北大学 机械工程与自动化学院, 辽宁 沈阳 110819
收稿日期：2016-09-21
基金项目：中央高校基本科研业务费专项资金资助项目(N140304005，N160313003);中国博士后科学基金资助项目(2014M551105，2015T80269)。
作者简介：刘杨(1982 -)，男，辽宁沈阳人，东北大学副教授；闻邦椿(1930-)，男，浙江温岭人，东北大学教授，博士生导师，中国科学院院士。

摘要：针对由于横向常开裂纹引起的裂纹故障转子系统, 根据Paris理论, 运用由应力强度因子及应变能密度函数推导得到的裂纹柔度系数矩阵, 建立了裂纹转子系统动力学模型.由于含有横向裂纹故障的转子系统的刚度矩阵会发生相应变化, 因此, 在已知柔度系数矩阵的基础上, 通过对比无故障转子系统与裂纹转子系统的动力学方程, 可以推导出裂纹转子系统的剩余量动力学方程.分析裂纹所在轴段两端的受力情况, 结合剩余量动力学方程, 从振动信号中提取出常开裂纹发生的位置及裂纹深度信息.结果表明, 本文运用的常开裂纹信息提取方法理论上只需要采集两个不同转速下, 三个节点的振动信号进行诊断分析, 简便可行, 并可为转子系统、悬臂梁、桁架结构等常开裂纹故障诊断提供一定的理论依据.
关键词：常开裂纹柔度系数裂纹位置裂纹深度Paris理论
Research on Fault Diagnosis Information of Open Crack for Rotor System
LIU Yang, SHI Tuo, YAO Hong-liang, WEN Bang-chun
School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: LIU Yang, E-mail:lyandky@sina.com
Abstract: According to the rotor system with crack fault caused by transverse open crack, the dynamic model of crack rotor system was established using the crack compliance coefficient matrix that is derived from the stress intensity factors and strain energy density function based on the Paris theory. Since the stiffness matrix of rotor system with transverse crack fault will have corresponding changes, the surplus dynamics equation of cracked rotor system can be deduced by comparing the dynamics equations of the fault-free rotor system with that of the crack rotor system based on the flexibility coefficient matrix. Moreover, the information of open crack's location and depth can be extracted from the vibration signal by analyzing stress condition at both ends of the shaft segment where crack exists and combining with the surplus dynamics equation. The results show that the open crack information extraction method used in this paper only needs to collect the vibration signal of three nodes at two different speeds to make diagnostic analysis. This method is simple and feasible, which may provide certain theoretical basis for the open crack fault diagnosis of the rotor system, cantilever beam and truss structure.
Key Words: open crackflexibility coefficientcrack locationcrack depthParis theory
空分设备等大型机组转子系统出现疲劳裂纹故障的现象较为常见.长时间的运行、过高的温度以及加工误差引起的附加力矩等都可能导致转子系统在运行过程中出现横向裂纹.裂纹的进一步扩展会使转子系统产生突发性破坏, 进而引起转子系统运动失稳, 导致轴系破坏事故的发生.分析裂纹转子系统的振动特性、提取振动信号中的裂纹信息可以有效地判断裂纹发生的位置及深度.裂纹信息提取方法的研究已经成为研究裂纹转子问题的一个重要方向.
多年来, 国内外的学者及科研人员对裂纹转子系统中复杂的动力学现象作了大量相关研究, 取得了一些重要的研究成果.Dimarogonas等^[1]通过对裂纹转子系统的研究, 指出裂纹对转子振动特性的影响主要体现在转子刚度矩阵的变化上, 并且使用裂纹柔度函数来表示这种变化.Gounaris等^[2]基于断裂力学的方法, 提出了欧拉伯努利梁模型的裂纹单元.郑吉兵等^[3]分析了裂纹转子在非线性涡动影响下的动力学行为, 特别是系统响应的分叉与混沌特性, 发现许多周期3解随初值的改变而变为其他周期数解的情形.Darpe等^[4]使用非线性呼吸裂纹模型, 研究裂纹转轴纵向振动、横向振动和扭转振动耦合作用下的响应情况.综上所述, 目前对于裂纹故障转子系统的研究主要集中在裂纹故障模型的细化及裂纹故障转子模型的建立上, 尤其在裂纹故障转子系统动力学特性分析方面是当前研究的一个热点方向, 而在裂纹故障信息(裂纹位置和深度等)提取方面的研究还相对较少.本文利用裂纹柔度矩阵建立了裂纹转子系统动力学模型, 通过比较健康转子系统与裂纹转子系统的动力学方程, 得到了裂纹转子系统的剩余量动力学方程.基于剩余量动力学方程, 总结出了常开裂纹发生位置及裂纹深度信息的判定方法.
1 故障转子系统的动力学模型1.1 裂纹模型假设一段给定刚度特性的轴段单元上有一个深度为a的横向裂纹, 在裂纹轴段单元的两端受到力P₁~P₁₂的作用.如图 1a所示, P₁与P₇为轴向力, P₂, P₃, P₈, P₉为剪切力, P₄, P₅, P₁₀, P₁₁分别为绕x, y轴的弯矩, P₆与P₁₂为扭矩.裂纹轴段单元截面半径为R, 裂纹宽度为2b, α, h, dx的含义如图 1b所示.
图 1(Fig. 1)

图 1 裂纹模型参数Fig.1 Crack model parameters (a)—裂纹轴段；(b)—裂纹断面.

根据Paris理论^[5], 可以得到在i方向上由裂纹引起的附加位移u_i:

(1)

式中:J(α)为应变能密度函数, 具体形式参见文献[6]; P_i为i方向上的相关载荷.
宽度为2b的裂纹柔度系数为

(2)

由柔度系数可以组集成裂纹柔度矩阵, 柔度矩阵的阶数取决于轴段上节点的自由度数, 当仅考虑弯曲振动时, 裂纹柔度矩阵为
裂纹柔度矩阵C_c^e中各相关的柔度系数为

(3)

其中:v为材料剪切系数; E为轴段材料的弹性模量; 其余变量的具体形式见附录.
裂纹对转子系统的影响可以通过裂纹轴段的刚度变化来表示.假设无裂纹转子轴段的柔度矩阵为

(4)

其中:I为轴段的截面惯性矩; E为轴段材料的弹性模量; l为轴段的长度.
含有裂纹的转子轴段整体柔度矩阵为

(5)

相应地, 裂纹转子轴段的刚度矩阵为

(6)

其中, T_c为裂纹轴段刚度转换矩阵.
1.2 裂纹转子系统动力学模型转子系统可分为若干个轴段, 一个轴段有2个节点, 当仅考虑弯曲振动时, 一个节点有4个自由度, 转子系统的广义坐标可以定义为
无裂纹转子系统的质量、刚度、阻尼矩阵、陀螺矩阵和激励向量分别为M, K, C, G, F, 则无裂纹转子系统的运动微分方程为

(7)

若将转子系统分为N段, 假设在第R段存在裂纹, 则此时系统的质量、阻尼、陀螺力矩均不发生变化, 只有刚度矩阵K_c变化为
式中, ΔK是ΔK^e=T_c(C_c^e)^-1T_c^T扩展后的形式.
裂纹转子系统的运动微分方程为

(8)

式(8)减去式(7)可得裂纹转子系统的剩余量动力学方程:

(9)

式中, Δq=q_c-q, 为转子系统的剩余振动量^[7].
2 裂纹位置参数的确定由于裂纹仅在第R段存在, 则ΔK的具体形式为

(10)

式中0_i×j为i×j维零矩阵.
因为式(10)中仅ΔK^e部分存在非零元素, 因此式(9)的右侧可变化为

(11)

式中F_c^e为裂纹轴段两端由裂纹引起的附加力.
常开裂纹的各柔度因子在转动中不发生变化, 根据式(11)可知F_c^e是简谐变化的.可将式(9)转换为复数表达形式:

(12)

设Δq=ΔAexp(jωt), 由式(12)可得

(13)

根据轴段两端的受力条件, F_c^e可以用裂纹轴一侧的力和力矩分量来表示^[8-9], 则式(13)变为

(14)

考虑到式(11)中F_c的具体形式, 矩阵E(jω)中只有第4R-3列到第4R+4列具有实际意义.当已知轴段上节点i, j的x和y方向以及节点k任一方向(x或y)的剩余振动量后, 由式(14)可以得到

(15)

(16)

根据式(14)和式(16)可以推出节点k在x或y方向剩余振动量的表达式(以y方向为例):

(17)

(18)

式中, E_kR(jω)=[E_(4k-3)(4R-3)(jω) … E_(4k-3)(4R+4)(jω)], 则裂纹的位置可以通过式(19)进行判断:

(19)

式中, Δy′_k为仿真测得的节点k的响应.
根据式(19), 遍历R=1, 2, …, N, 求得各轴段的λ_R值, 则λ_R值最小的轴段即是裂纹发生的轴段.
3 裂纹深度信息的提取根据式(3)可知, 裂纹深度a的信息包含在裂纹柔度系数中, 所以, 求得转子系统的裂纹柔度系数是解决裂纹深度问题的关键^[10].
由式(11)可以推导出

(20)

式中:q_c^e为转子系统裂纹所在轴段两端的响应矢量; F_c^e可根据式(15)求得.
由式(20)可得

(21)

整理可得

(22)

式中(T_c)^-1=(T_c^TT_c)^-1T_c^T, 则由式(22)可得

(23)

假设[a₁ a₂ a₃ a₄]^T =(T_c^TT_c)^-1T_c^TF_c^e, [b₁ b₂ b₃ b₄]^T=T_c^Tq_c^e, 则由式(23)可得

(24)

在确定了裂纹轴段的位置R后, q_c^e是已知的, 进而, a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄均是可以求得的.式(24)中4个方程却有5个未知的裂纹柔度系数, 为了使裂纹柔度系数能够求解, 取2个不同的转动频率ω, 通过式(8)求解其对应的响应q_c^e.根据式(3)可知裂纹柔度系数不随转动频率ω的变化而变化, 因而在2个不同频率下可相应地得到8个方程去求解5个未知数, 这是一个不相容方程组, 可以通过最小二乘法来求出裂纹柔度系数c₂₂, c₃₃, c₄₄, c₄₅和c₅₅.
图 2为裂纹柔度系数随无量纲裂纹深度变化的曲线图.从图中可以看出, 每一个裂纹柔度系数与无量纲裂纹深度a均一一对应, 将求得的c₂₂, c₃₃, c₄₄, c₄₅和c₅₅分别代入式(3), 可求得对应的无量纲裂纹深度a₂₂, a₃₃, a₄₄, a₄₅, a₅₅, 可利用取平均值的方式求出无量纲裂纹深度a, 进而可以推算出裂纹深度a.
图 2(Fig. 2)

图 2 裂纹柔度系数Fig.2 Crack flexibility coefficient

4 数值仿真以单跨单盘转子为例, 裂纹转子系统经过离散化的有限元模型如图 3所示, 图中小点表示节点, 数字表示轴段号.盘直径为60 mm, 位于轴段3;左、右两个轴承分别在轴段1和13上, 2个支承轴承简化为弹簧阻尼器.转子系统各轴段单元的参数如表 1所示.
图 3(Fig. 3)

图 3 单跨单盘转子模型Fig.3 Model of single-span and single-disk rotor

表 1(Table 1)

表 1 转子系统各轴段单元参数Table 1 Parameters for each unit shaft segment of rotor system mm 参数 1 2 3 4 5 6 7 8 9 10 11 12 13 长度 15 40 40 40 40 40 20 40 40 40 40 40 15 直径 10 10 10 10 10 10 60 10 10 10 10 10 10

表 1 转子系统各轴段单元参数 Table 1 Parameters for each unit shaft segment of rotor system

转子系统弹性模量为2.1×10¹¹ Pa, 密度为7 850 kg/m³, 阻尼选用比例阻尼, 具体数值为2×10³ N·s/m, 轴承刚度为2×10⁶ N/m, 不平衡量为156×10^-6 kg·m.假设裂纹位置在第6轴段的中央, 预设故障裂纹深度为2.5 mm, 参考附录中的无量纲形式, 可知无量纲裂纹深度为0.5.
常开裂纹的柔度系数不随转角的变化而变化, 且响应中仅有1倍频, 需要利用式(9)即剩余量动力学方程来确定裂纹位置.首先得到仿真数据中节点4和节点11的x(水平)和y(竖直)方向位移, 以及节点5的y方向位移, 根据式(19)可以分别求得各轴段的λ_R值.如图 4所示, λ_R值最小处即为裂纹发生处, 可以判断出裂纹发生在轴段6.
图 4(Fig. 4)

图 4 裂纹位置判定Fig.4 Determination of crack location

在确定了裂纹位置后, 通过求得转速为6 000，6 600 r/min时转子系统裂纹轴段两端的响应, 利用式(24)求得c₂₂, c₃₃, c₄₄, c₄₅和c₅₅, 以及其对应的无量纲裂纹深度a₂₂, a₃₃, a₄₄, a₄₅, a₅₅, 如表 2所示.
表 2(Table 2)

表 2 裂纹柔度系数Table 2 Crack flexibility coefficient c22 c33 c44 c45 c55 0.257 5 0.322 4 0.387 2 1.371 6 3.832 7 a22 a33 a44 a45 a55 0.51 0.56 0.5 0.48 0.49

表 2 裂纹柔度系数 Table 2 Crack flexibility coefficient

通过对a₂₂, a₃₃, a₄₄, a₄₅, a₅₅求平均值可以得出无量纲裂纹深度a=0.508, 与预设的无量纲裂纹深度0.5相比偏差极小.参考附录中的无量纲形式, 又因为轴的半径为5 mm, 可以进一步确定出轴的实际裂纹深度为2.54 mm.
5 稳健性测试为了验证该方法的稳健性, 仍使用上述模型, 假设裂纹位置在第6轴段的中央, 无量纲裂纹深度为0.1, 转速为6 000 r/min(100 Hz).假设由于实际工况下的噪声干扰, 测得的高次谐波响应信号服从正态分布, 其方差为0.05.利用转子系统剩余振动量方程进行裂纹位置诊断, 进行20次仿真, 结果如图 5所示.从图可以看出, 该方法有较好的稳健性.
图 5(Fig. 5)

图 5 稳健性测试Fig.5 Robustness test

6 试验验证本试验应用Bently转子试验台进行, 试验台如图 6所示.转子系统转轴长450 mm, 各轴段单元具体参数如表 1所示.在轴段单元6处预制切口, 利用疲劳机扩展形成裂纹, 裂纹深度为3 mm.位移传感器安装于距联轴器始端距离为135 mm(节点5)、355 mm(节点11)处, 运用B & K3560B采集仪, 对振动信号进行采集、滤波、放大、A/D转换, 采样频率3.2 kHz.
图 6(Fig. 6)

图 6 Bently转子试验台Fig.6 Bently rotor test rig

分别采集5 000，5 500 r/min时转子系统节点5和节点11的x和y方向位移, 以及节点6的y方向位移.继而, 仍以剩余振动量方程作为提取裂纹位置的方法, 将试验数据代入式(19), 通过计算可以分别求得各轴段的λ_R值, 如图 7所示, λ_R值最小为0.472×10^-6 m, 即裂纹位置发生在轴段6处.
图 7(Fig. 7)

图 7 裂纹位置判定Fig.7 Determination of crack location

求得裂纹位于轴段6处后, 测量轴段5两端节点(节点5、节点6)在转速5 000，5 500 r/min时的响应, 利用式(24)求得c₂₂=0.405 9, c₃₃=0.505 0, c₄₄=0.775 5, c₄₅=2.405 6和c₅₅=5.934 6；结合式(3)求得无量纲裂纹深度为0.623, 进一步确定出轴的实际裂纹深度为3.115 mm.试验提取到的裂纹位置与裂纹深度信息和实际情况基本相符.
7 结论1) 基于裂纹柔度矩阵建立了裂纹转子系统动力学模型及剩余量动力学方程.通过分析裂纹轴段两端的受力条件, 并结合剩余量动力学方程, 从转子系统振动信号中提取常开裂纹发生的位置及裂纹深度信息.
2) 应用常开裂纹信息提取理论可以为转子系统、悬臂梁、桁架结构等常开裂纹故障诊断提供一定的理论依据.
3) 此方法需要采集两个不同转速下三个节点的振动信号进行诊断分析, 简便可行.但同时该方法也存在不足, 主要体现在理论仿真研究中使用的是简支刚度, 对故障信号没有干扰, 而实际当中轴承支撑刚度和工作环境较为复杂, 增大了故障信息提取难度, 需要在振动信号处理上多研究.
附录:
参考文献

[1]	Dimarogonas A D, Papadopoulos C A. Vibration of cracked shafts in bending[J].Journal of Sound and Vibration, 1983, 91(4): 583–593.DOI:10.1016/0022-460X(83)90834-9
[2]	Gounaris G, Dimarogonas A D. A finite element of a cracked prismatic beam in structural analysis[J].Computers & Structures, 1988, 28: 309–313.
[3]	郑吉兵, 孟光. 考虑非线性涡动时裂纹转子的分叉与混沌特性[J].振动工程学报, 1997, 10(2): 190–197. ( Zheng Ji-bing, Meng Guang. The nonlinear influences of whirl speed on bifurcation and chaos of a cracked rotor[J].Journal of Vibration Engineering, 1997, 10(2): 190–197.)
[4]	Darpe A K, Gupta K, Chawla A. Coupled bending, longitudinal and torsional vibrations of a cracked rotor[J].Journal of Sound and Vibration, 2004, 269(1/2): 33–60.
[5]	Tada H, Paris P C, Irwin G R. The stress analysis of cracks handbook[M]. Hellertown, Pensylvania: Del Research Corporation, 1973.
[6]	Papadopoulos C A, Dimarogonas A D. Coupled longitudinal and bending vibration of rotating shaft with an open crack[J].Journal of Sound and Vibration, 1987, 117(1): 81–93.DOI:10.1016/0022-460X(87)90437-8
[7]	Bachshmid N, Pennacchi P, Vania A. Identification of multiple faults in rotor systems[J].Journal of Sound and Vibration, 2002, 254(2): 327–366.DOI:10.1006/jsvi.2001.4116
[8]	Saavedra P N, Cuitino L A. Crack detection and vibration behavior of cracked beams[J].Computers & Structures, 2001, 79(16): 1451–1459.
[9]	谢平, 杜义浩. 基于信息熵的裂纹转子动力特征分析与诊断方法[J].机械工程学报, 2009, 45(1): 195–199. ( Xie Ping, Du Yi-hao. Crack rotor dynamic feature analysis and diagnosis method based on information entropy[J].Journal of Mechanical Engineering, 2009, 45(1): 195–199.)
[10]	Chondros T G, Dimarogonas A D, Yao J. Vibration of a beam with a breathing crack[J].Journal of Sound and Vibration, 2001, 239(1): 57–67.DOI:10.1006/jsvi.2000.3156