汽轮机转子等效应力的计算方法对比及预测分析

李汪繁， 徐佳敏， 李思琦

(上海发电设备成套设计研究院有限责任公司， 上海 200240)

我国现阶段能源革命仍需立足以煤为主的基本国情，抓好煤炭的清洁高效利用。在我国提出“双碳”目标后，国务院、发展和改革委员会、国家能源局等部委配套发布了“十四五”节能减排综合工作方案、全国煤电机组改造升级实施方案等指导性文件，明确要科学发挥煤电的兜底保障作用和灵活调节能力，并将灵活性改造作为存量煤电机组的工作重点之一。在此背景下，煤电机组参与电网调峰尤其是深度调峰时，起停及变负荷将会更频繁、更快速，使得进入汽轮机的蒸汽的热力参数和流量迅速变化，转子与汽缸等关键部件产生较大的温度梯度，热应力随之增大[1]。目前，大型煤电机组汽轮机通常采用双层缸设计来减小汽缸内部温度梯度，可在一定程度上减缓热应力加剧[2]。但对于转子而言，径向温度梯度难以有效大幅减小，且转子在承受热应力的同时还受到离心力作用，因此，其应力计算及监测对于汽轮机的安全稳定运行尤为重要。

国外对汽轮机转子应力的分析研究及监测应用可追溯至上世纪六七十年代[3]，尤其在美国Gallatin电厂2号机组发生转子断裂事故后[4]，转子结构安全性问题更是引起了汽轮机设计制造企业及相关研究院所的关注。我国自上世纪八十年代初开始引进300 MW和600 MW等级的煤电机组，引入阿尔斯通、西屋、西门子等企业基于应力或温差裕度开发的转子应力计算监测模型，逐步自主开展了汽轮机转子应力计算方法、传热系数计算方法、在线监测模型等方面的研究[5-8]。目前行业普遍采用基于模型的转子应力在线监测系统，其主要依据是基于非稳态导热方程的应力简化算法、基于现场数据的温度探针法和基于控制理论的惯性环节法等[9]。随着有限元软件及人工智能的发展，不少学者将有限元软件模拟结果作为精确解，并在完成样本数据的清洗筛选或聚类[10-11]、引入应力修正系数[12-13]等处理后，开展基于神经网络、支持向量机等机器学习算法的转子在线应力监测方法研究。

笔者梳理了汽轮机转子应力场和温度场的简化算法，以某300 MW级汽轮发电机组高压转子为例，开展冷态起动工况下应力场和温度场的简化计算和有限元软件模拟，对比分析了最危险部位的有关计算结果，引入修正系数后采用支持向量回归方法进行应力训练和预测，以期为汽轮机转子应力场和温度场的准确计算及实时监测提供技术支撑。

1 应力场简化算法

应力场简化算法是将汽轮机转子简化为无限长圆柱体，基于热弹性理论得到转子关键部位的等效应力。在起停、变负荷等工况下，转子主要承受热应力和离心应力[14]。

1.1 热应力

热应力是由转子在起停、变负荷等工况下存在的温度梯度引起的。根据转子温度分布和热应力计算公式，可得到转子任意半径r处的径向热应力σrth、切向热应力σθth和轴向热应力σzth[14]。

(1)

式中：E为转子材料弹性模量，MPa；β为转子材料线膨胀系数，1/℃；μ为转子材料泊松比；tM为转子体积平均温度，℃；t为转子任意半径和时间的温度，℃；ri、ro分别为转子内径和外径，m。

1.2 离心应力

转子在起停、变负荷等工况下自身旋转会产生离心应力。假设转子外表面和中心不受其他外力作用，在任意半径r处的径向离心应力σrω、切向离心应力σθω和轴向离心应力σzω计算公式[15]为：

(2)

式中：ω为转子旋转角速度，rad/s；ρ为转子材料密度，kg/m3。

1.3 等效应力

转子在热应力和离心应力的共同作用下，根据力的叠加原理和Von Mises公式，通过式(3)可求得转子等效应力σeq、径向应力σr、切向应力σθ和轴向应力σz。

(3)

2 温度场简化算法

转子热应力计算必须先求得转子温度场，温度场求解为非稳态导热问题，可对几何模型、数学方程、边界条件等作相应处理，采用有限差分法和解析法进行计算[10,16]。相比于有限元商用软件，简化算法实时性较好，所以常被应用于在线监测中。

