基于快速預測模型的輪盤瞬態變形影響

馬曉健，黃敬杰，徐如雪，智紹強

（中國航發沈陽發動機研究所，沈陽110015）

0 引言

為了達到設計點的工作狀態，航空發動機需要優化加速或減速路徑，從而引入了瞬態性能的概念。在瞬態條件下，主流道參數急劇變化，氣體與固體之間發生劇烈熱量交換，稱為熱浸潤。例如從慢車到最大狀態的瞬態加速，發動機機體必須適應新的穩態工作溫度，通常會吸收30%的燃料能量[1]。熱浸潤的作用結果之一是使發動機結構尺寸如葉尖及封嚴間隙發生瞬態變化[2]，進而影響發動機性能。例如，高壓渦輪葉尖間隙增加1%會導致其效率降低2%，進而導致發動機燃油消耗率增加。此外，與穩態相比，在瞬態加速過程中發動機工作線更容易接近喘振線，帶來了較大的喘振風險[3-4]。轉子在離心力作用下，通常加速時更容易與機匣涂層發生碰摩，而過度碰磨會造成葉片開裂、涂層脫落等故障，增大發動機失效風險[5]?？梢娙绻軌颢@得瞬態葉尖間隙變化規律，可以為優化發動機性能特性和降低發動機故障風險提供一種途徑[6]。

獲得葉尖間隙變化規律的方法包括實時監測法和數值計算方法，前者主要應用于發動機產品研制后期，后者則多用于研制中前期。數值計算方法中功能最為強大的是NASA的NPSS計劃[7]，該計劃基于準確的3維部件模型以及完整的CFD和FEA數據資源，缺點在于需要劃分網格、添加載荷和熱邊界條件，因而費時費力；有些則不基于物理模型，直接根據FEA模擬結果、間隙實測數據等，利用公式擬合發展出經驗性預 測方 法[8]；Fiola[9]、Kypuros等[10]和Melcher等[11]則采用簡單幾何模型以及經驗公式，模擬換熱過程及機械變形，達到快速預測葉尖間隙變化的目的。

航空發動機研制經驗表明，在發動機方案設計過程中發現葉尖間隙變化帶來的不利影響，可以大大減少發動機設計后期試驗投入[12]。為了在方案階段實現對葉尖間隙的快速評估，本文根據高壓渦輪的結構特點，吸收了現有快速預測方法的優點，改進了渦輪盤瞬態變形快速預測方法。

1 方法改進及模型假設

1.1 方法改進

針對E3-GE-FPS發動機第1級高壓渦輪建立了考慮換熱作用的葉尖間隙模型[13]，如圖1所示。該模型中高溫燃氣流經葉片、輪盤盤緣以及機匣內襯內表面，引自壓氣機出口的空氣用于冷卻輪盤側面以及機匣，自機匣外側垂直引入的沖擊氣流用于實現主動間隙控制[14]。

圖1 E3-GE-FPS第1級高壓渦輪[13]

從圖中可見，影響渦輪盤變形的形式有熱變形和離心力帶來的變形。通過對Kypuros和Melcher模型[10]輸入參數進行修正，基于有限體積法發展了輪盤變形預測模型。具體改進如下：

（1）原模型采用等厚圓盤表示輪盤，為了更接近真實的等強度設計輪盤剖面，本文采用參數化模型定義變厚度輪盤。

（2）原模型只考慮溫度及轉速影響，實際不同冷熱氣流量對于輪盤換熱也有影響，本文據此對換熱模型進行改進，增加流量考慮。

（3）原模型對于輪盤外緣進行絕熱處理，為接近真實情況，本文增加了高溫燃氣對輪盤盤緣的加熱考慮。

（4）原模型將輪盤看作單一單元，本文采用有限體積法對輪盤分網，并按次序對各單元進行換熱計算，模擬換熱過程。

（5）原模型任意時刻輪盤具有惟一溫度及材料屬性（彈性模量及熱膨脹系數），本文增加了不同單元不同屬性的考慮。

