考虑杨氏模量随轴向温度分布变化的转子有限元建模方法研究

何 鹏，刘占生，刘镇星

(哈尔滨工业大学 能源学院，哈尔滨 150001)

汽轮机和燃气轮机在工作过程中，受到工作介质及环境的影响，转子上存在一定的温度分布。在温度作用下，转子发生热弯曲，同时转子材料的力学性能有所变化，影响转子的振动特性。

在温度作用下，转子材料发生热膨胀，引起转子热变形。转子周向温度分布不均匀时，将产生热弯曲，影响转子振动状态。转子的热弯曲振动问题在工程中时常出现［1－3］。转子热弯曲成因可以分为周向介质及上下缸温差［4］、动静碰摩产生的温差［5－8］和转子在滑动轴承中涡动引起的Morton效应［9－10］。对于转子热弯曲所产生的振动问题已经有深入研究［11－12］，并且已经形成了有效的分析诊断方法和处理措施［13－15］。

温度不仅能引发转子热变形，还会引起转子材料杨氏模量的改变。不同材料的杨氏模量随温度变化规律不同，文献［16］通过测量不同温度下振动频率的变化，得到了含镍钢材的杨氏模量随温度变化曲线。文献［17］基于固体力学理论，推导了常用合金材料杨氏模量随温度的变化规律。由温度引起材料杨氏模量的变化对转子振动特性也有影响。晏水平等［18］引入剪切模量随温度的变化规律，分析了轴向温度分布对汽轮机转子轴系扭转振动的影响。朱向哲等［19－20］通过引入热应力效应，考虑温度作用下转子热变形引起刚度矩阵变化，分析了温度分布对转子振动特性的影响。刘少权等［21］引入杨氏模量随温度的变化，采用有限元软件分析了温度场对燃气轮机拉杆转子临界转速的影响。虽然对于温度导致材料杨氏模量变化引起的振动问题已经进行了一些研究，目前尚未有考虑温度引起杨氏模量变化的转子有限元模型，以分析在不同轴向温度分布下转子振动特性。

传统的转子动力学有限元模型［22］无法考虑转子轴向温度分布的影响，本文主要研究考虑轴向温度分布的有限元建模。首先拟合了杨氏模量随温度的变化关系，结合有限元方法，推导了考虑单元内温度分布的有限元刚度矩阵，建立了温度分布对转子振动特性的影响模型。通过数值仿真分析了在温度均匀分布、轴向线性分布和轴向二次分布情况下转子振动特性的变化。

1 杨氏模量随温度的变化规律

材料的杨氏模量是材料原子间作用力的体现，反映了材料应力与应变的比例关系，是材料重要的力学性能参数。对于某一材料，杨氏模量只与材料的温度有关。

不同材料的杨氏模量随温度的变化规律不同。文献［18］中，转子材料剪切模量随温度变化规律为G=G0(1－c(t－t0)/100)，其中 c为常数。在文献［19－20］中，杨氏模量随温度的变化规律为E(t)=207.8－0.069t－3.8 ×10－5t2。在文献［21］中，给出了所计算压气机和涡轮转子材料在不同温度下的杨氏模量值。通过曲线拟合，可以得到文献［21］中压气机转子材料杨氏模量随温度的变化规律为E(t)=E0(－3.862×10－7t2－1.199 ×10－4t+1.002)，其中拟合的均方根误差为0.001923;涡轮转子材料杨氏模量随温度变化规律为E(t)=E0(－5.213×10－7t2－8.347×10－5t+1.002)，其中拟合的均方根误差为0.003989。在文献［16］中，作者通过测量不同温度下叉状物的振动频率，得到了不同温度下含镍钢材杨氏模量随温度的变化曲线。通过曲线拟合，可得该材料杨氏模量随温度变化规律为:

可以看出，不同材料杨氏模量随温度的变化规律略有不同，但大多转子材料的杨氏模量随温度的变化规律可以用二次多项式形式表示:

其中，对于确定的材料，γ，β和α为常数。

2 考虑杨氏模量变化的转子模型

根据杨氏模量随温度的变化规律，结合转子轴向温度分布，可以建立转子有限元单元内杨氏模量的变化趋势，得到考虑轴向温度分布的单元刚度矩阵及转子有限元模型。

2.1 单元内杨氏模量的变化规律

本文主要考虑转子轴向温度分布引起杨氏模量的变化对转子振动特性的影响，故假设转子各个截面温度分布均匀，并且不考虑温度分布引起的转子的热变形。

转子有限元单元如图1所示。假设单元内温度轴向线性分布，两节点温度分别为Ti和Ti+1，则单元内任意位置处温度可以表示为:

需要说明这里可以根据实际转子轴向温度分布情况，假设单元内温度具有高阶多项式分布形式。

结合材料杨氏模量随温度的变化规律，将式(3)代入式(1)中，可以得到单元内杨氏模量的轴向变化规律:

其中f(x)与单元温度分布及轴向坐标有关:

2.2 考虑温度分布单元刚度矩阵

对于如图1所示转子有限元单元，假设轴向方向为z轴，则节点i和节点i+1各具有四个自由度［22］，即{x，y，－ θx，θy}T。在 XOZ 平面内，单元的自由度为{xi，θyi，xi+1，θyi+1}T。单元刚度矩阵为［22］:

其中，E为材料杨氏模量，I为单元截面惯性矩，N为单元型函数。

图1 考虑温度分布的单元模型Fig.1 Element with temperature

结合单元内杨氏模量随温度分布的变化规律，将式(4)代入式(5)，可以得到考虑杨氏模量随温度变化的单元刚度矩阵:

根据f(x)的表达式，可以统一写为如下形式:

其中，ci与单元节点温度和单元内的温度分布相关;n与单元内温度分布形式及杨氏模量随温度变化规律相关。根据杨氏模量随温度二次变化规律和单元内温度线性分布的假设，这里n最大取2。

将式(7)代入式(6)中，可得考虑温度分布的转子单元的刚度矩阵:

各元素的通用表达式为:

一般转子横向刚度相等，在YOZ平面内，单元的刚度矩阵与XOZ平面一致。综合两个平面内的振动，转子单元的横向振动自由度为{xi，yi，－ θxi，θyi，xi+1，yi+1，－ θxi+1，θyi+1}T，单元的刚度矩阵为:

其中:

由于本文只考虑杨氏模量随温度变化对刚度矩阵的影响，不考虑单元的质量矩阵、陀螺力矩矩阵受温度的影响，可采用常规方法得到［22］。

2.3 考虑杨氏模量随温度变化的转子动力学方程

根据单元的质量矩阵、刚度矩阵和陀螺力矩矩阵，可以合并形成转子系统的质量矩阵、刚度矩阵和陀螺力矩矩阵，进而可得转子系统的动力学方程:

其中［M］为转子系统的质量矩阵，［C］为包括陀螺力矩的转子阻尼矩阵，［Ki］为考虑转子系统温度分布的刚度矩阵，{Q}为转子节点载荷向量，{u}为转子节点自由度向量。

3 结果分析

为了体现转子轴向温度分布对转子振动特性的影响和便于验证模型，下面采用简单转子模型，分析轴向温度均匀分布、线性分布和二次分布情况下转子的振动特性。转子的结构如图2所示，转子长0.5 m，轴颈0.04 m，两端刚性支承。转子材料密度为7.85×103kg/m3，常温下杨氏模量为205 GPa，杨氏模量随温度的变化规律采用式(1)函数表示。

图2 模型转子结构Fig.2 The structure of the rotor model

3.1 轴向温度均布转子振动特性变化

假设转子温度均匀分布，可以采用本模型分别计算不同温度下转子的临界转速。针对结构简单，刚性支承转子，临界转速可以采用如下解析式计算［23］:

其中:i为临界转速阶数。