2.1 有限差分法

常用的有限差分法有一维有限差分和二维有限差分。一维有限差分计算简单快速，转子几何形状和边界条件的适应性较好，其应用也更广泛[16]。一维有限差分是将转子简化为无限长圆柱体模型，忽略轴向和周向换热，将三维导热非稳态微分方程简化为如下一维非稳态导热微分方程：

(4)

式中：τ为时间，s；a为转子材料导温系数，a=λ/ρc，m2/s，其中λ为转子材料导热系数，W/(m·K)，c为转子材料比热容，J/(kg·K)。

求解上述方程需建立初始条件和边界条件[14]：汽轮机高压和中压转子冷态起动温度场计算的初始条件为汽轮机冲转时内缸金属测点的温度，汽轮机低压转子冷态起动温度场计算的初始条件为汽轮机冲转时凝汽器真空对应的饱和温度；温态起动、热态起动和极热态起动温度场计算的初始条件为停机一定时间后的转子温度分布；转子内表面的边界条件为绝热；转子外表面的边界条件为第三类边界条件，其中表面传热系数采用经验公式计算，表面温度采用对应部位的蒸汽温度，表达式为：

(5)

式中：α为转子表面传热系数，W/(m2·K)；tf为转子外表面蒸汽温度，℃；tw为转子外表面温度，℃。

表1 一维非稳态导热差分方程和稳定性条件

2.2 解析法

采用解析法求解汽轮机转子温度，也是将转子视为无限长圆柱体，假设其初始温度处于均匀状态，内外表面的边界条件仍为绝热和第三类边界条件，且将转子材料各项物性和转子表面传热系数视为常数，通过拉普拉斯变换求解一维非稳态导热微分方程[18]，得到转子任意半径处温度t分布为：

(6)

(7)

式中：t0为转子初始温度，℃；ηz为第z段蒸汽温度的变化率，K/s；τz为第z段蒸汽温度的时间，s；Bio为转子外半径Biot数，Bio=αro/λ；βk为βkJ1(βk)-BioJ0(βk)=0的第k个正根，取前2项(k=1，2)已满足精度要求；J0、J1为第一类零阶和一阶Bessel函数。

3 有限元算法

有限元算法已被广泛应用于汽轮机转子应力场和温度场的计算中，商用数值模拟软件的计算精度也得到了实验数据验证[10,19]。因此，可将有限元软件模拟结果作为精确解，与简化算法的计算结果进行对比分析。

3.1 数学模型

由于汽轮机转子为轴对称体，有限元软件求解温度场时，需满足二维非稳态导热方程(8)，初始和边界条件与上述有限差分法一致。对求解区域离散化，进行泛函极值求解和变分计算，可得到n阶线性代数方程组(9)，以求解n个节点的温度[20]。

(8)

Kt-P+N∂t/∂τ=0

(9)

式中：K为刚度矩阵；N为变温矩阵；t为温度向量；P为热载荷向量；∂t/∂τ为温度变化率向量。

在求得温度场的基础上，通过平衡方程、几何方程和物理方程，可求得转子应力场[21]。

3.2 有限元模型及边界

以汽轮机转子实际结构尺寸为依据，建立二维轴对称几何模型，并划分网格，建立有限元模型；以当地蒸汽温度、压力等参数为依据，采用经验公式计算转子轮缘、叶轮、汽封等外表面传热系数[14,22]，并施加第三类边界条件，转子端面施加第二类边界条件，轴径处施加第一类边界条件，再计算转子温度场；施加离心力、压力等力载荷边界，并导入转子温度场，得出转子各部位等效应力计算结果。

4 算例分析

考虑到与温态起动、热态起动、极热态起动等工况相比，冷态起动工况下转子所受到的温度冲击更为剧烈，因此以某300 MW级汽轮发电机组高压转子为例，采用上述有限差分法和解析法计算转子冷态起动工况下最危险部位的温度场和应力场，并与有限元模拟结果进行对比分析。

4.1 初始条件和边界条件设定

在采用有限元软件进行数值模拟时，假设冷态起动工况的初始温度为50 ℃，在转子外表面施加以冷态起动曲线为依据计算得到的转子表面传热系数和蒸汽温度，材料物性以转子实时温度为依据设定，计算得到转子冷态起动工况下的最危险部位为第1级动叶右侧根部，其应力分布如图1所示。

