基于DYNAMICS软件的透平膨胀机临界转速计算与分析

2018-07-13斌赖欢吕波侯

风机技术 2018年3期

马斌赖欢吕波侯予

（1.中国空气动力学研究与发展中心；2.西安交通大学）

0　引言

随着工业技术的发展，要求旋转机械的工作转速和效率不断提高。转子-支承系统的稳定性是设计转子系统的关键。当系统工作转速接近系统的横向振动频率时，转子会产生强烈的振动和噪声，严重时会造成系统损毁[1]。因此进行转子动力学分析成为系统设计阶段的重要内容之一。现代的转子动力学计算方法可以分为两大类：传递矩阵法和有限元法。本文针对基于传递矩阵算法的商业软件DYNAMICS进行了介绍，对两种不同结构的透平膨胀机转子进行了分析，与通用有限元软件ANSYS计算结果和实验结果进行了对比研究。

1　转子动力学分析的基本理论

一个机械系统的运动微分方程可表达为

式中：M，C，K为系统的质量、阻尼和刚度矩阵；Z为系统的广义坐标矢量；F是作用在系统上的广义外力。对于转子系统，还需要考虑转动带来的陀螺效应、转动刚度和阻尼。因此转子系统的运动微分方程可表达为

式中：C是阻尼矩阵，非对称阵；G是陀螺矩阵，反对称阵；K是刚度矩阵的对称部分；S是它的不对称部分。上述矩阵均是与转速ω有关的函数。

传递矩阵法和有限元法均基于上述理论，求解式（2）的特征值来求解转子的临界转速。传递矩阵法是对各个轴段的截面状态向量进行逐段推算，直到计算到转轴的另一端，主要特点是：矩阵的阶数不随系统的自由度数增大而增加，因而编程简单，占内存少，运算速度快[2]。有限元法是将连续系统离散化以变分原理为基础形成的一种数值近似解法[3]，随着转子有限元模型的不断完善以及计算机的发展，使得有限元可以对大型复杂结构进行数值计算[4]。但系统复杂时会导致自由度数特别大，耗费计算时间[5]。

2　DYNAMICS临界转速计算

DYNAMICS由俄罗斯顶级透平机械工业领袖——Alfa-Tranzit公司研发，其核心技术团队来自国际久负盛名的莫斯科航空航天学院（AMI），合作单位涉及俄罗斯SU系列及米格系列战机引擎的制造商。DYNAMICS是一款针对于旋转机械设计与分析的转子动力学软件，致力于高性能旋转机械转子系统的振动、平衡和稳定科学研究，可以实现转子系统的结构建模、仿真分析和故障诊断。它的分析设计对象从简单的Jeffcott转子到多轴系转子、齿轮系统，再到复杂的机匣-多转子系统，能够模拟简单的刚性支承、正交弹性支承和复杂的滚动轴承、滑动轴承、挤压油膜阻尼器（SFD）支撑等，同时能够提供长短轴承理论、半膜假设、全膜假设等理论，此外还可以对外界阻尼、材料内阻、间隙、碰撞等因素加以考虑[6]。图1为DYNAMICS界面。

图1　DYNAMICS界面图Fig.1　Interface of DYNAMICS

2.1　前处理

DYNAMICS的数据结构由子系统、转配体和子模型三级组成，具有多级转子结构完整模型的文本与图形并存的前处理能力；提供标准的建模单元库，包括梁、壳、集中质量、耦合、支撑、齿轮、超单元等；提供标准类型的材料库；可以施加的载荷包括：不平衡力、一般形式的动态力、叶片飞出、基础运动等；支持转子模型的2D和3D可视化。

2.2　求解计算

软件使用传递矩阵法进行转子动力学求解，计算速度快、精度能够满足工程需求。可以选择的求解模块包括：自振频率、临界转速、不平衡响应、稳定性、瞬态响应等。

2.3　后处理

软件可以输出位移、速率、加速度、反作用力和力矩、能量等计算结果，允许对时间信号进行充分的分析，可以得到平均值、轨迹、傅立叶变换、瀑布图等计算结果。

3　计算结果

本文以箔片动压气体轴承支撑的直径Φ1=25mm双盘转子透平膨胀机和静压气体轴承支撑的Φ2=17mm单盘转子透平膨胀机为研究对象，转子参数见表1，转子结构如图2所示。动压轴承采用铍青铜（QBe1.7）和不锈钢（0Cr18Ni9）两种材质的双层鼓泡箔片动压轴承，轴承长度L1=36mm，轴承名义直径D1=25.04mm，设计半径间隙C1=20μm，刚度与阻尼根据文献[7]计算结果选取，如图3所示。静压轴承采用双排环面切向小孔供气形式，每排8个小孔均布，轴承长度L2=27mm，轴承直径D2=17.04mm，设计半径间隙C2=20μm，节流孔直径d=0.3mm，节流孔距轴承端面距离l=6.5mm，供气压力0.45MPa，刚度和阻尼如图3所示。从图3可以看出，切向供气具有较小的交叉刚度和较大的主阻尼，因此具有较好的稳定性[8]。

