裂纹角影响下的裂纹转子非线性分析

刘桂珍，于 影，马晓君，李小海，闻邦椿

（1.佳木斯大学机械工程学院，黑龙江佳木斯154007；2.东北大学机械工程与自动化学院，辽宁沈阳110819）

裂纹故障是旋转机械常见的故障之一。裂纹产生的原因主要是由于转子材料本身缺陷或长时间的服役而产生了疲劳裂纹，裂纹故障是一种后果严重、诊断困难的常发性故障［1］。裂纹的存在使转轴刚度周期性变化，其潜在危害性远超一般故障，因此转轴裂纹故障的研究特别引起工程界的重视［2］。自20世纪60年代末期在汽轮机轴上发现裂纹以来，各国学者纷纷开展裂纹转子振动特性的研究工作［3］。裂纹周期性的张开与闭合引起转轴刚度发生相应的变化，是大型旋转机械运行不稳定的重要因素［4］。裂纹角是裂纹模型中的一个重要参数，对裂纹转子的动力学行为有着重要的影响［5］。

本文以带有裂纹故障的转子系统为研究对象，利用拉格朗日方程建立了非稳态油膜力作用下的转子系统裂纹故障的4质量8自由度非线性动力学模型；基于裂纹转子系统的强非线性和对初值的敏感性，应用数值分析研究裂纹角作为唯一控制参数时系统的动态响应；提取了转子系统随裂纹扩展的动态过程，对该系统响应的非线性行为和故障机理进行分析，从而为该类转子系统的故障诊断和系统的安全运行提供理论依据。

1 力学模型与微分方程

1.1 拉格朗日方程描述

设有n个质点组成的质点系，受完整的理想约束，具有N个自由度，其位置可由N个广义坐标方程来确定。则有

式中：L为拉格朗日函数，L=T-V；T为系统的动能函数；V为系统的势能函数；R为与系统阻尼相对应的耗散函数；Qi为作用在系统上的广义力；qi为系统独立的广义坐标；N为系统的总自由度个数。

1.2 刚度影响系数的计算

设转轴半径为R、长度为L1的无质量弹性圆轴，在转轴中央有一深度为a的弓形横向裂纹，如图1所示。

图1 裂纹截面示意图Fig.1 Crack section of the rotor

如果只考虑裂纹处弯矩的作用，根据断裂力学理论，由于裂纹的存在将在裂纹局部产生附加角位移。设在η、ξ方向弯矩作用下的无量纲局部柔度系数分别为［6］

式中：v为泊松比；E为杨氏模量；局部裂纹深度为最大裂纹深度；h=

几何修正因子F2(η h)、F1(η h)分别为

对式（2）进行数值积分，得到无量纲的局部柔度系数与无量纲裂纹深度的关系曲线，如图2所示。

图2 柔度与无量纲裂纹深度关系Fig.2 Relation between flexibility and the crack deepness

ε、δ为仅与裂纹深度a有关的相关刚度参数，其表达式如下：

式中：L1为两端支承间轴的长度；R为轴颈；v为泊松比。

1.3 非稳态油膜力模型［7］

非线性油膜力模型采用短轴承假设下的Ca-pone 非线性油膜力模型，该模型有较好的精度和收敛性。在短轴承油膜力假设条件下的无量纲雷诺方程为

式中：h为无量纲油膜厚度为油膜厚度；C为轴承径向间隙；z为无量纲轴向位移，z=zˉL；zˉ为轴向位移为无量纲油膜压力；pˉ为油膜压力；μ为油膜黏性系数；x、y分为无量纲轴颈中心x、y方向的位移；x˙、y˙为无量纲轴颈中心x和y方向上的速度分量。

由式（5）可得无量纲油膜压力为

式中：D为轴承直径。

无量纲非线性油膜力最终可以表示为

式中：

1.4 运动微分方程

图3 所示为转子-定子-轴承系统中间轴上具有一横向裂纹的力学模型，两端由滑动轴承支撑，滑动轴承直径为D，长度为L。两轴承之间为一无质量弹性轴，其半径为R，长度为L1，转轴中央有一对称布置的圆盘，以及深度为a的弓形横向裂纹。O1为转子的几何中心；Oc为转子质心；mi（i=1，2，3，4）分别为转子、轴承、定子和转轴在轴承处的半集中质量，单位kg；ki（i=1，2，3）分别为转轴、轴承支撑处和定子基础的刚度系数，单位为N·m-1；ci（i=1，2，3）分别为转轴、轴承支撑处和基础对定子的阻尼系数，单位为N·s·m-1。

图3 非稳态油膜力的转子-轴承系统裂纹故障的力学模型Fig.3 Mechanical model of the rotor-bearing system with a transverse crack

设转盘在水平方向、垂直方向移动的位移为x1、y1，左端轴承在水平方向、垂直方向移动的位移为x2、y2，定子在水平方向、垂直方向移动的位移为x3、y3，右端轴承座处水平方向、垂直方向移动的位移为x4、y4。

采用耗散能量形式的拉格朗日方程推导出对应裂纹故障的转子-定子-轴承系统的质点运动微分方程：

