趙合陽 白廣忱 費成巍 邢偉震

(北京航空航天大學 能源與動力工程學院，北京100191)(中國人民解放軍空軍93705部隊，遵化064200)

葉尖間隙直接影響著航空發(fā)動機的效率、可靠性、安全性和耗油率［1］.其合理設計和控制一直是研制高性能、高可靠性航空發(fā)動機亟待解決問題之一［2］.渦輪盤徑向變形對葉尖徑向間隙有重要影響，盤的合理設計分析對葉尖間隙設計和控制具有重大意義.發(fā)動機一個工作循環(huán)中，輪盤徑向變形受機械載荷和熱載荷等方面諸多因素影響，隨工作狀態(tài)的不同而變化，合理選擇徑向變形的分析方法是進行渦輪盤設計和葉尖徑向運行間隙設計與控制的基礎［3－4］.

目前國內(nèi)外許多研究機構和學者對渦輪盤徑向變形進行了數(shù)值分析［1－8］.這些研究成果都基于確定性分析，即采用確定的參數(shù)來確定發(fā)動機工作時的渦輪盤徑向變形.由于該方法沒有考慮影響渦輪盤等構件變形的各方面因素的隨機性，只靠留出裕度的方法來保證葉尖間隙的安全性，因而具有很大的盲目性.事實上，影響渦輪盤變形的諸多因素都具有很明顯的隨機性，應該考慮各種隨機因素，從概率的角度研究非線性動態(tài)概率分析，才能更客觀描述渦輪盤變形規(guī)律，改善葉尖徑向間隙設計和控制的合理性.

因此，基于有限元和響應面方法［9－11］，考慮材料屬性的非線性、熱載荷和離心載荷的動態(tài)性，嘗試對輪盤徑向變形進行非線性動態(tài)概率分析.

1 模型建立及計算方法

1.1 有限元模型

選取某型航空發(fā)動機一級高壓渦輪盤為研究對象，輪緣前端與渦輪軸以圓柱表面定心，用徑向銷釘鏈接，工作時前后均有空氣對輪盤側面進行冷卻.渦輪盤的幾何形狀是軸對稱的，可認為渦輪盤承受的外載荷和約束情況是軸對稱的.簡化榫槽、銷釘孔等結構，建立渦輪盤有限元模型，如圖1所示.

圖1 渦輪盤有限元模型

根據(jù)發(fā)動機冷卻系統(tǒng)特點和發(fā)動機實測溫度，輪盤外緣A1、輪盤前側與封嚴圈環(huán)結合處A2、輪盤中心A3等處的溫度邊界可作為第1類邊界條件處理，渦輪盤前、后側面邊界B1，B2，B3可用熱交換的第3類邊界處理.

1.2 概率分析方法

渦輪盤非線性動態(tài)概率分析的基本思想為:

1)選取某航空發(fā)動機典型飛行剖面參數(shù)［2］，將發(fā)動機從地面啟動—慢車—起飛—爬升—巡航這一段［0，215 s］作為計算范圍，取12個關鍵點作為計算點，其飛行載荷譜如圖2所示.

2)考慮材料屬性的非線性和載荷的動態(tài)性，進行渦輪盤變形熱-固耦合分析，計算出渦輪盤徑向變形規(guī)律，并將渦輪盤最大變形量(即全局最大值)作為概率設計的輸出響應.

3)合理選取隨機變量，對隨機參數(shù)進行確定性擬合試驗，再對渦輪盤徑向變形輸出響應進行精確求解，獲得足夠多的樣本值，利用這些樣本值擬合輸出響應與隨機變量之間的響應面模型.

4)將響應面模型代替有限元模型，基于Monte Carlo法(MC法)對渦輪盤徑向變形進行概率分析和靈敏度分析.

圖2 計算采用的載荷譜

假設每一次渦輪盤徑向變形的非線性動態(tài)分析的輸出響應 Y 與輸入?yún)?shù) X=［X1，X2，…，Xr］之間的響應面模型為

式中，a0，bi，cij(i=1，2，…，r;j=1，2，…，r)為待定系數(shù).

假設系統(tǒng)要求渦輪盤最大徑向變形量為δ，則極限狀態(tài)函數(shù)為

由式(2)可知:H＜0時為失效模式;反之為安全模式.若各隨機變量相互獨立且均值和方差矩陣為

則

其中，E(·)和D(·)分別表示均值和方差函數(shù).

若 H(X)符合正態(tài)分布，則可靠度［11］為

由式(5)可得各個隨機參數(shù)的靈敏度:

