任意斜裂纹转子的耦合振动研究1)

焦卫东 蒋永华 李 刚 蔡建程

(浙江师范大学工学院,浙江金华 321004)

引言

裂纹是旋转机械的常见故障,对设备运行安全的潜在危害很大.转子在旋转过程中,裂纹面承受拉、压应力的交替作用,裂纹开/合(或呼吸)行为弓起转轴刚度的周期性变化,导致转子复杂的耦合振动[1-6].对裂纹转子的耦合振动机理与特性进行研究,是转子裂纹诊断的基础.

对于最简单的横裂纹,即裂纹面同时垂直于转轴轴线与基面,Darpe 等[7]深入研究了裂纹转子刚度变化的机理,分析了裂纹转子的纵向、弯曲与扭转耦合振动的特征.随后,又将横裂纹推广至更一般的横-斜裂纹情形,即裂纹面垂直于基面但不垂直于转轴轴线.从刚度变化机理与耦合振动特性两方面,与横裂纹进行了对比研究.相比于横裂纹,横-斜裂纹导致转子更多刚度参数发生耦合,弓起弯曲、扭转甚至纵向耦合振动[8].纵观国内、外现有的研究,主要集中于横裂纹或横-斜裂纹转子的振动机理与特性问题[9-16].但在某些特殊工况下,转子裂纹会呈现出更加复杂的几何形态,例如裂纹面既不垂直于基面也不垂直于转轴轴线,即任意斜裂纹.例如,在大扭矩和强弯矩载荷作用下裂纹会沿着螺旋方向扩展,从而形成螺旋裂纹或斜裂纹[17];由齿轮啮合力所导致的大扭矩,弓发裂纹的斜向扩展[18].转子旋转过程中,张开型裂纹承受恒定方向的拉应力作用,导致转子刚度发生定值削弱.不同于张开型裂纹,呼吸型裂纹承受拉、压应力的交替作用,激起转轴刚度的周期性变化[1,19].因此当转子包含有呼吸型任意斜裂纹时,由于裂纹面两个方向角的交互作用,其刚度参数的变化规律及其交叉耦合机理明显不相同于横裂纹与横-斜裂纹转子,不同方向上的刚度参数发生广泛、强烈且复杂的耦合,弓起具有不同特征的多自由度耦合振动[20].目前,这方面的研究还不够深入.

本工作采用计算机仿真方法,基于力学以及数值分析理论,研究包含不同类型裂纹特别是任意斜裂纹转子的变刚度机理,以及由此弓发的不同故障激励作用下的耦合振动特性,以期为转子裂纹诊断提供参考数据.

1 裂纹转子的耦合振动机理

1.1 裂纹转子的刚度[7-8,20]

横裂纹、横-斜裂纹以及任意斜裂纹的坐标系统及其空间变换关系如图1 所示.

图1 中,xoy面为基面.图1(d)描述了各坐标系统的空间变换关系.显然,从最基础的横裂纹(坐标系xyz)出发,其裂纹面绕z轴旋转一定角度θ1(0°＜θ1＜90°),即可得到横-斜裂纹(坐标系x′y′z′);横-斜裂纹的裂纹面绕其坐标轴y′再旋转一定角度θ2(0°＜θ2＜90°),即可得到任意斜裂纹(坐标系x′′y′′z′′).

图1 裂纹转子的坐标系统Fig.1 The coordinate systems of cracked rotor

图1 裂纹转子的坐标系统(续)Fig.1 The coordinate systems of cracked rotor(continued)

在材料疲劳裂纹与断裂性能分析中,基于应变能理论的有限元方法得到广泛应用[21-24].例如,文龙飞等[21]重点研究了动载荷作用下扩展裂纹尖端应力强度因子的求解方法;王晓明等[22]将表征能量耗散的变量弓入到应变能函数中,形成新的弹性势的显式表达,从而得到精确匹配实验数据的数值模拟结果.本文采用Timoshenko 梁单元对转子单元进行建模,考虑如图1 所示的纵向、弯曲及扭转所有6 个方向自由度.根据卡斯蒂利亚诺定理,裂纹单元的柔度参数表达为[25-26]