图1 冷态起动工况下转子局部应力分布

采用有限差分法和解析法时，将冷态起动工况的初始温度也设定为50 ℃，并假设转子材料物性和表面传热系数不随温度和蒸汽热力参数的变化而变化，取值以计算结果偏于安全为原则。热应力计算式(1)中热应力系数Eβ/(1-μ)随温度的变化而变化，如图2所示，当温度为400 ℃时，热应力系数最大，因此选取400 ℃为转子材料定性温度。基于参考文献[14]中的传热系数计算公式，取6 000 W/(m2·K)作为转子该部位的表面传热系数。

图2 转子材料特性和热应力系数随温度的变化

4.2 温度计算结果及分析

为便于分析，将冷态起动工况按照起动曲线的温升率变化分为4个阶段，第Ⅰ、Ⅳ阶段的温升率为零，第Ⅱ阶段的温升率略大于第Ⅲ阶段的温升率。

采用3种格式的有限差分法、解析法和有限元软件，得到转子最危险部位在冷态起动工况下的外表面及其对应中心温度的计算结果如图3所示。

图3 冷态起动工况下外表面及中心温度计算结果

从图3可得以下结论：

(1) 4种简化算法计算得到的转子外表面温度较为接近(最大相对偏差仅为0.09%)，但均高于有限元软件模拟结果，最大偏差出现在起动第Ⅰ阶段，绝对偏差为37.8 K(相对偏差为11.17%)，这一方面是因为起动第Ⅰ阶段的负荷较低，有限元软件施加的表面传热系数较小，转子外表面换热不明显，另一方面是因为在有限元软件模拟中，该部位除表面换热外，还会向轴向低温部位导热，造成其温度偏低。

(2) 虽然在起动第Ⅰ阶段，4种简化算法计算得到的转子中心温度略有差别，但整体偏差不大(最大绝对偏差为4.96 K，相对偏差为2.73%)，也均大于有限元软件模拟结果，两者差值在起动第Ⅰ和第Ⅳ阶段出现极值，绝对偏差分别为16.1 K和15.5 K(相对偏差分别为5.90%和3.03%)。起动第Ⅰ阶段出现较大偏差是由于转子外表面温度在该时期存在较大偏差，导致有限元软件模拟得到的中心温度偏小；起动第Ⅳ阶段出现较大偏差是由于起动后期转子温度上升，有限元软件模拟中设置的转子材料导温系数小于简化算法设定常数，转子导热较慢，同时其存在轴向低温导热，使得转子中心温度偏低。

4种简化算法和有限元软件得到的转子冷态起动工况下，外表面温度和中心温度的差值如图4所示。由图4可知，由于冷态起动工况下进入汽轮机的蒸汽的初始温度远高于转子初始温度，导致起动第Ⅰ阶段转子外表面与中心的温差较大，但随着转子温度的均匀化而逐渐减小；起动到达第Ⅱ阶段后，蒸汽温度开始上升，转子外表面与中心的温差逐渐开始增大；而在第Ⅲ和第Ⅳ阶段，由于蒸汽温升率逐渐减小为零，转子径向温度趋于均匀，温差随之减小。

图4 冷态起动工况下外表面和中心温度温差计算结果

4.3 应力计算结果及分析

根据温度计算结果，分别采用应力场简化算法和有限元软件得出转子最危险部位的外表面在冷态起动工况下的等效应力，计算结果如图5所示。

图5 冷态起动工况下外表面等效应力计算结果

从图5可得以下结论：

(1) 在冷态起动第Ⅰ阶段，由于蒸汽对转子产生热冲击，导致转子外表面温度明显高于内部温度，产生较大的径向温度梯度和热应力，使转子外表面等效应力出现最大值(有限元软件模拟结果为521.66 MPa)，并随着转子温度的均匀化而逐渐减小。同时起动到达第Ⅱ阶段后，蒸汽温度开始上升，转子外表面等效应力随之增大。在起动第Ⅲ和第Ⅳ阶段蒸汽温度上升速率逐渐减缓，转子外表面等效应力随之减小。