表1　转子参数表Tab.1　Parameters of different rotors

图2　转子结构图（单位：mm）Fig.2　structure of different rotors（unit:mm）

图3　轴承刚度和阻尼图Fig.3　Variations of the stiffness coefficients and the damping coefficients versus rotational speed

图4为DYNAMICS分析获得不锈钢双层鼓泡箔片轴承支撑的Φ1=25mm转子系统Campbell图。可见，该系统呈现出高频涡动频率与低频涡动频率两极分化的现象。低频涡动共4阶，转速均在1.5×104rpm（250Hz）以下，而高频涡动（本文只给出了其中的前2阶）的转速都在18×104rpm（3000Hz）以上。图5为Φ1=25mm转子的振型图，8 499rpm，10 182rpm，186 466rpm对应反向涡动，11 698rpm，14 131rpm，208 858rpm对应正向涡动。前4阶为刚体运动模态，5，6两阶为弯曲模态。造成以上结果的原因是转子刚度远大于轴承刚度，使转子弯曲模态的频率远高于转子刚体运动模态的频率。低频的前4阶临界转速与高频的5，6阶临界转速之间较宽的区域意味着转子系统具有较宽的安全工作转速范围。可见1.5×104rpm到18×104rpm之间都是该转子系统的安全工作转速。18×104rpm这一转速已经远远高于透平实际工作转速，因此在理论上透平只需要顺利跨越前四阶临界转速便能在较宽的转速范围内稳定工作。

图4　不锈钢气浮轴承支撑的Φ1=25mm转子系统Campbell图Fig.4　Campbell diagram of the rotor system supported by stainless steel foil bearings(Φ1=25mm)

图5　不锈钢轴承支撑的Φ1=25mm转子振型图Fig.5　Mode sharp of the rotor system supported by stainless steel foil bearings(Φ1=25mm)

图6为DYNAMICS分析获得静压轴承支撑的Φ2=17mm转子系统Campbell图。该系统同样呈现出高频涡动频率与低频涡动频率两极分化的现象。低频涡动共4阶，转速均在2.5×104rpm（416.7Hz）以下，而高频涡动（本文给出了其中的前2阶）的转速都在15×104rpm（2 500Hz）以上。图7为Φ2=17mm转子的振型图，其中17 343rpm，23 482rpm，156 145rpm对应反向涡动，18 053rpm，24 786rpm，173 205rpm对应正向涡动。前4阶为刚体运动模态，5，6两阶为弯曲模态。造成以上结果的原因同样是转子刚度大于轴承刚度，使转子弯曲模态的频率高于转子刚体运动模态的频率。与Φ1=25mm转子系统相比，静压轴承刚度和阻尼大于动压轴承，而其安全工作转速在2.5×104rpm到15×104rpm之间，明显小于Φ1=25mm转子系统。这是由于Φ2=17mm转子为单止推盘偏置，属于细长结构，刚度较低。文献[9、10]分析结果表明相同轴径（Φ=25mm）条件下，这种转子低阶临界转速对轴承刚度和阻尼敏感度较大，稳定性较差，弯曲模态的频率较低，工作范围较小。

图6　静压轴承支撑Φ2=17mm转子系统Campbell图Fig.6　Campbell diagram of the rotor system supported by hydrostatics gas bearings(Φ2=17mm)

图7　Φ2=17mm转子振型图Fig.7　Mode sharp of the rotor system supported by hydrostatics gas bearings(Φ2=27mm)

4　理论计算结果对比

表2和表3分别给出了不同轴承支撑的Φ1=25mm转子前6阶临界转速ANSYS与DYNAMICS的计算结果，可以看出前4阶临界转速的计算结果差异较大，临界转速提高差异变小，第5、6阶临界转速的计算结果差异小于1%。ANSYS计算时间大于100s，DYNAMICS计算时间小于1s。造成低阶临界转速计算结果差异较大的原因是由于ANSYS与DYNAMICS的插值方式不同，当转速小于40 000rpm时，ANSYS采用40 000rpm时的刚度和阻尼，当转速大于75 000rpm时，ANSYS采用75 000rpm时的刚度和阻尼；当转速超出给定范围时，DYNAMICS采用内部拟合值。因此低转速时ANSYS所选刚度和阻尼较大，高转速时ANSYS所选刚度和阻尼较小。因为转子刚体运动模态频率取决于支撑的刚度，轴承刚度增加，转子的临界转速提高[11,12]，而转子弯曲模态频率依赖转子刚度，所以出现上述差异。同理，由于不锈钢轴承刚度大于铍青铜轴承刚度，因此不锈钢轴承支撑转子的前4阶临界转速大于相对应铍青铜轴承支撑转子的前4阶临界转速，而两者的第5、6阶临界转速基本相同。

