汽輪機轉子等效應力的計算方法對比及預測分析

李汪繁， 徐佳敏， 李思琦

(上海發電設備成套設計研究院有限責任公司， 上海 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%，可滿足應力監測的工程實踐需求。