基于三維熱彈性接觸的汽輪機調節級葉片樅樹型葉根輪緣結構的優化

張明輝，張 荻，謝永慧

（西安交通大學 能源與動力工程學院，西安710049）

隨著經濟的快速發展，電力需求持續增長，對發電機組的可靠性要求也不斷提高.葉片作為汽輪機的重要零部件，主要功能是將熱能轉化為機械能，其安全可靠性直接影響機組的正常運行［1］.

調節級位于汽輪機通流部分之首，由于大負荷、高溫及部分進汽等因素影響，在大功率汽輪機所有零部件中，調節級的工作條件是最苛刻的［2］.樅樹型葉根具有裝卸方便、葉根輪緣承載面接近等強度、承載能力高等優點，因而被廣泛應用于大功率汽輪機調節級葉片上［3］.目前的樅樹型葉根輪緣形狀多是在常規狀態下設計的，而調節級葉片工作在高溫蒸汽中，由于葉片和轉子的材料特性存在差異，線膨脹系數不同，隨著溫度的提高，葉根和輪緣存在很大的熱變形，導致葉根輪緣接觸不良，葉根齒的承載及應力狀況發生很大的變化.因此，有必要在考慮高溫熱變形的影響下，對調節級樅樹型葉根輪緣結構進行優化，提高葉片的安全性.

近年來，一些學者開展了通過改進葉根結構提高葉片性能的研究［4-14］，獲得了一些成果.如西安交通大學渦輪機教研室［4］采用平面應力有限元分析程序對樅樹型葉根輪緣應力進行了研究，得到了此類葉根輪緣合理結構設計的一些初步結論.李佳其等［5］、Papanikos等［7］和Meguid等［8］分別對叉型葉根、燕尾型葉根和樅樹型葉根結構進行了分析，初步討論了若干設計尺寸對葉根安全性的影響.2002年，Song等［9］以經濟性為目標，對二維樅樹型葉根進行多變量優化設計，獲得了較好的優化效果.2005年，趙海［11］通過UG 軟件對渦輪榫頭／榫槽三齒結構進行了參數化建模，分析了主要設計參數及溫度載荷對結構強度的影響.2009年，林香等［14］利用Ansys軟件進行優化分析，降低了航空發動機燕尾型榫連接件的最大應力和接觸應力.

國內外研究者在葉根輪緣的優化方面已經進行了一些研究，但大多集中于特定尺寸參數對目標函數的影響或對實際模型進行簡化后的二維分析上，針對基于三維有限元分析的樅樹型葉根輪緣結構優化，尤其是考慮高溫影響下的優化研究仍然較少.近年來，隨著加工工藝的提高及優化理論和優化算法的成熟，合理地設計葉根輪緣尺寸、提高葉片的強度性能變得切實可行且非常必要.

基于以上討論，筆者利用Ansys軟件建立了調節級樅樹型葉根輪緣的參數化模型，在考慮高溫影響下，以減小葉根輪緣處最大等效應力為目標，對葉根輪緣的關鍵尺寸進行了優化設計，最終獲得了該葉根輪緣的優化型線.

1 熱應力問題及結構優化方法

本文研究中考慮了高溫熱應力對汽輪機葉片的影響，并采用現代結構優化方法對調節級樅樹型葉根輪緣結構進行了優化.

1.1 彈性熱應力問題的有限元方程

物體溫度發生變化時，由于熱變形將產生線應變.如果物體各部分熱變形不受任何約束，則物體發生變形而不產生應力.當存在約束或者溫度變化不均勻，熱變形不能自由進行時，物體中就會產生應力，即熱應力［15］.當彈性體的溫度場已知時，可以求解出彈性體各部分的熱應力.

物體由于熱膨脹只產生線應變，所以剪應變為零，這種由于熱變形產生的應變可以看作是物體的初應變：

式中：α為材料的線膨脹系數；φ為結構的穩態或瞬態溫度場；φ0為結構的初始溫度場.