2 計算結果與討論

2.1 確定性分析

本文考慮離心載荷和熱載荷的動態(tài)性(如圖2所示)和導熱系數(shù)和膨脹系數(shù)的非線性(如表1所示)，利用熱-固耦合的分析方法對渦輪盤徑向變形進行非線性動態(tài)概率分析.首先，結合導熱系數(shù)對其進行非線性動態(tài)熱載荷分析;然后，將熱分析換為結構分析，將計算得到的溫度場作為溫度載荷加入結構分析模型中，在熱載荷和機械載荷的耦合作用下進行徑向變形量計算.經(jīng)計算得渦輪盤徑向變形隨時間變化規(guī)律如圖3所示.

表1 不同溫度下的導熱和膨脹系數(shù)(參考溫度20℃)

圖3 輪盤徑向變形量隨時間變化規(guī)律

由圖3可以看出:在發(fā)動機啟動—慢車—起飛爬升過程中，輪盤的變形量呈增大趨勢;在加速爬升時(即 t∈［165 s，190 s］時)輪盤的變形量達到了最大值;由起飛進入巡航狀態(tài)時，渦輪盤變形量又開始稍微減小.經(jīng)分析，當t=178 s時(燃氣溫度最高、轉(zhuǎn)速最大的時刻)渦輪盤變形量達到全局最大值:1.236 mm，此時輪盤中各節(jié)點的最大變形量同時達到最大，其中渦輪盤溫度分布和徑向變形云圖如圖4和圖5所示.

圖4 t=178 s時渦輪盤溫度分布云圖

圖5 t=178 s時輪盤徑向變形量云圖

由圖4可知:渦輪盤的溫度場不但沿徑向變化，也沿軸向變化.輪緣處的溫度最高，沿徑向渦輪盤溫度逐漸降低，在盤心處取得最小值;由圖5可知:盤的徑向變形在銷釘固定處最小，在盤緣處最大.

2.2 概率分析

選取渦輪盤不同位置的溫度和對流換熱系數(shù)(根據(jù)各對象與燃氣溫度的換熱特點［2］計算)，以及轉(zhuǎn)速ω、材料密度ρ、彈性模量E和泊松比υ為隨機變量，其數(shù)字分布特征如表2所示.假設所有參數(shù)均服從正態(tài)分布且相互獨立.

將表2中的隨機變量分布特征導入有限元模型中，用MC法對各隨機參數(shù)進行抽樣，得到一系列的樣本點，再用這些樣本點對有限元模型進行非線性動態(tài)確定性分析，得到281組試驗樣本，其中輸出響應Y的仿真歷史見圖6.

再利用Box Behnken矩陣法［11］擬合響應面函數(shù)式(1)，得到響應面函數(shù)系數(shù)，進而建立響應面模型，如式(7)所示(忽略系數(shù)小于10－5的項).

表2 輪盤的隨機變量選取

圖6 有限元模型輸出響應Y仿真歷史

將渦輪盤徑向變形的響應面模型代替有限元模型，利用MC法進行1萬次抽樣仿真，其中Y仿真歷史曲線圖和直方圖如圖7和圖8所示.由圖8可知:Y滿足正態(tài)分布，均值為1.2361 mm，方差為 0.030131 mm.

根據(jù)式(2)～式(6)可得渦輪盤徑向變形量在不同設計值δ下的可靠度值，如表3所示.

由表3可知:當δ=1.33mm時，失效數(shù)為15，可靠度R=0.99852，基本滿足工程需要.因此，在渦輪盤變形設計時，建議δ=1.33 mm.

圖7 響應面模型輸出響應Y仿真歷史

圖8 輪盤徑向變化直方圖

表3 不同設計值δ下的可靠度

同時，也可以對盤的變形進行逆概率分析，通過逆概率分析，其結果如表4所示.

表4 逆概率分析結果(部分)

由表4可以看出，通過逆概率分析，可以得到每個可靠度下的各個隨機參數(shù)的極限值，為渦輪盤設計與優(yōu)化提供量化參考.

2.3 靈敏度分析

根據(jù)式(5)和式(6)可計算出各隨機變量對渦輪盤徑向變形的靈敏度和影響概率，如圖9和表5所示.

圖9 隨機變量對渦輪盤徑向變形的靈敏度分析

表5 渦輪盤徑向變形靈敏度分析結果

在圖9中，最重要的隨機輸入變量(靈敏度最大)在最左邊，其他依次向右排列.在表5中，靈敏度有正負之分，正表示隨機輸出變量隨輸入?yún)⒘空兓?，負則表示反變化.由圖9和表5可得以下結論:

轉(zhuǎn)子轉(zhuǎn)速ω對輪盤徑向變形的影響最大，靈敏度達到0.9，影響概率為0.536.符合本文采用的二代發(fā)動機輪盤的工程實際情況;輪盤溫度Ta1，Tb2，Ta3和 Tb1對渦輪徑向變形也有非常重要的影響.泊松比、彈性模量、導熱系數(shù)等對輪盤變形的影響很小.所以，在對渦輪盤徑向變形設計與優(yōu)化中，首要考慮轉(zhuǎn)速、輪緣溫度的影響.同時，也為葉尖徑向運行間隙設計與控制提供了量化參考.

3 結論

本文考慮材料屬性的非線性和熱載荷、離心載荷的動態(tài)性，對輪盤徑向變形進行了概率分析，得出如下結論:

1)通過渦輪盤徑向變形的確定性分析，得到了輪盤徑向變形量及其隨時間的變化規(guī)律，并且當t=178 s時(燃氣溫度最高、轉(zhuǎn)速最大時刻)渦輪盤變形量達到全局最大值:1.236 mm.

2)通過概率分析和逆概率分析，建立了輪盤變形概率分析的響應面模型;得到了輪盤徑向變形分布特征:均值為 1.236 1 mm，方差為0.030131 mm，且服從正態(tài)分布;還得到不同設計值δ下的徑向變形的失效數(shù)和可靠度，建議δ=1.33 mm，此時可靠度 R=0.998 52，符合設計要求.另外還得到不同可靠度下的各隨機參數(shù)的極限值，為渦輪盤徑向變形設計和優(yōu)化，以及葉尖徑向運行間隙的設計和控制提供了依據(jù).

3)通過靈敏度分析，得到了各隨機參數(shù)的靈敏度和影響概率，找到了影響徑向變形變化的主導因素(轉(zhuǎn)速 ω)和主要因素(溫度 Ta1，Tb2，Ta3和Tb1)，為渦輪盤徑向變形設計和優(yōu)化，以及葉尖徑向運行間隙的設計和控制提供了量化參考.

References)

［1］Hennecke D K，Trappmann K.Turbine tip clearance control in gas turbine engines［R］.NASA N83-29254，1983

［2］Lattime S，Steinetz B.Turbine engine clearance control systems:current practices and future directions［R］.AIAA 2002-3790，2002

［3］漆文凱，陳偉.某型航空發(fā)動機高壓渦輪葉尖間隙數(shù)值分析［J］.南京航空航天大學學報，2003，35(1):63－67 Qi Wenkai，Chen Wei.Tip clearance numerical analysis of an aero-engine HPT［J］.Journal of Nanjing University of Aeronautics ＆ Astronautics，2003，35(1):63－67(in Chinese)

［4］Pilidis P，Maccallum N R L.Models for predicting tip clearance changes in gas turbines［R］.NASA N83-29258，1983

［5］NSSA Glenn Research Center.HPT clearance control［R］.NASA/CR-2005-213970，2005

［6］豈興明，樸英，祝劍虹.某型航空發(fā)動機高壓渦輪葉頂間隙三維數(shù)值分析［J］.航空動力學報，2008，3(5):904－908 Qi Xingming，Pu Ying，Zhu Jianhong.3-D numerical analysis of the tip clearance of an aero-engine high pressure turbine［J］.Journal of Aerospace Power，2008，3(5):904－908(in Chinese)［7］Eeom Y S，Yoo K S，Park J Y.Reliability-based topology optimization using a standard response surface method for three-dimensional structures［J］.Structural and Multidisciplinary Optimization，2011，43(2):287－295

［8］李昌，韓興.基于響應面法齒輪嚙合傳動可靠性靈敏度分析［J］.航空動力學報，2011，26(3):711－715 Li Chang，Han Xing.Analysis of reliability sensitivity for gear engagement based on response surface methods［J］.Journal of Aerospace Power，2011，26(3):711－715(in Chinese)

［9］Kypuros J A，Melcher K J.A reduced model for prediction of thermal and rotational effects on turbine tip clearance［R］.NASA/TM-2003-212226，2003

［10］Annette E N，Christoph W M，Stehan S.Modeling and validation of the thermal effects on gas turbine transients［J］.Journal of Engineering for Gas Turbines and Power，2005，127(3):564－572

［11］Hakan E.Design sensitivity analysis of structures based upon the singular value decomposition［J］.Computer Methods in Applied Mechanics and Engineering，2002，191(5):3459 －3476