表2　铍青铜轴承支撑Φ1=25mm转子计算结果Tab.2　Calculation result of the rotor system supported by beryllium bronze foil bearings

表3　不锈钢轴承支撑Φ1=25mm转子计算结果Tab.3　Calculation result of the rotor system supported by stainless steel foil bearings

表4给出了Φ2=17mm转子前6阶临界转速ANSYS与DYNAMICS的计算结果，可以看出各阶临界转速差异均小于2%。ANSYS计算时间为70s，DYNAMICS计算时间小于1s。对比表2～4可知，ANSYS与DYNAMICS计算结果，动压轴承支撑转子系统前4阶临界转速差异较大，静压轴承支撑转子系统前4阶临界转速差异较小。从图3（a）可以看出，动压轴承主刚度Kxx，Kyy随转速变化较大，静压轴承主刚度Kxx，Kyy随随转速变化较小，这是因为动压轴承刚度主要依赖转子转速，静压轴承刚度主要依赖于轴承供气压力。如前文所述，ANSYS与DYNAMICS的插值方式不同，低转速时，静压轴承刚度插值差异小，因此静压轴承支撑转子系统低阶临界转速计算结果差异较小。

表4　静压轴承支撑Φ2=17mm转子计算结果Tab.4　Calculation result of the rotor supported by hydrostatics gas bearings

5　理论计算与实验结果对比

文献[7]对于铍青铜轴承和不锈钢轴承支撑的Φ1=25mm转子系统的自然频率预测值分别为9546rpm（159.1Hz）和11 388rpm（189.8Hz），铍青铜轴承支撑转子非同步涡动振动频率实验值分布在150Hz（9 000rpm）～200Hz（12 000rpm）之间，不锈钢轴承支撑转子非同步涡动振动频率实验值分布在190Hz（11 400rpm）～210Hz（12 600rpm）之间。本文铍青铜轴承和不锈钢轴承支撑转子系统3，4阶临界转速计算结果9 394rpm、11 289rpm和11 698rpm、14 131rpm与非同步涡动振动频率实验结果基本吻合。1，2两阶临界转速缺乏实验数据，高阶临界转速已超出转子工作范围，无法验证。

图8为Φ2=17mm转子轴心轨迹变化图，可以看出随着转速提高，转子轴心轨迹变化明显。转子转动方向为逆时针，转速为20 660rpm时，转子涡动方向为顺时针（反进），转速为21 124rpm时，转子轴心轨迹出现扭曲，是因为转子涡动发生转向。其他转速时转子涡动方向为逆时针（正进）。转速超过21 949rpm后，转子轴心轨迹形状保持稳定，表明转子跨过临界转速进入稳定运行。

图9为Φ2=17mm转子振幅、相位随转速变化图，从图9（a）中可以看出转子X，Y方向的振幅和相位从10 000rpm到30 000rpm时变化剧烈，转速超过30 000rpm后振幅和相位保持稳定。从图9（b）中可以看出在转速19 561rpm、22 948rpm时转子X方向振幅有极大值，转速15 263rpm、20 660rpm和21 656rpm时转子Y方向振幅有极大值。X方向振动相位在19561rpm～20660rpm之间跨过90°。Y方向振动相位分别在15 263rpm～17 770rpm、19 561rpm～20 660rpm、21 124～21 656rpm 之间跨过180°。在转速为20660rpm时，X方向振动相位大于Y方向振动的相位，因此出现图7中涡动与转子转动反向的现象。由文献[13-15]可知，发生相位偏转的转速为转子的临界转速，但本文中转子为涡轮驱动，启动阶段转速很难稳定在特定转速，因此无法确定相位偏转速度。综合振幅与相位判断转子的临界转速为：17 000rpm、19 500rpm、20 500rpm、21 500rpm。计算结果分别为 17 343rpm、18 053rpm、23 483rpm、24 786rpm，计算结果略大于实验结果，造成差异的原因是计算时低转速的刚度和阻尼由插值得到，与轴承实际刚度和阻尼有差距，想要得到准确的临界转速需要根据实际的轴承刚度和阻尼计算。出现X，Y方向振动幅值与相位差异的原因一方面是轴承本身为各向异性，另一方面是由于Y方向有重力作用，造成轴承在X，Y方向的受力以及刚度和阻尼差异较大，从而造成X，Y方向的振动特性有较大区别。