从表达式可以看出，解析式将转子当作简支梁处理，没有考虑转子陀螺力矩的影响。为了便于对比分析，在计算温度均布转子临界转速时，也不考虑陀螺力矩影响，结果如表1所示。

表1 模型转子前两阶临界转速随温度变化Tab.1 Variation of model rotor critical speeds with temperature

图3 模型转子前两阶临界转速随温度的变化趋势Fig.3 Variation of critical speeds with the rotor temperature

转子临界转速随温度的变化趋势及与解析解的对比也可以从图3看出。

从表1和图3看出，随着温度的增加，转子材料杨氏模量减小，转子临界转速降低，本文模型计算结果和解析解比较接近。转子在500℃时，转子材料杨氏模量变为常温值的83.39%，本文模型和解析解的前两阶临界转速均降低8.68%，这验证了本文模型的正确性。

3.2 轴向温度线性分布转子振动特性变化

假设转子温度沿着轴向线性增加，其中低温端固定为20℃，分别分析模型转子在不同轴向温差下的临界转速，结果见表2。

在轴向温差作用下，转子临界转速的变化规律如图4所示。从图中可以看出，随着轴向温差的增加，转子的两阶临界转速是以抛物线规律降低。在转子轴向温差为500℃时，转子一阶临界转速降低了3.41%，二阶临界转速降低了3.57%。

表2 模型转子前两阶临界转速随轴向温差变化Tab.2 Variation of model rotor critical speeds with axial temperature difference

图4 转子前两阶临界转速随线性分布轴向温差变化趋势Fig.4 Variation of the first two critical speeds with axial temperature difference of linear distribution

在轴向温差作用下，转子在高温部位的刚度降低较多，低温部位的刚度相对较大，这可能会使得转子的振型发生变化。为此，对比分析了转子在不考虑温度分布和轴向温差为500℃的线性分布两种情况下转子前两阶振型，结果如图5所示。从图中可以看出，在轴向温差500℃的线性温度分布下，转子的振型与不考虑温度情况之间的差别很小，故可以不考虑轴向温度分布对转子振型的影响。

图5 转子前两阶振型的变化Fig.5 Variation of the first two modal shape

3.3 轴向温度二次多项式分布转子振动特性变化

假设转子的温度在轴向以二次多项式分布，如图6所示为不同的二次多项式温度分布，采用有限元方法分别计算各种二次多项式温度分布下转子的临界转速。

图6 不同温差的转子温度轴向二次分布Fig.6 Several temperature difference quadratic distribution of rotor axial temperature

表3 不同温差二次分布下转子临界转速Tab.3 Rotor critical speeds under several temperature difference quadratic distribution of axial temperature

转子轴向温度在不同二次多项式分布下的前两阶临界转速如表3所示，临界转速的变化趋势如图7所示。转子在轴向温度二次分布下，临界转速变化比线性分布大。在最大温差为500℃时，转子的一阶临界转速降低了7.11%，二阶临界转速降低了5.46%。可见转子轴向温度分布不同，对转子临界转速的影响差别较大。

图7 不同二次多项式温度分布下转子临界转速变化趋势Fig.7 Variation of first two critical speeds under different axial temperature quadratic distribution

从上面三种不同形式轴向温度分布情况下转子临界转速的变化可以看出，本文的模型能够充分考虑转子轴向温度分布的不同，体现转子轴向温度分布对转子材料杨氏模量的影响，为转子热弹耦合动力学分析奠定了基础，可以提高转子在实际工作中临界转速的计算精度。

需要说明的是，本文仿真中所采用的转子轴向温度分布都是理想化的，实际机组中转子轴向温度分布是十分复杂的。例如汽轮机在工作过程中，转子的轴向温度变化规律在各级之间是不同的，同时轴向温度分布也受到汽轮机启动过程的影响。为此，在应用本文模型计算实际转子临界转速时，需要首先确定转子的温度分布，得到转子有限元模型各节点温度。此外，本模型基于转子材料杨氏模量随温度的变化规律，在计算中需针对不同材料，采用不同的变化规律。