當物體中存在初應變時，應力與應變的關系可表示成

式中：D為彈性矩陣.

將式（2）代入虛位移原理的表達式，則可得到包括溫度應變在內求解熱應力問題的最小位能原理表達式，其泛函表達式如下：

式中：Πp為結構總勢能；u為結構位移；Ω 為體積域；Γ為表面域；f為體積力；T為表面力.

在求解域Ω 內進行有限元離散，得到有限元求解方程為：

式中：K為結構 剛度矩陣；a為節點位移向量；P為載荷.

式中：Pf為體積載荷項；PT為表面載荷項；Pε0為溫度應變引起的載荷項.

式中：B為應變矩陣.

由此可見，熱應力問題與無熱載荷的應力分析問題相比，除了增加一項以初應變形式出現的溫度載荷Pε0外，其他完全相同.筆者充分考慮了汽輪機運行溫度場對調節級葉片及輪緣的影響，建立了葉片與輪緣的三維熱彈性接觸有限元模型.

1.2 結構優化方法

現代結構優化理論由有限元計算和優化方法兩部分組成［16］.目前為止，有限元分析方法已經廣泛應用于航空航天、建筑和機械制造等領域，而優化方法也已形成了較完善的體系.

一般而言，基于參數化有限元分析的優化設計包括3個基本要素［17］（1）設計變量：設計過程中需要不斷調整賦值的變量參數；（2）狀態變量：設計過程需要滿足的約束條件變量，是設計的因變量，設計變量的函數；（3）目標函數：設計中極小化的變量參數必須是設計變量的函數.

最優化問題的數學模型是找到一組設計變量x＝［x1，x2，…，xn］T，滿足

式中：x為設計變量，x＝［x1，x2，…，xn］T；f（x）為目標函數；gi（x）和hj（x）為約束條件或狀態變量.

即在滿足約束條件的前提下，找到一組最優的設計變量，使目標函數值達到最小.

結構優化的主要步驟包括建立數學模型和優化迭代控制［18］.建立數學模型即確立包括設計變量、約束條件和目標函數的分析模型；優化迭代主要是綜合考慮優化效果和求解時間，選擇一定的數學優化模型（即優化方法）.

筆者基于有限元分析軟件Ansys建立調節級樅樹型葉根輪緣的數學模型，采用模式搜索算法進行優化控制，完成結構優化工作.

模式搜索算法是一種直接優化方法，不依賴目標函數的導數信息求解最優化問題，在解決不可導函數或求導異常麻煩的函數的優化問題時非常有效［19］.本文使用的模式搜索算法由Ansys的參數化設計語言APDL編寫，并通過對一個多峰函數求極值來驗證其準確性.

對于函數

其中，-1.204≤x1≤1.204，-1.204≤x2≤1.204.理論解為：x1＝x2＝0時，f取得極大值4.3.

采用本文編寫的算法，以x1和x2為設計變量，-1.204≤x1≤1.204，-1.204≤x2≤1.204為狀態變量，-f為目標函數進行優化，獲得的優化結果為x1＝3.815×10-7，x2＝2.289×10-7時，f取得極小值-4.300，所得結果與理論解非常吻合，從而驗證了該算法的準確性.

2 計算模型

本文的研究對象是考慮高溫影響并具有復雜接觸狀態的汽輪機調節級葉片葉根和輪緣.計算模型包括葉片和對應的轉子部分，圖1給出了計算模型的三維造型、整體網格及葉根輪緣接觸處的網格.葉片和輪緣材料特性見表1.

圖1 計算模型及網格Fig.1 The computational model and grid

表1 葉片和輪緣材料特性Tab.1 Material properties of the blade and rim