式中,ui和Pi分别为沿着第i个坐标方向的节点位移与节点力.U0为无裂纹单元的应变能,Uc为裂纹导致的外加应变能.

无裂纹单元的弹性应变能表达式为

式中,A=πR2为转子轴横截面积,E为杨氏模量,G为刚性模量,I为转轴截面面积矩,I0为截面极惯矩,αs为Timoshenko 梁剪切系数.

由裂纹导致的外加应变能为

式中,E′=E/(1-v),ms=1+v,v为泊松比.

裂纹面的位移可以用张开、滑移与剪开3 种模式来描述.基于三向应力分析,分别推导3 种位移模式的应力密度因子进而通过面积积分计算Uc,最终得到裂纹单元的各个柔度参数gij,i,j=1,2,···,6.考虑裂纹单元各节点位移qi,i=1,2,···,12 的静平衡条件,有[q1-12]T=T[q1-6]T,T为变换矩阵.从而,裂纹单元的刚度矩阵为Kc=TG-1TT,其中柔度矩阵G=[gi j].

1.2 裂纹转子的振动响应

在全局静态坐标系q下,裂纹转子的运动学方程以矩阵向量形式表达为

式中,Ms,Cs以及Ks分别为总体质量、阻尼以及刚度矩阵,fs则为总体激励力向量.

总体质量矩阵Ms由单元质量矩阵组装而成,单元质量矩阵的计算采用一致性方法[27],考虑六自由度Timoshenko 梁的剪切变形与转动惯量效应.同理,总体刚度矩阵Ks也需要由所有单元(包括无裂纹单元与裂纹单元)的刚度矩阵借助合适的方法组装生成[28].对于总体阻尼矩阵Cs,采用文献[7]或[8]所建议的比例阻尼进行估算,估计式为Cs=αdMs+βdKs.此外,在构建总体激励力向量fs时,需要全面考虑转子有限元模型中各个节点所受到的外部激励作用,例如转子质量不平衡离心力、动-静碰摩力、非线性油膜力或外加的扭转激励等[29-34].

相比于无裂纹转子,裂纹转子的振动响应直接受裂纹单元刚度值的影响,而由裂纹弓起的刚度参数变化则需要借助总体应力密度因子值的符号由振动响应进行估算.实际仿真过程中,假设总体质量与阻尼矩阵Ms与Cs保持不变,只有Ks因裂纹的呼吸行为而不断变化,其值在转子每旋转一度后被更新一次.具体实施时,以节点力估计总体的应力密度因子值,其符号用于确定裂纹闭合线的位置,进而确定裂纹柔度系数积分运算的积分限,获得裂纹单元柔度矩阵与刚度矩阵的估计.再经过静态坐标系变换,即可组装为总体刚度矩阵Ks.Ks连同根据转子最新位置更新的fs一起用于估算下一个旋转角度的转子振动响应.如此不断重复,即可获得裂纹转子的振动响应.

2 仿真分析

考虑一个两端支撑单盘转子系统,裂纹位于转子盘右端靠近盘的位置.转轴长度L=0.7 m,直径D=0.015 m,转子盘质量m=1 kg.整个转轴共划分为14 个单元,包含15 个节点.裂纹几何参数包括裂纹深度a与裂纹面方向角(θ1与θ2),裂纹的类型主要取决于后者.

在随后的仿真计算中,所用到的仿真参数如裂纹单元柔度系数的积分计算参数αi与δi以及纽马克-β 数值算法中的时间步长Δt等,参照相关文献进行选取.

2.1 裂纹单元的变刚度特性分析

裂纹单元长度l=0.7/14=0.05 m,裂纹深度比,裂纹面方位角设置为θ1=45°,θ2=60°.如图2 所示,包含横裂纹(T)、横-斜裂纹(TS)以及任意斜裂纹(AS)单元的刚度参数分别以线型“-----”、“–-–”和“—”表示.图2 中横坐标“CCL Position”意为裂纹闭合线(crack closure line,CCL)位置.CCL 概念由Darpe 等提出,用于精确求解式(3)所示的面积积分[7].

图2 三类不同裂纹的交叉耦合刚度参数(ki j,i, j=1,2,···,6)Fig.2 The cross-coupled stiffness coefficients kij,i, j=1,2,···,6 of three different types of crack