(2) 简化算法得到的转子外表面等效应力基本一致，最大绝对偏差为5.25 MPa，相对偏差为0.88%，出现在起动300 s时。在起动第Ⅰ和第Ⅲ阶段，通过简化算法得到的等效应力与有限元软件模拟结果绝对偏差的极值分别为75.74 MPa和-21.96 MPa，相对偏差分别为12.68%和29.29%，出现在300 s和11 100 s时。这主要是由于简化算法和有限元软件模拟的温度场分布存在一定差异，说明冷态起动工况下转子外表面的等效应力以热应力为主，外表面和中心的温差对等效应力有一定影响。

基于四点隐式有限差分法和有限元软件模拟的结果，选取300 s、6 750 s和11 100 s这3个时刻，探讨转子径向温度梯度与等效应力的关系，如图6所示。由图6可知，300 s时的转子径向温度梯度明显大于6 750 s和11 100 s时的径向温度梯度，导致该时刻等效应力简化计算结果和有限元软件模拟结果均大于其他2个时刻的计算结果；同时，在300 s时通过简化算法得到的转子径向温度梯度大于有限元软件模拟结果，导致该时刻通过简化算法得到的转子外表面等效应力较大。可见转子径向温度梯度对转子外表面等效应力有着重要影响。

图6 有限差分法和有限元软件计算的转子径向温度变化

4.4 应力预测分析

在研究简化算法和有限元软件模拟结果的基础上，得到冷态起动工况下该转子最危险部位的等效应力修正系数和安全因子，如图7所示。其中，修正系数为等效应力有限元精确解与简化算法结果的比值；安全因子为等效应力与当前温度下材料屈服强度的比值。对于汽轮机转子冷态起动工况，安全因子小于2即认为结构强度安全[14]。

图7 冷态起动工况下等效应力修正系数和安全因子

由图7可知，在起动第Ⅰ阶段前期，转子等效应力较大，修正系数在0.88～1.15之间，其中最大等效应力点(有限元解为521.66 MPa)的修正系数为0.88，简化算法偏向安全考量；在起动第Ⅱ和第Ⅲ阶段，修正系数为1.2～1.41，其中最大等效应力点(有限元解为102.23 MPa)的安全因子为0.26，对转子安全性影响较小；在起动第Ⅰ阶段的后期和第Ⅳ阶段，等效应力相对较小，修正系数较大，但最大等效应力点(有限元解为32.88 MPa)的安全因子仅为0.09，不会引起较严重的安全性问题。在工程实践中，可清洗筛选样本数据后再建立等效应力预测模型。

基于汽轮机转子简化算法和有限元软件模拟结果，以主蒸汽压力、主蒸汽温度、金属初温、转速、功率、功率变化率、温度变化率等作为特征值输入，冷态起动工况下该转子最危险部位的应力修正系数为目标值，形成60组样本数据，分别选取45组和15组作为训练集和测试集。采用基于交叉验证的支持向量回归(CV-SVR)方法[23-24]预测，得到的计算结果如图8所示。

由图8可知，测试集验证时平均相对偏差为1.82%，最大相对偏差出现在等效应力为72.17 MPa时，此时的安全因子仅为0.20，对转子安全性影响较小；所有等效应力大于100 MPa的样本点中，最大相对偏差仅为0.54%，可满足在线应力监测的精度要求。

图8 转子外表面等效应力训练和预测结果

5 结 论

(1) 有限元软件模拟得到的转子温度场和应力场，与简化算法所得的计算结果存在一定偏差：外表面和中心温度最大绝对偏差均出现在起动第Ⅰ阶段，分别为37.8 K和16.1 K；外表面等效应力在起动第Ⅰ和第Ⅲ阶段均存在较大偏差，分别为75.74 MPa和-21.96 MPa。在工程实践中，可通过清洗筛选样本数据、引入应力修正系数、设置安全因子监测范围等处理手段优化并建立转子等效应力的预测模型。

(2) 简化算法和有限元软件得出的结果均表明，转子外表面和中心温差对转子外表面等效应力有一定影响，且转子外表面径向温度梯度与转子外表面等效应力呈正相关关系，可作为转子关键部位等效应力衡量的依据。

(3) 采用CV-SVR方法训练并预测汽轮机转子等效应力，等效应力较大时(大于100 MPa)，预测结果的最大相对偏差仅为0.54%，可满足应力监测的工程实践需求。