有限元模型共有166 618單元、162 113節點.計算中對整個有限元模型設置了汽輪機工作轉速3 000r／min的離心力.另外，根據相關研究［20-21］，汽輪機穩定運行時，調節級處沿軸心至轉子外表面的徑向溫度梯度較小，尤其是轉子頂部與葉根直接接觸的輪緣部分的溫度基本不變，為調節級的工作溫度.因此，筆者依據相關透平制造企業經驗，對計算模型設置調節級正常運行下的工作溫度為534 ℃.計算模型的邊界條件設置如下：（1）葉根輪緣承載設置接觸邊界面條件；（2）輪緣對稱面節點設周期對稱邊界條件；（3）輪緣底部一列節點設徑向約束；（4）輪緣底部中間節點設軸向約束.

3 結構優化過程及參數設置

3.1 結構優化過程

葉根輪緣結構優化是通過改變模型的某些關鍵尺寸參數完成的，在建立有限元模型的過程中，設計變量和約束變量必須以參數的形式出現［22］.因此，本文利用Ansys的參數化設計語言APDL 設置設計變量、狀態變量（約束條件）和目標函數，完成葉根輪緣多參數計算模型的建立以及有限元計算的循環進行.

圖2給出了葉根輪緣結構優化的流程圖.

具體過程如下：

（1）為設計變量賦值，判斷是否滿足狀態變量.如不滿足，通過優化控制程序重新為設計變量賦值，如滿足，則創建參數化模型.

（2）加載求解，進行有限元分析.

（3）由有限元分析結果中提取目標函數值，判斷是否收斂.如果收斂，則退出優化循環，得出結果，否則通過優化控制程序重新為設計變量賦值并進行循環分析.

（4）得到最優方案后，對結果進行后處理.

圖2 優化過程流程圖Fig.2 Flow chart of the optimization process

3.2 優化參數設置

（1）設計變量

綜合考慮葉根輪緣各尺寸對應力分布的影響，選取8個關鍵尺寸作為設計變量進行研究，分別為葉根展開角度A、接觸面傾角B、齒面間距Tp、齒面寬度Cl、葉根齒內圓角半徑R1、葉根齒外圓角半徑Rc1、輪槽外圓角半徑R2和輪槽內圓角半徑Rc2.圖3為變量示意圖.

（2）狀態變量

為了保證優化過程中創建的模型在幾何上合理，避免怪異的形狀出現，模型尺寸必須滿足一定的幾何約束，即優化過程中的狀態變量.如圖3所示，對模型設置了5個狀態變量：葉根輪緣非接觸面的傾角C＞0；葉根齒面寬度Cl與輪緣承載面寬度Clr之差E的絕對值小于1，即-1＜E＝Cl-Clr＜1；葉根第三齒齒頸寬度L12＞4；葉根底部寬度L34＞1；輪緣底部寬度L56＞5.如果不滿足以上幾個條件，就會出現不合理的葉根輪緣結構，甚至導致葉根和輪緣實體部分的交叉以及網格劃分時出錯.因此，為了保證計算過程順利進行，必須先對狀態變量進行判定，在得到合理的幾何結構之后，再進行建模和分析計算.

圖3 變量示意圖Fig.3 Schematic diagram of critical variables

（3）目標函數

以葉根和輪緣處最大等效應力的較大者Fmax作為目標函數.

經過最初幾次試算，逐步排除不合理的葉根輪緣尺寸，得到了各變量較合適的取值區間，見表2.

表2 變量的取值范圍Tab.2 Range of variables

4 結果與討論

4.1 原始葉根輪緣計算結果

首先對原始尺寸的葉根輪緣模型進行常溫下（20 ℃）和高溫下（工作溫度534 ℃）的應力分析.圖4和圖5分別為葉片和輪緣在低溫和高溫下的等效應力云圖.由圖4和圖5可見，常溫下葉根處最大等效應力為154 MPa，輪緣處最大等效應力為169 MPa；高溫下葉根處最大等效應力為233MPa，輪緣處最大等效應力為198 MPa.高溫下葉根和輪緣的最大等效應力較常溫下分別增大了51.3% 和17.2%，增幅較大.

