汽轮机转子系统中共振组织的全局拓扑规律分析

徐 璐， 褚衍东， 杨 琼， 李险峰

(1.兰州交通大学 机电工程学院，兰州 730070； 2.兰州交通大学 数理学院，兰州 730070)

《世界能源展望》中提到可再生能源的全球份额在发电领域迅猛增长,现有常规发电厂的目标之一是补偿剩余负荷，即使在将来的能源情景中，汽轮机依旧起着重要作用。因此，获取更多汽轮机转子系统动力学信息，对其机理进一步地认识，最终进行非线性动力学设计，具有重要的现实意义和广泛的工程背景。

研究大型转子系统的稳定性与分岔理论具有重要意义。描述系统动力学问题的微分方程具有非线性特性，在一定条件下必然发生分岔和混沌现象[1-5]。刘小峰等[6]利用分岔图、Poincaré截面图等，揭示了交叉刚度对碰摩转子系统运行状态的影响。徐文标等[7]基于Paris裂纹扩展规律，在参数平面内分析了裂纹扩展对含拉杆裂纹的转子系统动力学特性的影响，发现了多种分岔现象。张金剑等[8]利用单参数数值工具探讨了质量偏心、励磁电流和气隙动静偏心对水电机组碰摩转子-轴承系统的影响。我国对于汽流激振的研究起步较晚，是从1980年开始的。柴山等[9-11]对汽轮机的汽流激振力做了较详细的研究，包括叶片优化设计、汽流激振力计算公式的推导及计算程序的设计等。接着，文献[12-17]充分利用分岔图、Poincaré截面图、轴心轨迹图、频谱图等单参数数值工具讨论了高压或低压汽轮机转子系统的分岔和混沌行为。最近，曹丽华等[18-19]以耦合了非线性汽流激振力和油膜力的超超临界汽轮机转子系统为研究对象，基于分岔图、最大Lyapunov指数图等，揭示了非线性汽流激振力对系统的影响。翁雷等[20]对汽轮机转子系统做了较全面的分析，分别探讨了叶尖汽流激振对汽轮机转子、基础松动转子[21]、裂纹转子[22]、碰摩转子[23]及裂纹-碰摩[24]转子系统的影响。罗跃刚等[25]分析了汽流激振和偏心量对密封-转子系统的影响。

综上所述，尽管关于转子系统非线性动力学行为的研究有丰硕的成果，但是都是通过单参数分岔分析实现，导致动力学信息有限且不完整。这表明我们对大型旋转机械的双参数结构的理解还很远，更不用说更高维的参数空间。因此，针对大型旋转机械系统的非线性动力学行为的研究，真正的挑战是发现并总结二维参数平面内分岔行为的拓扑普遍性。

基于以上原因，本文以油膜支撑，在汽流激振力、非线性油膜力和不平衡离心力共同作用下的汽轮机转子系统为研究对象，利用高清的双参数平面周期稳定相图、分岔图、Lyapunov指数图和时间序列图，分析了转子系统内周期吸引子的全局拓扑规律，目的是补充和完善已得到的研究成果，更全面地掌握汽轮机转子系统的非线性动力学特性，揭示更多隐藏在系统中的动力学信息，为大型旋转机械实际设计及故障预测提供有价值的依据。

1 系统力学模型及运动方程

1.1 非线性间隙汽流激振力计算公式

图1为系统受力简图，其中fax，fay是汽轮机非线性汽流激振力Fa在x和y方向的分力。图2为沿环向的两个横剖面展开图。

图1 系统受力图Fig.1 Loading diagram of the system

图2 汽流在静、动叶片间流动示意图Fig.2 The schematic diagram of air-exciting force moving between stator and rotor blades

则非线性汽流激振力的无量纲模型为：

Fa=A1δE+A3δ3E3

(1)

其中

1.2 油膜力计算公式

忽略扭转振动和陀螺力矩，在图3中：O1,O2分别为轴承内瓦和转子几何中心；O3为转子质心；转子两端由半径为R，长度为L的滑动轴承支撑；m1,c1和m2,c2分别为转子在轴承处和圆盘处的等效集中质量和结构阻尼；Fx,Fy为非线性油膜分量；h为平均油膜厚度。

图3 油膜支撑转子系统示意图Fig.3 A schematic diagram of the rotor system supported by oil film

图3是油膜支撑转子示意图，利用短轴承模型[26]，则设x，y为轴承位移，则无量纲油膜力在x和y方向的分力为