4 结论

结合转子材料杨氏模量随温度的变化规律，推导了考虑转子温度分布的单元刚度矩阵，建立了考虑转子轴向温度分布的有限元模型。仿真分析了模型转子在轴向温度均布、线性分布和二次多项式分布下，转子临界转速的变化趋势，结果表明模型能够较好考虑转子轴向温度分布对转子临界转速的影响，提高高温工作转子临界转速计算精度。

［1］陆 山，赵 明，任平珍，等.某型发动机转子热弯曲变形及其影响数值分析［J］.航空动力学报，1997，12(3):243－246.LU Shan，ZHAO Ming，REN Ping-zhen，et al.A numerical analysis of thermal bending deformation and its influence on rotor［J］.Journal of Aerospace Power，1997，12(3):243 －246.

［2］庄景菁.转子热弯曲引起汽轮机振动的分析和处理［J］.华东电力，2002，30(7):61－64.ZHUANG Jing-qing.Analysisand treatmentofturbine vibration caused by thermal bending of rotor［J］.East China Electric Power，2002，30(7):61 －64.

［3］吴子斌.汽轮发电机转子突发性振动原因［J］.东北电力技术，2003，24(1):28－29.WU Zi-bin，The cause of burst vibration of the rotor of turbogenerator［J］.Northeastern Electric Power Technology，2003，24(1):28 －29.

［4］戚先萍，胡璧刚，许都纯，等.转子热弯曲瞬态温度场及热弯曲变形试验与分析［J］.燃气涡轮试验与研究，1997，10(2):32 －36.QI Xian-ping，HU Bi-gang，XU Du-chun，et al.Experiment and analysis of the transient temperature distribution and the thermal bending of rotor［J］.Gas Turbine Experiment and Research，1997，10(2):32 －36.

［5］Goldman P，Muszynska A，Bently D E.Thermal bending of the rotor due to rotor-to-stator rub［J］.International Journal，2000，6(2):91 －100.

［6］黄葆华，杨建刚.摩擦热弯曲对转子稳定性的影响［J］.振动工程学报，2001，14(1):23－28.HUANG Bao-hua，YANG Jian-gang.Influence of thermal bow of the stability of rubbing bearing-rotor system［J］.Journal of Vibration Engineering，2001，14(1):23－28.

［7］田永伟，杨建刚.基于特征值分析的摩擦转子热稳定性研究［J］.振动与冲击，2010，29(7):66－68.TIAN Yong-wei，YANG Jian-gang.Thermal stability of a rubbing rotor based on eigenvalue analysis［J］.Journal of Vibration and Shock，2010，29(7):66 －68.

［8］贺 威，袁惠群，朱向哲.瞬态传递矩阵法分析热弯转子系统的热碰摩故障［J］.东北大学学报(自然科学版)，2009(6):861－864.HE Wei， YUAN Hui-qun， ZHU Xiang-zhe. Transient transfer matrix method for thermal rub-impact faults in hot bending rotor system［J］.Journal of Northeastern University(Natural Science)，2009(6):861－864.

［9］Keogh P S，Morton P G.Journal bearing differential heating evaluation with influence on rotor dynamic behaviour［J］.Proceedings oftheRoyalSocietyofLondon.SeriesA:Mathematical and Physical Sciences，1993，441(1913):527－548.

［10］Keogh P S，Morton P G.The dynamic nature of rotor thermal bending due to unsteady lubricant shearing within a bearing［J］.Proceedings of the Royal Society of London.Series A:Mathematical and Physical Sciences，1994，445(1924):273.

［11］任平珍，柴卫东，胡壁刚，等.航空发动机转子热弯曲稳态响应计算方法研究［J］.燃气涡轮试验与研究，1996(3):27－32.REN Ping-zhen， CHAI Wei-dong， HU Bi-gang， et al.Research on the computation method of the steady state response of the thermal bending aeroengine rotor［J］.Gas Turbine Experiment and Research，1996(3):27－32.