圖4 葉片的等效應力圖Fig.4 Equivalent stress contours of the blade

圖5 輪緣的等效應力圖Fig.5 Equivalent stress contours of the rim

另外，分別統計了常溫和高溫下葉根各齒面的承載狀況，見表3.由表3可知，常溫下壓力面葉根齒和吸力面葉根齒的總承載在50%左右，兩個面承載均勻，并且壓力面和吸力面各齒的承載比例均在1／3左右，說明在常溫下該樅樹型葉根各齒接近等強度，設計合理；而在考慮了高溫引起的熱變形后，壓力面和吸力面各齒的承載極為不均，第3齒的承載大大升高，而第1齒的承載則下降較多，如壓力面第1齒承載變為10.05%，而第3齒承載變為55.91%，第3齒的承載已經超過了壓力面總載荷的一半，并達到了第1 齒承載的5 倍以上，葉根的承載狀況惡化.

表3 葉根齒面承載分布Tab.3 Load distribution on root teeth %

由以上數據可以看出，對于調節級葉片，常溫下設計良好的樅樹型葉根輪緣，在高溫下由于熱變形的影響，葉根各齒的承載極為不均，導致產生大應力.葉根處的這種應力分布和承載狀況為葉片的正常工作帶來很大的安全隱患，因此有必要在考慮高溫熱變形影響的前提下對調節級樅樹型葉根輪緣結構進行優化.

4.2 基于模式搜索算法的優化結果

4.2.1 優化過程

設置模式搜索算法的第一次迭代網格尺寸為1.當迭代次數超過100、網格尺寸小于1×10-4或目標函數殘差小于1×10-4三者滿足其一時，優化過程結束.

圖6～圖8顯示了模式搜索的具體過程.圖6給出了每次迭代得到的最佳函數值.由圖6可以看出，經過大約4次成功搜索后，目標函數值已經接近最佳值.從第15步開始，函數值開始趨于平穩，逐步逼近最優值185.58MPa.圖7為每次迭代的網格尺寸.每次成功的迭代后網格尺寸增加，失敗的迭代后網格尺寸減小，本次計算最終的網格尺寸為6.1×10-5，滿足收斂條件，優化結束.圖8給出了每次迭代對應的計算次數，本次分析共計算了306次.可以看到，當目標函數趨近于最優值時，由于網格尺寸越來越精細，需要多次調整才能找到合適的網格尺寸，所以每次迭代的計算次數增加.

圖6 目標函數變化曲線Fig.6 Variation curve of objective function

4.2.2 結果分析

圖7 網格尺寸變化曲線Fig.7 Variation curve of mesh size

圖8 每次迭代計算次數變化曲線Fig.8 Calculation counts in each iteration

圖9 優化前后葉根和輪緣等效應力圖Fig.9 Equivalent stress of root and rim before and after optimization

圖9給出了優化前后葉根輪緣的等效應力云圖.對比圖9（a）和圖9（b）可以看到，優化前后葉根輪緣的應力分布大體相同，但最大等效應力值有較大變化.優化前，葉根最大等效應力為232.84MPa，位于壓力面第三齒圓角處，優化后，葉根最大等效應力為186.64 MPa，同樣位于壓力面第三齒圓角處，但比優化前降低了19.84%，優化效果明顯；輪緣處最大等效應力在優化前為198.46 MPa，優化后為185.58 MPa，均位于壓力面第三齒圓角處，優化后比優化前降低了6.48%，優化效果也較好.同時，優化后葉根和輪緣處的最大等效應力相近，承載狀態改善，有利于葉片安全工作.