(2)

其中

(3)

(4)

(5)

(6)

1.3 系统动力学方程

设x1，y1为转子左端轴承处的径向位移；x2，y2为转盘处的径向位移，则系统的无量纲运动方程为

(7)

式中，τ=ωt。式(7)中的参数意义及取值，如表1所示。

表1 汽轮机转子系统参数意义及取值Tab.1 The parameter meanings and values of the system

2 计算方法

在本文，双参数平面内的数值结果主要使用isoperiodic图[27]呈现。Isoperiodic图不仅显示了一个周期内尖峰数量的积累，即波形的复杂性，以及在参数平面内尖峰数量积累的地方以及积累的快慢，而且也可以有效地用于处理试验数据。对于大型旋转机械系统，制作isoperiodic图是一项非常耗时的任务。因此，采用GPU并行计算技术。

为了生成isoperiodic稳定相图，借助GPU并行计算的优势，将感兴趣的参数窗口划分为N×N的等距点网格，通常可达2 000×2 000=4×106。对于每个点，通过对方程进行数值积分来获得时间演化。对方程组式(7)使用标准四阶Runge-Kutta算法。为了进行积分，从任意选择的初始条件开始，从左到右水平扫描参数，通过“跟随吸引子”向右进行。在所有情况下，都会放弃起初2×105步以避免瞬态响应。在对时间步长进行后续80×105积分后，可以获得isoperiodic稳定相图。每一张isoperiodic图会使用19种颜色的调色板显示每个周期的尖峰数量。通过循环19个基本颜色“模19”来绘制超过19个尖峰的振荡，即通过为它们分配一个颜色索引。黑色用于表示“混沌或拟周期”，即缺乏数字可检测的周期性。

同时，本文采用一种新的方法，即根据运动方程初始小的扰动副本(“克隆”)的时间演化，计算模型的李雅普诺夫指数谱，其中雅可比矩阵的复杂计算替换为线性化版本的运动方程的时间演化。“克隆法”李雅普诺夫指数的简要并行计算代码请参考文献[28-29]。

3 具体数值结果分析

在图4中偏心量e和旋转速度ω同时变化，进汽速度v保持200不变，刚度k=2.5×107。通过右侧的颜色棒很容易区分周期窗口的周期数、大小及形状。将图4(a)的周期序列划分为两组，这两组周期结构在图4(b)和图4(c)中进一步放大详述。这些周期区域具有自相似性。

图4 汽轮机转子系统在(ω,e)参数平面内的稳定周期相图Fig.4 The isoperiodic diagrams of the steam turbine rotor system in (ω,e) parameter plane

在图4(b)中，13条周期数不同的“细丝”嵌入在大片的黑色区域内。为了更加深入地理解这组周期序列的演化细节，沿着蓝色直线做单参数分岔图。在图5(a)中，水平扫描e，随着e的增大且旋转速度ω固定在2 000，13个周期窗口依次锁定在非周期区域内：P-3→P-17→P-14→P-11→P-8→P-13→P-18→P-5→P-17→P-12→P-7→P-16→P-9。图5(b)呈现了通过“克隆”法求得的最大Lyapunov指数图，可见堆积在其中的周期区域将拟周期区域瓜分，这是一个周期解与拟周期解频繁转换的过程。拟周期解随参数变化在某些参数域锁相到周期解，因此这些周期区域称为共振组织或锁模结构。

图5(a)中A点处(e=0.060 8×10-3)的状态传递矩阵的特征值为：λ=[-0.760 7±0.131 8i,0.359 6±0.658 1i, 1.012 1, 0.356 8]，主导Floquet乘子|λ|max=1.012 1>1，即|λ|max通过正实数穿过单位圆，因此在A点处发生了鞍结分岔(SN)。在B点处(e=0.070 5×10-3)的状态传递矩阵的特征值为：λ=[-0.542 7±0.086 0i, 0.446 1±0.702 0i, 1.052 7, 0.211 8]，主导Floquet乘子|λ|max=1.052 7>1，即|λ|max通过正实数穿过单位圆，则在B点处发生了鞍结分岔(SN)。

由此可得：周期解通过鞍结分岔进入拟周期解，即两根“细丝”之间的边界为鞍结分岔曲线。