交叉耦合刚度参数kij,i≠j,是造成裂纹转子不同方向振动耦合的内因[35].对于AS 型裂纹转子,刚度参数交叉耦合的现象更明显.如图2,在水平剪切-垂直剪切(k23)、水平剪切-扭转(k24)、水平剪切-垂直弯曲(k25)以及垂直剪切-水平弯曲(k36)等方向均出现强烈的交叉耦合现象.

保持裂纹面方位角θ1=45°不变,θ2在30°到90°之间变化.具有不同方位角的AS 型裂纹单元的刚度特性曲线如图3 所示.图中分别以线型“-----”、“–-–”、“—”、“–+–”以及“–o–”线型按照θ2递增的顺序加以对比描述.

由图3 可见,随着方向角θ2的增大,在纵向(k11)、垂直剪切方向(k33)、扭转方向(k44)以及垂直弯曲方向(k55)的刚度值单调下降.但是,在水平弯曲方向(k66)以及水平剪切方向(k22)则不存在这种变化趋势,刚度参数曲线彼此出现了明显交叉,表明AS 型裂纹的裂纹面方向角θ1与θ2之间的交互作用效应.显然,在整个旋转周期的不同位置,这种交互效应对刚度的影响是不同的,特别是在水平弯曲(k66)和水平剪切(k22)方向,呈现出明显的非线性特征.

图3 方位角(θ2)对包含任意斜裂纹的转子单元的变刚度特性影响Fig.3 The influence of oriented angle(θ2)on stiffness variation of the rotor element including an arbitrary slant crack

接下来,考虑一个包含AS 型裂纹的转子轴段,保持裂纹面方位角θ1=45°,θ2=60°不变,裂纹深度a以1/10 倍的转轴直径D(0.015 m)为增量从0.001 5 m 均匀增大到0.007 5 m,即裂纹深度比(/D)从0.1 均匀增大到0.5,裂纹转子的刚度特性如图4所示.对应于不同裂纹深度的刚度特性曲线,采用与图3 相同的线型按照a或递增的顺序进行描述.

由图4 可以看到,随着裂纹深度a或的增大,AS 型裂纹转子在所有6 个自由度方向上的正刚度参数kii,i=1,2,···,6 均呈现出单调减小的趋势,即裂纹深度越大,裂纹转子的刚度值就越小.而且,这种刚度值单调变化趋势是整体性的、连续性的,发生在转子整个旋转周期的不同转角位置(以裂纹面闭合线位置“CCL position”来确定,见图4 横坐标).虽然裂纹转子的刚度参数值相对于转角位置的变化曲线明显是非线性的,但是对应于不同裂纹深度的转子变刚度特性曲线几乎具有相同的走向或变化方向.

2.2 裂纹转子的耦合振动特征分析

图4 裂纹深度a 对裂纹转子变刚度特性的影响Fig.4 The influence of crack depth a on stiffness variation characteristics of cracked rotor

考虑无外加的扭转激励情况.不平衡质量偏心距为1.6×10-5m,转子转速22 rad/s(3.5 Hz),约等于1/10 倍的弯曲自然频率(35.2 Hz).当无裂纹转子只受到不平衡激励作用时,振动谱图中以旋转频率(或基频)为主,水平与垂直两个方向的弯曲振动基频的幅值水平差不多;当转子包含不同空间方向(θ1=45°且θ2=30°,60°和90°)的任意斜裂纹时,裂纹导致转子的柔度增大,基频分量的幅值水平相比于无裂纹情况也显著增加.在水平与垂直两个方向的弯曲振动谱中,还包括二倍频和三倍频谐波分量.而且,从纵向与扭转方向的振动谱中可以发现明显的基频、二倍频以及微弱的四倍频谐波成分,而此时裂纹转子系统只受到弯曲方向的不平衡激励,并无外加的纵向与扭转方向激励,说明此时系统的振动行为主要受弯曲与扭转耦合机理的支配.对θ1=45°,θ2分别为30°,60°和90°三种情况进行对比研究,发现水平、垂直弯曲方向振动谱的基频、二倍频、三倍频以及扭转方向振动谱的基频、二倍频、四倍频谐波分量幅值水平与裂纹面方向角θ2呈负相关关系,即θ2越大,幅值水平越低;纵向振动谱则相反,其基频、二倍频、四倍频谐波分量幅值水平与θ2呈正相关关系,即θ2越大,幅值水平越高.由于篇幅限制,无外加扭转激励情况的仿真图未给出.文献[8]尽管以TS 型裂纹转子作为研究对象,但在只有不平衡激励的情况下,转子的振动响应特性与这里研究的AS 型裂纹转子具有较大的相似性,可以作为参考,特别是该文献给出的图8(无裂纹转子的不平衡振动响应)和图9(裂纹转子的不平衡振动响应).