改進后的輪盤模型如下：通過8個半徑（R1~R8）以及4個寬度（W1~W4）定義參數化模型，如圖2（a）所示。對輪盤進行分網：軸向分為3列，徑向分成n行，從而獲得由3n個圓環組成的輪盤離散模型。定義冷熱氣溫度、流量以及對流換熱路徑從而建立輪盤熱浸潤模型，如圖2（b）所示。

圖2 輪盤模型

1.2 模型假設

本文研究策略既不同于需要龐大數據資源、費時費力的NPSS計劃，又不同于靠數學擬合形成的經驗方法，而是基于輪盤模型，遵循傳熱學、力學定律，針對發動機方案設計特點和需求，通過合理假設發展一種輪盤變形快速預測方法。所做假設如下：

（1）在發動機方案設計階段，尚未形成完整的輪盤特征，因此不考慮葉片安裝槽以及前后連接結構等細節特征是合理的。

（2）事實上輪盤表面的對流換熱系數受流體類型、環境以及離心力作用而呈現區域各異分布，但真實的對流換熱系數有賴于試驗數據，為簡化計算，本文參考Kypuros和Melcher模型，采用恒定的對流換熱系數。

（3）在真實發動機中，對流換熱與熱傳導同時發生，本文為適應程序算法需要，將其假設為：在1個時間間隔內，對流換熱先進行，穩定后再進行熱傳導直到穩定。

（4）任意時刻輪盤沿軸向均存在溫度梯度，本文假設對于相同半徑的3個圓環單元，取均值作為等效溫度，用于機械變形計算。

（5）為簡化計算，假設任意時刻，溫度在各單元圓環內部均布，因而材料屬性（彈性模量及熱膨脹系數）對于各圓環只有單一值。

2 渦輪盤瞬態模擬方法

2.1 輪盤溫度計算方法

2.1.1 熱傳導機理

固體的1維熱流由傅里葉傳導定律給出[15]。將此方法擴展到多個維度，可模擬1個單元同時受四周所有單元的熱傳導作用?；谔摂M穩態溫度假設，提出了單位時間Δt內某單元受四周單元傳導的熱量變化[8]

式中：K'為單元與某一相鄰單元綜合傳導率；Tn、Tv分別為相鄰單元或目標單元溫度；SumKT為相鄰單元溫度與傳導率乘積之和；SumK為相鄰單元傳導率之和；m為固體質量；Cpf、Cps、Cp為固體等壓比熱容。

2.1.2 對流換熱機理

液體流經固體表面的對流換熱模型如圖3所示。假設單位時間內液體質量由流量Wf確定。考慮固體表面至中心發生的熱傳導作用，可以獲得單位時間Δt內液體與固體中心發生的熱量交換[8]

圖3 熱對流

式中：Z為整個系統的熱特性

式中：A為接觸面積；h為對流換熱系數；k為固體熱傳導率；X為固體表面距中心距離；Cpf，、Cps分別為液體、固體的等壓比熱容；Tf0、Ts0為液體、固體初始時刻溫度。

2.1.3 瞬態溫度計算流

將上述方程應用到渦輪盤熱浸潤模型（圖2（b））中：式（1）用于輪盤所有單元熱傳導，式（2）用于模擬氣流與輪盤外部單元的換熱。

基于此可以建立輪盤的熱浸潤計算流程。對于某一瞬態操作，冷熱氣與輪盤之間產生溫差，輪盤外部單元吸收或釋放，溫度變化，在輪盤內部形成1個新的溫度場；由于輪盤內部單元之間存在溫差，必然發生熱傳導，最終獲得新的溫度分布。將各單元溫度傳遞給下一時刻，從而建立瞬態條件下的輪盤溫度計算命令流。

2.2 輪盤變形預測方法

2.2.1 應力計算

在方案設計階段，輪盤應力分析仍然可以在有限等厚圓環近似模型上進行，2個相鄰圓環A和B的定義如圖4所示。圖中：σ為應力；r為半徑；w為圓環厚度；下標i和o表示內、外表面，r和h表示徑向和周向。

圖4 輪盤簡化計算模型

在相鄰圓環公共面，根據力的平衡及周向變形協調，徑向及周向應力差值分別為