旋转数不断变化，结果形成了“魔鬼楼梯”[32](如图5(c))。随着e的不断增加，r/s不是线性增长，而是一系列的跳跃和不同大小的台阶，不同的台阶代表着不同的锁模频率，锁模的大小则由台阶的宽度来反映。同时，这些台阶具有分形的自相似性。沿着台阶，整个“魔鬼楼梯”平稳地单调递增，表明周期运动被锁定在关于e的一个不断增加且有限范围内的某个方向上。

图5 沿着图4(b)中蓝色直线的分岔图、指数图及“魔鬼楼梯”Fig.5 The bifurcation diagram,Lyapunov exponent diagram and “Devil’s staircase” along the blue line of Fig.4(b)

图4(c)放大了第二组共振组织。正如图6所示，有7个周期窗口嵌入在拟周期区域内：P-11→P-18→P-7→P-17→P-10→P-13→P-16。通过计算在C点处(ω=1 812)的状态传递矩阵的特征值：λ=[0.530 3±0.721 0i,-0.382 4±0.508 4i,0.130 5±0.495 2i]，由3对共轭复数组成，主导Floquet乘子|λ|max=1.002 5>1，即|λ|max通过其中一对共轭复数穿过单位圆，因此在C点处发生了Hopf分岔。由于Hopf分岔的存在导致了共振组织的出现。

图6 沿着图4(c)中白色直线的分岔图及“魔鬼楼梯”Fig.6 The bifurcation diagram and “Devil’s staircase” along the white line of Fig4.(c)

这组共振组织的旋转数的分布也满足法里树序列的拓扑规律，且形成了单调递增的“魔鬼楼梯”。

图8(c)～图(i)显示了图7(a)中被淡蓝色的A,B,...,G点所标记的参数组合下的强度振荡的时间演化过程，并通过波形变形说明了新尖峰的产生，揭示了一种按有规律的方式所组织的周期振荡。绿色垂直线标记振荡周期T1。图8(c)～图8(i)中的时间演化表明了在转子模式的稳定复杂性中隐藏的规律：每一族似乎都是一些准相同模式的固定组合的串联，尖峰越多，波形就越多。

图7 k×v参数平面内稳定振荡和混沌振荡的曲折分布Fig.7 The distribution of periodic and chaotic oscillations in k×v parameter plane

图8 沿着图7(a)中黄色直线的分岔图、指数图及时间历程图Fig.8 The bifurcation diagram,Lyapunov exponent diagram and time series diagrams along the yellow line of Fig7.(a)

在图7(b)中，随着v的增大，顺着水平方向扫描图像，有8个周期窗口堆积在黑色的拟周期区域内。通过仅仅调整一个参数而得到的单参数分岔图9(a)，相应地，图9(b)中有8个负指数被锁在虚直线以下及图9(c)中单调递减的“魔鬼楼梯”。

图9 沿着图7(b)中蓝色直线的分岔图、指数图及“魔鬼楼梯”Fig.9 The bifurcation diagram,Lyapunov exponent diagram and “Devil’s staircase” along the blue line of Fig7.(b)

图10呈现了参数空间k×ω，k×e，v×e内共振组织的全局拓扑规律。很容易识别出这些不同的双参数平面具有相同的拓扑规律，即法里树序列。因此，满足法里和的共振组织分布是该系统的共性。在旋转机械中，转子的质量偏心会产生振动频率为ω的周期性激振力。若转子弹性力有非线性就会出现许多非线性共振现象。

图10 k×ω，k×e，v×e平面内共振组织的全局拓扑机制Fig.10 The global topological mechanism of periodic oscillations of k×ω，k×e，v×e parameter planes

4 结 论

针对大型旋转机械系统，对其中的动力学现象做详细的分类是一个费时费力的工作，由于GPU并行计算方法的卓越性能和算法的不断改进，本文从二维参数角度揭示了汽轮机转子系统中共振组织排列的拓扑规律，可以得到：满足法里和序列的共振组织分布是该系统的共性。从非线性动力学角度来看，二维参数平面相图显示了波形的复杂性，以及周期吸引子在控制参数空间中积累的速度和位置。因此，这一共性有利于预测在更宽的参数空间动力学行为的演变；从参数匹配的角度来看，周期振荡的展开会表明系统对某些参数的敏感程度，因此实际试验也可以从以上的图表中受益，并提供需要什么样的精度才能区分相图中呈现的周期和混沌窗口的复杂交替。从实际应用角度来看，拓扑规律的提取是将复杂的机械现象转化为简单的数学问题，为优化设计、参数匹配及故障预测提供了可视化、可比较的图表，有利于对转子系统非线性动力学特性更进一步的理解。