考虑存在外加的扭转激励情况.此时,不平衡激励与扭转激励同时作用.采用谐波扭转激励Tsin(ωet),其中T=10 N·m,扭转频率ωe=35 Hz(约等于弯曲自然频率ω0).在不同裂纹面方向角(θ1=45°,θ2=30°,60°,90°)下弯曲(y与z)、纵向(u)与扭转(θ)方向振动响应的时域与频域波形如图5～图7 所示.需要注意的是,当θ2=90°时实际上对应的是横-斜裂纹.

从图5～图7 可以观察到,不平衡激励叠加扭转激励下裂纹转子的变刚度特性明显不同于无外加扭转激励情况.当只有不平衡激励作用时,裂纹转子的刚度参数与转轴旋转频率同步变化且趋势平缓;不平衡激励与扭转激励共同作用时,裂纹转子的刚度参数与扭转激励频率的变化趋势剧烈.

图5 不平衡激励与扭转激励共同作用下任意斜裂纹转子的不平衡振动响应(θ1=45°,θ2=30°)Fig.5 Unbalance response of the arbitrary slant crack rotor with torsional excitation:θ1=45°,θ2=30°

图6 不平衡激励与扭转激励共同作用下任意斜裂纹转子的不平衡振动响应(θ1=45°,θ2=60°)Fig.6 Unbalance response of the arbitrary slant crack rotor with torsional excitation:θ1=45°,θ2=60°

图7 不平衡激励与扭转激励共同作用下任意斜裂纹转子的不平衡振动响应(θ1=45°,θ2=90°)Fig.7 Unbalance response of the arbitrary slant crack rotor with torsional excitation:θ1=45°,θ2=90°

扭转激励作用下裂纹转子刚度的突变特性,是导致拍振与振动调制的根本原因[7].其中,振动调制表现为以扭转激励频率ωe为中心、转轴旋转频率ω 为半带宽的对称边带调制现象.在水平和垂直两个方向上,转轴旋转频率ω 以及扭转激励频率ωe两侧的边带频率的幅值近似相等;弯曲与扭转振动的耦合,使扭转激励频率出现在弯曲振动谱中,而且被旋转频率及其高次谐波分量所调制.扭转激励频率ωe与转轴旋转频率ω 及其高次谐波分量mω 之间的相互作用,导致围绕ωe两侧的边带频率ωe±mω 的出现.在图5～图7 的原始振动响应谱图中,由于各特征频率的幅值彼此差异较大,导致一些围绕扭转激励频率ωe的高级调制边带成分无法被观察到.通过调整原始谱图的纵坐标对其进行局部放大,可以清楚地显示出这些特征频率的存在.例如图5 右上角的局部放大图所示的ωe±mω 成分.在纵向上,与只有不平衡激励情况相比,不平衡激励叠加扭转激励下裂纹转子的振动谱出现明显变化.扭转激励的弓入,导致谱图中高频分量的产生,并在时域波形中出现显著的削波现象,其机理见文献[7],不再赘述.

最后,将不平衡激励与谐波扭转激励联合作用下任意斜裂纹转子的弯曲振动谱特征参数(ω,2ω,ωe与ωe±ω)的幅值对方向角θ2的敏感性综合于图8.横-斜裂纹转子的弯曲振动谱特征参数对方向角θ1的敏感性在图9 中也对比给出.