式中：xi、yi(i=1，2，3，4)分别为应各质量的位移坐标；b为轴颈的偏心量，单位为mm；Fx、Fy分别为水平、铅垂方向的油膜力，单位为N；ε、δ为仅与裂纹深度a有关的相对刚度参数，如式（3）所示；F(ψ)为裂纹开闭函数，ψ=τ-φ1+β-arctany/x，β为裂纹方向与偏心之间的夹角，φ1为初相位，这里 设φ1=0；σ为 Sommerfeld 修 正 数 ，σ=为无量纲圆盘偏心量；τ为无量纲时间，τ=ωt。

2 数值模拟

运用四阶Runge-Kutta 法对数值进行求解。在计算中为了能够较快地得到稳定解，应将步长选得尽量小且周期足够多。为了消除瞬态响应的影响，舍弃前40 个周期，计算轨迹图时取后10～20 个周期。选取系统参数如下：m1=4.0 kg，m2=32.1 kg，m3=50.0 kg，m4=20.0 kg；c1=1.05 kN ·s ·m-1，c2=2.1 kN·s·m-1，c3=2.1 kN·s·m-1；k1=250 kN·m-1，k2=250 kN·m-1，k3=25 000 kN·m-1，kr=1 000 kN·m-1；R=0.025 m，L=0.570 m，δ2=0.2 mm；η=0.018 MPa；R1=0.015 m；f=0.2；通过计算得出转子系统的三级固有频率分别为f1=13.974 7 Hz，f2=43.586 4 Hz，f3=113.109 6 Hz。

图4（a）为裂纹角β=3/4π 时系统响应随激励角频率ω变化的分岔图。结合此状态下的庞加莱映射图，可以观察到，当激励角频率ω∈（100，120）rad/s时，系统响应为拟周期运动；当ω∈（150，210）rad/s时，系统响应为周期运动；当ω=240 rad/s时，系统响应为拟周期运动；当ω=275 rad/s时，系统响应为P4周期运动；当ω=300 rad/s时，系统响应周期运动；当ω∈（350，450）rad/s 时，系统响应为拟周期运动；当ω∈（500，550）rad/s时，系统响应为混沌运动；当ω=600 rad/s 时，系统响应为 P2 运动；当ω=650 rad/s时，系统响应2倍拟周期运动；当ω=700 rad/s时，系统响应为混沌运动；当ω=750 rad/s 时，系统响应为拟周期运动；当ω=800 rad/s 时，系统响应为混沌运动。

在裂纹角β=3π∕4 的情况下，随着激励频率的小幅变化，转子系统快速经历了“拟周期运动→周期运动→拟周期运动→P4→P1→拟周期运动→混沌运动→P2→2 倍拟周期运动→混沌运动→2 倍拟周期运动→拟周期运动→混沌运动”的复杂过程。

图4（b）是裂纹角β=π时系统响应随激励角频率随ω变化的分岔图。结合此状态下庞加莱映射图，可观察到：当ω=100 rad/s时，系统处于拟周期运动；当ω=120 rad/s 时，系统响应为混沌运动；当ω∈（150，210）rad/s时，系统转为周期运动；当ω∈（240，275）rad/s 时，系统为拟周期运动；当ω=300 rad/s时，系统为周期运动；当ω∈（350，400）rad/s 时，系统响应为拟周期运动；当ω=450 rad/s 时，系统响应为周期运动；当ω∈（500，550）rad/s时，系统处于混沌运动；当ω=600 rad/s 时，转子系统响应为P2 运动；当ω=650 rad/s 时，转子系统响应为2 倍拟周期；当ω∈（700，750）rad/s时，转子系统响应为混沌运动；当ω=800 rad/s 时，转子系统响应为拟周期运动。

图4 不同裂纹角β时系统响应随ω变化的分岔图Fig 4 Bifurcation diagram with ω changes in the different crack angle β

在裂纹角β=π 的情况下，随着激励频率的小幅变化，转子系统响应快速的完成“拟周期运动→混沌运动→P1运动→拟周期运动→P1运动→拟周期运动→P1 运动→混沌运动→P2 运动→2 倍拟周期运动→混沌运动→拟周期运动”的复杂过程。

3 结论

对带有裂纹故障的转子-定子-轴承系统进行非线性动力学研究，得出以下研究结论：

（1）当裂纹转角β参数不变，伴随激励频率的变化，数值仿真结果表明带有裂纹故障的转子系统以“周期运动→拟周期运动→混沌运动”为主线，呈现出丰富的非线性行为。

（2）激励频率比较低时，裂纹角对转子动力学特性影响不大；当转子工作转速较高时，伴随着裂纹转角大小的变化，转子系统呈现出复杂的非线性动力学行为，以激励频率ω=275 rad/s 时为例，裂纹角β=3π∕4 为系统响应为P4 周期运动；裂纹转角β=π时，则系统转为拟周期运动。

（3）当转子系统具有裂纹故障时，随着转轴刚度的变化，系统运动在较短的速度变化范围内重复出现“周期运动→拟周期运动→混沌运动→周期运动→拟周期运动→混沌运动”的交替运动过程，且变化敏感的非定常运动的范围非常小；因此，应对具有裂纹故障的转子系统进行诊断时需多加注意。