此外，為了解葉根型線上應力的整體變化趨勢，給出了優化前后葉根中截面曲線a-b（圖10）上的等效應力分布（圖11）.同時，考慮到葉根最大應力常出現在壓力面第三齒圓角處，給出了優化前后葉根壓力面第三齒圓角處曲線c-d（圖10）上的等效應力分布（圖12）.在圖11中，按由a至b的路徑，曲線中共出現5處峰值，除中間小峰值為葉根底部的應力外，其余四處依次表示吸力面第二齒圓角、吸力面第三齒圓角、壓力面第三齒圓角、壓力面第二齒圓角處的應力值.從圖11可以看出，優化前后沿葉根型線的應力分布總體趨勢不變，在葉根齒圓角處應力值較大，第三齒圓角處應力較第二齒圓角處大.優化前，壓力面第三齒圓角處應力明顯高于其他部位；優化后，壓力面第三齒圓角處應力大幅減小，第二齒圓角處應力也有一定程度減小，吸力面第二齒和第三齒圓角處應力變化較小，另外各齒圓角處的應力差別減小，趨于均勻化，滿足樅樹型葉根等強度設計的要求.從圖12c-d曲線的應力分布可以看出，優化后壓力面第三齒圓角處整體應力水平明顯降低.

圖10 曲線a-b 和曲線c-d 位置示意圖Fig.10 Diagram of curves a-b and c-d

圖11 優化前后a-b 曲線上等效應力分布Fig.11 Equivalent stress distribution along curve a-b before and after optimization

圖12 優化前后c-d 曲線上等效應力分布Fig.12 Equivalent stress distribution along curve c-d before and after optimization

整體看來，經過優化，葉根輪緣處的應力水平大幅降低，承載狀況得到改善，優化效果明顯.

表4給出了優化前后設計變量的數值.

表4 優化前后設計變量數值Tab.4 Values of design variables before and after optimization

5 結 論

（1）通過有限元分析軟件Ansys的參數化語言APDL，建立調節級樅樹型葉根輪緣結構的參數化分析模型，使用優化算法控制優化過程，進行葉根輪緣的結構優化，能獲得相關尺寸的最優值，大大降低葉根輪緣的應力水平.此套方法可行且有效，具有良好的工程實用價值.

（2）對于調節級葉片，常溫下設計良好的樅樹型葉根輪緣由于在高溫下受熱變形的影響，其應力分布和齒面承載狀況變得惡劣，存在安全隱患.

（3）對于本文的調節級樅樹型葉根輪緣模型，選取8個特征尺寸作為設計變量、葉根輪緣處最大等效應力為目標函數，采用模式搜索算法進行優化控制，對調節級樅樹型葉根輪緣結構進行優化分析，獲得了葉根輪緣的優化型線.優化后，葉根最大等效應力下降了19.84%，輪緣最大等效應力下降了6.48%，葉根各齒圓角處應力趨于平穩.整體優化效果良好，提高了葉片運行的可靠性.

［1］王立清，蓋秉政.含表面裂紋T 型葉根應力強度因子的數值 計 算［J］.動 力 工 程，2008，28（2）：169-171，220. WANG Liqing，GAI Bingzheng.Numerical calculation of stress intensity factors for blade’s T-root part with a surface crack ［J］.Journal of Power Engineering，2008，28（2）：169-171，220.

［2］范小平，周振杰，李曦濱.新型135MW 汽輪機調節級葉片開發設計［J］.東方電氣評論，2010（1）：33-35. FAN Xiaoping，ZHOU Zhenjie，LI Xibin.Design of new governing stage blade for 135 MW steam turbine［J］.Dongfang Electric Review，2010（1）：33-35.

［3］吳厚鈺.透平零件結構和強度計算［M］.西安：西安交通大學出版社，2007.

［4］渦輪機教研室.用有限元方法討論樅樹形葉根的合理結構［J］.西安交通大學學報，1977，11（2）：59-74. Turbine Teaching and Research Section.Research on the proper structural of a turbine blade with fir-tree root by finite element method［J］.Journal of Xi’an Jiaotong University，1977，11（2）：59-74.

［5］李家其，賈立言，徐光輝，等.叉型葉根的改進設計研究［J］.汽輪機技術，1992，34（4）：20-23. LI Jiaqi，JIA Liyan，XU Guanghui，etal.Research on improvement design of fork-shaped blade root［J］.Turbine Technology，1992，34（4）：20-23.