从图9 可以清楚地看到,横-斜裂纹转子的弯曲振动谱中旋转频率ω 以及扭转激励频率ωe的第一级边带分量ωe±ω 的幅值对方向角θ1相当敏感.随着θ1从10°增大到90°,这些分量的幅值迅速减小;而对于其他的谱分量如旋转频率ω 的二次谐波2ω、扭转激励频率ωe及其第一级边带分量ωe±ω,情况则有所不同.例如,在水平方向,随着θ1增大,旋转频率ω 的二次谐波2ω 的幅值增大,扭转激励频率ωe的幅值减小;而在垂直方向,两者均增大.当θ1=90°时,情况出现了变化,即扭转激励频率ωe的幅值在水平方向上增大而在垂直方向上减小.实际上,此时的裂纹类型已经由横-斜裂纹转变为横裂纹了.

对于如图8 所示的任意斜裂纹转子振动,这些特征频率的幅值与方向角θ1与θ2两者均有着密切的关系,呈现出不同的变化模式.当θ1很小时,例如θ1=30°的情况,随着θ2从10°增大到90°,弯曲振动谱中旋转频率ω、扭转激励频率ωe的第一级边带频率ωe±ω 的幅值均单调减小;当θ1中等大小时,例如θ1=45°的情况,这些特征频率幅值先增大然后减小,趋势变化拐点近似在θ2=36°左右;当θ1较大时,例如θ1=60°或75°的情况,这些特征频率幅值先减小后增大,然后又减小,变化拐点分别近似在θ2=[25.71°,45°]和θ2=[30°,45°]左右.

图8 任意斜裂纹转子振动谱特征参数对方向角的敏感性:θ1=[30°,45°,60°,75°],θ2=10°～90°Fig.8 Sensitivity of spectral parameters of the arbitrary slant crack to orientation angle:θ1=[30°,45°,60°,75°],θ2=10°～90°

图9 横-斜裂纹转子振动谱特征参数对方向角的敏感性:θ1=10°～90°,θ2=90°Fig.9 Sensitivity of spectral parameters of the transverse slant crack to orientation angle:θ1=10°～90°,θ2=90°

横-斜裂纹与任意斜裂纹转子弯曲振动谱中特征频率幅值的不同变化特性,体现了两类裂纹对裂纹面方向角参数敏感性的差异.应该注意的是,在θ1较大的情况,任意斜裂纹转子弯曲振动谱中ω 和ωe±ω两者的幅值在第一个减小阶段的变化趋势相比于图9 所示的横-斜裂纹情况,突变性更加强烈,这说明任意斜裂纹比横-斜裂纹对裂纹面方向角的敏感度更高.此外,与横-斜裂纹相比,如图8 所示的任意斜裂纹具有明显的方向敏感性,暗示了任意斜裂纹转子所具有的不同振动特征.文献[8]已导出一些用于裂纹类型(横或横-斜)辨识的有效准则.本仿真分析所获得的一些结果,则有助于揭示任意斜裂纹转子的振动特征,便于对振动检测所获取的谱特征甚至裂纹参数进行评估,例如弯曲振动谱中旋转基频与二倍频分量、扭转激励频率及其边带成分以及裂纹面方向角等.

3 结论

由于裂纹面两个方向角的交互作用,任意斜裂纹转子在不同方向上刚度参数的变化规律及其交叉耦合机理明显不同于横裂纹与横-斜裂纹.这种交互效应不仅体现在刚度参数的量值上,还体现在刚度曲线的变化趋势上,具有明显的非线性特征.任意斜裂纹所具有的复杂变刚度特性,导致裂纹转子发生复杂的非线性振动,并在纵向、弯曲与扭转多个自由度方向上发生强烈耦合.无论是在不平衡激励还是在扭转激励的作用下,弯曲振动与扭转振动的幅度都更大,危害也更严重.

任意斜裂纹转子相比于其他两类裂纹转子,其振动谱特征也不同.从弯曲振动谱特征参数(如旋转基频与二倍频、扭转激励频率及其边带成分的幅值)的敏感性分析结果来看,这些参数对裂纹面方向参数相当敏感.从而,通过比较弯曲振动谱中垂直与水平方向的扭转激励频率的幅值大小,可以对裂纹面的方向参数进行估计.此外,在只有不平衡激励情况下,扭转振动的时域与频域波形显示:当裂纹类型从任意斜裂纹转变为横-斜裂纹时,特征频率的幅值呈现单调变化,这一结果对于裂纹类型的辨识也是有帮助的.