對于已知尺寸的自由旋轉圓環，其內、外表面應力存在固有關系。因此如果輪盤盤心或盤緣應力已知，利用式（4）、（5）可逐步計算各圓環公共面應力值，獲得輪盤應力分布。

基于此建立輪盤應力的迭代運算流程。假設盤心應力（盤心通常是自由表面，因此徑向應力為0），逐級向外計算應力直至盤緣，對比計算的盤緣徑向應力與由于葉片離心力帶來的徑向應力是否一致，否則調整盤心周向應力，反復迭代直到找到最佳值。

2.2.2 變形計算

假定輪盤材料為彈性材料，不考慮塑性變形，同時不考慮輪盤沿軸向變形。根據輪盤近似模型可知，將各圓環的變形量從最小半徑處累加直至最大直徑圓環，即可獲得輪盤累積半徑變化，即盤緣徑向變形，表示為

對于最小直徑的圓環單元，已知內表面徑向應力為0，外表面由于受輪盤其他部分的離心作用而產生向外的拉應力。其外表面變形的計算可以通過計算自由旋轉圓環的變形疊加外部受載的靜止圓環變形獲得[16]

式中：ρ為密度；ω為角速度；ν為泊松比；E為彈性模量；q為外部壓力。

對于其他圓環，徑向應變可表示為

式中：ε為應變；T為圓環等效溫度；α為某溫度下的熱膨脹系數。

假設單個圓環的徑向應變沿半徑方向呈線性變化，對半徑進行積分，則單個圓環徑向的變化為：

式中：Ac，Bc為描述圓環應變線性變化的常數。

3 計算結果及分析

3.1 計算流程圖

根據前述分析，結合渦輪盤熱浸潤以及變形的主要要素，形成了渦輪盤瞬態變形計算流程，如圖5所示。其中，棕色模塊代表性能分析，需要定義基準發動機以及瞬態歷程并輸出氣動參數及轉速；在熱浸潤模擬（藍色模塊），通過建立數學物理模型模擬換熱作用，獲得輪盤上的不同區域溫度分布；綠色模塊表示力學計算過程，輸入轉速、溫度、材料屬性以及葉片離心力，通過不斷迭代計算應力分布獲得最優解，進而預測輪盤盤緣變形。

圖5 渦輪盤瞬態變形流程

基于該流程，在Matlab中編制了預測代碼，通過設定模型參數、氣動輸入可以快速獲得溫度、應力、變形的預測結果。

3.2 發動機基準模型

引用文獻s[10]中輪盤、葉片的相關參數用于本文模型算法。為適應算法特點，對部分參數進行了調整：在保持盤緣寬度、內外徑不變的條件下，將等厚剖面改為典型輪盤剖面；初始輪盤溫度為高壓排氣溫度θ3的1/2；考慮材料屬性隨溫度的變化，所采用的Inco?nel718材料屬性隨溫度的變化趨勢[17]及其線性表達如圖6所示；此外根據算法定義計算時間間隔為0.5 s，假設葉片數量為50。綜合上述參數假設為基準參數設定。

圖6 In718材料彈性模量E和熱膨脹系數α隨溫度變化

建立以該渦輪盤為對象的單軸發動機性能基準模型，其設計點轉速對應狀態性能參數，見表1。利用性能設計軟件Gasturb[18]模擬油門由慢車保持（換算轉速為0.7），然后增大到設計點轉速對應狀態（換算轉速為1.0）并繼續保持的瞬態過程，獲得關鍵位置氣動數據隨時間的變化關系。

表1 單軸發動機設計點轉速對應狀態主要性能參數

3.3 變形預測結果

采用等厚空心輪盤模型和基準參數設定，相應的輪盤盤緣變形歷程如圖7所示。分為3個變化階段：在慢車和設計點的穩定狀態，輪盤變形緩慢增加，由于轉速并未變化，變形來源于輪盤溫度升高，表明本文方法模擬了輪盤不斷從環境中吸收熱量；當轉速增大時（計算步長150起始），輪盤在離心力作用下變形加快，瞬態歷程結束時變形值最大，為1.45 mm，與文獻[11]所采用方法獲得的結果（最大2.11 mm）對比可知，初始變形相同，整體變形趨勢相似，而本文方法預測的變形結果較低；將時間間隔從基準（0.5 s）調整到1.5、2.5 s，變形趨勢保持不變，而變形數值增大，接近或超過文獻[11]的預測結果。