［12］任平珍，陆 山，赵 明.转子热弯曲变形及其影响的数值分析方法［J］.机械科学与技术，1997，16(2):279－282.REN Ping-zhen，LU Shan，ZHAO Ming.Numerical analysis method of thermal deflection and its affection on vibration response of rotor［J］.Mechanical Science and Technology，1997，16(2):279 －282.

［13］陆 山，杜生广，赵 明，等.降低航空发动机转子热弯曲及其影响方法研究［J］.航空发动机，1997(3):23－27.LU Shan，DU Sheng-guang，ZHAO Ming，et al.Research on reducing the thermal bending and its influence of aeroengine rotor［J］.Aeroengine，1997(3):23 －27.

［14］张连祥.航空发动机转子热弯曲引发的典型故障分析［J］.振动与冲击(增刊)，2008，27(S):7－9.ZHANG Lian-xiang.Analysis of typical aeroengine vibration fault related to thermal bowed rotor［J］.Journal of Vibration and Shock，2008，27(S):7－9.

［15］张丽娟，杨新华，马呈霞，等.转子热弯曲引起发电机组振动的诊断与处理［J］.大电机技术，2009(4):13－16.ZHANG Li-juan，YANG Xin-hua，MA Cheng-xia，et al.Remedy on vibration of turbogenerator due to thermal-caused bending of rotor［J］.Large Electric Machine and Hydraulic Turbine，2009(4):13－16.

［16］Kimball Jr A L，Lovell D E.Variation of young's modulus with temperature from vibration measurements［J］.Physical Review，1925，26(1):121 －124.

［17］潘留仙，焦善庆，杜小勇.高温下常用合金材料线胀系数，杨氏模量与温度的关系［J］.湖南师范大学自然科学学报，2000，23(2):47－51.PAN Liu-xian，JIAO Shan-qing，DU Xiao-yong.The relation between linear expansion coefficient and yang's modulus of alloy materials ＆ temperature at high temperature［J］.Journal of Natural Science of Hunan Normal University，2000，23(2):47 －51.

［18］晏水平，黄树红.汽轮发电机组转子温度分布对其扭转振动的影响［J］.中国电机工程学报，2000，20(11):10－13.YAN Shui-ping，HUANG Shu-hong.Impact of temperature distribution on the torsional vibration of turbogenerator rotors［J］.Chinese Society for Electrical Engineering，2000，20(11):10－13.

［19］朱向哲，袁惠群，贺 威.稳态热度场对转子系统临界转速的影响［J］.振动与冲击，2007，26(12):113－116.ZHU Xiang-zhe，YUAN Hui-qun，HE Wei.Effect of steady thermal field on critical speeds of a rotor system［J］.Journal of Vibration and Shock，2007，26(12):113－116.

［20］朱向哲，贺 威，袁惠群.稳态温度场对转子系统振动特性的影响［J］.东北大学学报:自然科学版，2008(1):113－116.ZHU Xiang-zhe，HE Wei，YUAN Hui-qun.Effects of steady temperature field on vibrational characteristics of a rotor system［J］.Journal of Northeastern University(Natrual Science)，2008(1):113－116.

［21］刘少权，张艳春，杜兆刚，等.温度场对燃气轮机拉杆转子临界转速的影响［J］.燃气轮机技术，2011，24(2):20－23.LIU Shao-quan，ZHANG Yan-chun，DU Zhao-gang，et al.Prediction of the influence of temperature field on the critical speeds of a rod-fastened rotor［J］.Gas Turbine Technology，2011，24(2):20 －23.

［22］钟一谔，何衍宗，王 正，等.转子动力学［M］.北京:清华大学出版社，1987:176－185.

［23］GentaG. Consistent matricesin rotordynamic［J］.Meccanica，1985，20(3):235 －248.