［6］邢譽峰，諸德超.航空發動機渦輪盤榫槽的形狀優化設計［J］.航空學報，1995，16（4）：488-491. XING Yufeng，ZHU Dechao.Shape optimum design of engine’s turbodisk tenon-grooves［J］.Acta Aeronautica et Astronautica Sinica，1995，16（4）：488-491.

［7］PAPANIKOS P，MEGUID S A，STJEPANOVIC Z.Three-dimensional nonlinear finite element analysis of dovetail joints in aeroengine discs［J］.Finite Element in Analysis and Design，1998，29（3／4）：173-186.

［8］MEGUID S A，KANTH P S，CZEKANSKI A.Finite element analysis of fir-tree region in turbine discs［J］.Finite Elem and Design，2000，35（4）：305-317.

［9］SONG W B，KEANE A，REES J，etal.Turbine blade fir-tree root design optimization using intelligent CAD and finite element analysis［J］.Computer and Structures，2002，80（24）：1853-1867.

［10］馬全海.面向對象的隨機有限元與接觸問題的可靠性形狀優化設計研究［D］.南京：南京航空航天大學能源與動力學院，2002.

［11］趙海.渦輪榫頭榫槽結構設計方法研究［D］.南京：南京航空航天大學能源與動力學院，2005.

［12］姚利兵，莫蓉，劉紅軍，等.基于強度約束的葉片榫頭參數化設計［J］.航空制造技術，2007（8）：93-95. YAO Libing，MO Rong，LIU Hongjun，etal.Strength constraint based parameterized design of turbine blade serration［J］.Aeronautical Manufacturing Technology，2007（8）：93-95.

［13］楊敏超，孫蘇亞.渦輪榫頭榫槽的結構優化設計［J］.航空動力學報，2010，25（8）：1876-1882. YANG Minchao，SUN Suya.Structural optimization of turbine tenon／mortise［J］.Journal of Aerospace Power，2010，25（8）：1876-1882.

［14］林香，黃致建，郝艷華.燕尾型榫連接件的結構優化設計［J］.制造業信息化，2009（10）：61-63. LIN Xiang，HUANG Zhijian，HAO Yanhua.The optimization design of the dovetail tenon joints［J］.Mechanical Engineer，2009（10）：61-63.

［15］王勖成.有限單元法［M］.北京：清華大學出版社，2003.

［16］陳衛東，蔡萌林，于詩源.工程優化方法［M］.哈爾濱：哈爾濱工程大學出版社，2006.

［17］博弈創作室.APDL參數化有限元分析技術及其應用實例［M］.北京：中國水利水電出版社，2005.

［18］梁醒培，王輝.基于有限元方法的結構優化設計［D］.北京：清華大學出版社，2010.

［19］吳興遠.模式搜索法在最優化問題中的應用［J］.軟件導刊，2009，8（8）：122-123. WU Xingyuan.Application of pattern search in optimization problem［J］.Software Guide，2009，8（8）：122-123.

［20］鄭李鵬.基于APDL 的聯合循環汽輪機轉子熱應力計算及啟動優化的研究［D］.浙江：浙江大學能源工程學院，2011.

［21］馬環宇，楊敏強，李勇.300 MW 汽輪機調節級的溫度場分析［J］.東北電力大學學報，2007，27（2）：14-17. MA Huanyu，YANG Minqiang，LI Yong.Temperature distribution analysis of governing stage for 300 MW steam turbine［J］.Journal of Northeast Dianli University，2007，27（2）：14-17.

［22］許素強，夏人偉.結構優化方法研究綜述［J］.航空學報，1995，16（4）：385-393. XU Suqiang，XIA Renwei.Methods of structural optimization：an overview［J］.Acta Aeronautica et Astronautica Sinica，1995，16（4）：385-393.