圖7 輪盤變形計算結果對比

可見利用本文改進的方法可以較好地預測變形趨勢，而通過定義合理的參數如時間間隔，可以進一步對變形預測結果進行優化。

3.4 參數影響分析

對于渦輪盤熱浸潤模型，影響換熱的主要參數包括對流換熱系數、熱傳導率、時間間隔以及初始輪盤溫度，不同參數設定影響換熱效果，進而決定輪盤溫度。直接影響應力應變的參數則包括葉片數量、輪盤剖面以及材料屬性。在基準參數設定上，調整參數數值或類型（見表2），利用單變量法分析各參數對換熱及應力應變的影響。

表2 輪盤換熱及變形影響因素

3.4.1 換熱分析

設定基準參數，瞬態歷程結束時輪盤徑向溫度分布如圖8所示。從圖中可見，盤緣溫度較高，盤心溫度較低，中間呈現漸變的溫度梯度[19]，證明了模型中燃氣加熱輪盤盤緣的有效性。

圖8 對流換熱系數對溫度分布影響

從圖中還可見，在不同的對流換熱系數下獲得的溫度水平也不同。較低的換熱系數，表明流體對輪盤較難傳遞熱量，輪盤溫度水平較低，反之亦然。

對不同熱傳導率的換熱過程模擬結果如圖9所示。從圖中可見，高熱傳導率使溫度更高的盤緣更易于向溫度低的盤心傳導熱量，降低了徑向溫度梯度，使其變得更平。

圖9 輪盤熱傳導率對溫度分布影響

假設時間間隔變化時發動機單位時間間隔瞬態氣動參數保持不變。在此基礎上，將模型換熱的時間間隔從0.50 s調整為0.25 s和0.75 s，對溫度的分布影響如圖10所示。

圖10 計算時間間隔對溫度分布影響

從圖中可見，隨時間間隔的增大，溫度水平升高。這是由于在1個大的時間間隔，換熱更充分，進而導致溫度升高，反之亦然。

對流換熱是由流體與金屬固體之間的溫度差驅動的，因此固體的初始溫度對傳熱過程有很大影響。本文分析了不同初始溫度（壓氣機出口溫度θ3、θ3/2、以及ISA室溫（θ0=15℃））對溫度分布的影響，如圖11所示。

圖11 輪盤初始溫度對溫度分布影響

從圖中可見，輪盤初始溫度可以決定輪盤（尤其是盤心位置）的溫度水平，改變溫度梯度。在工程實際中，造成這一情況出現的主要源于發動機的起動時機，即是冷起動或是熱起動。

3.4.2 應力應變分析

本文在建立的方法中引入了葉片離心力作用。通過調整葉片數量（25、50、75）定義不同離心力，相應的應力分布如圖12所示。

圖12 葉片數量對周向應力分布影響

從圖中可見，不同葉片離心力使輪盤周向應力發生偏置，但差別較小。且葉片越多，離心力越大，周向應力也越大，反之略微降低。

典型剖面輪盤與等厚輪盤周向應力分布的對比如圖13所示。從圖中可見，典型剖面使輪盤周向應力水平顯著降低，尤其是在盤心處，降低比例為44%。

圖13 輪盤剖面類型對周向應力分布影響

模型采用的輪盤材料屬性（彈性模量和熱膨脹系數）具有溫度相關性。應變計算結果如圖14所示。與恒定的材料屬性設定對比可知，徑向應變整體有所降低。

圖14 材料屬性對徑向應變分布影響

4 變形預測與結果分析

4.1 不同換熱參數導致的輪盤變形

將表2中參數對應的輪盤溫度、應力應變計算結果作為輸入對輪盤變形進行預測，采用相同油門（轉速）變化歷程、考慮不同換熱參數（表2中第1～4項）獲得的輪盤變形結果如圖15所示。

圖15 不同換熱參數導致的輪盤瞬態變形

結合表2和圖15具體分析如下：

（1）與基準相比，改變輪盤初始溫度θi（4A，4C），使輪盤初始變形出現約±65%變化，原因是將初始溫度賦值給輪盤，自盤心到盤緣累積的熱變形導致初始變形發生改變；隨計算步長增加這種偏差逐漸減小，至瞬態歷程結束時變為±23%，主要是由于換熱受溫度差值影響，較低的初始輪盤溫度反而具有較高的熱量交換，結果使3種情況預測結果隨計算步長逐漸逼近。

（2）根據對流換熱機理可知，熱量與換熱系數h以及時間間隔Δt呈正相關，即在高換熱系數和大的時間間隔下，輪盤獲得的熱量更多，因而溫度更高（圖15）。將這2個參數各自增大50%，瞬態歷程結束時輪盤累積變形增大約10%（1C，3C），反之則有所減?。?A，3A）。

（3）由前述分析可知，輪盤材料熱傳導率k決定徑向溫度梯度，進而影響輪盤變形。而傳導率的改變通常使部分單元溫度升高，而另一部分單元溫度降低，2種單元的熱變形相互抵消，使總的變形差別并不明顯（15）。熱傳導率增加50 %（2C），瞬態歷程結束時輪盤累積變形僅增大約3%，反之則有所減?。?A）。

4.2 直接參數作用導致的輪盤變形

采用相同油門（轉速）變化歷程、考慮直接參數（表2中第5～7項）作用獲得的輪盤變形結果圖16所示。

圖16 直接參數作用導致的輪盤瞬態變形

結合表2和圖16具體分析如下：

（1）與基準相比，增加葉片數量NoB（5C），增大了輪盤承受由葉片提供的離心力，根據前述分析，使輪盤整體應力水平升高，結果輪盤產生更大的變形。使葉片數量改變50%，瞬態歷程結束時偏置量約為±3%（圖16）；除葉片數量外，單個葉片質量、葉片質心位置也可導致離心力變化，進而影響變形。

（2）改用等厚輪盤（6A），瞬態歷程結束時，使變形預測結果增加8.8%。這是由于采用等厚輪盤，輪盤應力水平較高；而采用將盤心加厚的等強度設計，可以降低整體應力水平，提高輪盤抗變形能力。

（3）相比可變材料屬性，恒定彈性模量E和熱膨脹系數α（7A）的設定將使預測結果增加8.3%。這是由于應力主要由離心力和溫度作用導致，因此2種情況下應力基本相同。但由于彈性模量影響彈性變形，熱膨脹系數決定熱變形，2項共同作用導致變形發生變化。

5 結論

（1）本方法是對現有簡化預測方法的發展，主要改進包括：可用于不同類型不同尺寸輪盤剖面；考慮了冷熱氣流量對換熱的影響；增加燃氣對盤緣加熱的模擬；對輪盤分網并分區域進行換熱模擬；考慮材料屬性隨溫度的變化。

（2）對流換熱系數、時間間隔影響輪盤整體溫度水平；熱傳導率影響盤心與盤緣間的溫度梯度；而初始溫度決定盤心及中部溫度。葉片數量變化使輪盤周向應力發生偏置；與等厚輪盤相比，典型剖面能顯著降低周向應力水平；可變材料屬性設定，使徑向應變有所降低。

（3）綜合影響變形的主要參數，其中初始輪盤溫度θi影響最為顯著，使輪盤初始變形出現約±65%偏置；換熱系數h以及時間間隔Δt使變形發生約10%偏差；采用等厚輪盤或恒定材料屬性的設定，導致預測結果出現8.8%和8.3%的偏差。葉片數量的調整可導致3%的變形變化；熱傳導率k改變對于預測結果影響較小。

（4）本方法較好地預測了加速過程中的渦輪盤盤緣變形瞬態歷程，通過定義參數可進一步優化預測結果。本方法基于結構模型，遵循傳熱學、力學基本定律，通過調整參數定義，可針對不同發動機、不同工作狀態建立滿足特殊需求的輪盤變形預測模型，具有較好的工程適用性。