起動過程渦輪葉片瞬態溫度場的準穩態處理研究

常 勝，杜 強，王 沛，柳 光，劉 軍，肖向濤

(1. 中國科學院 工程熱物理研究所，北京 100190；2. 中國科學院大學 工程科學學院，北京 100049；3. 中國科學院 輕型動力創新研究院，北京 100190；4. 中國科學院大學 航空宇航學院，北京 100049)

燃氣渦輪發動機起動、停車等過程(也稱過渡態)中來流溫度和壓力的急劇變化，將會引起燃氣與渦輪葉片換熱的大幅度動態變化，同時在金屬導熱影響下，葉片的溫度場也會處于非穩定狀態。過渡態下異常的溫度分布會導致渦輪轉子承受較大的局部熱應力及熱變形，進而使葉片發生熱疲勞現象，嚴重影響渦輪部件工作壽命，因此有必要開展發動機過渡過程渦輪葉片瞬態溫度場的研究，為葉片可靠性分析提供依據。

過渡過程涉及流動穩定和熱穩定兩類問題，通常，建立系統的熱穩定(即建立穩定的溫度場)所需的時間要遠遠大于建立流動穩定(即建立穩定的速度場)所需要的時間[1]。當高溫燃氣流過渦輪葉片通道后，其熱擾動穿透固體的深度有限，此時葉片的熱響應將慢于氣動響應，應考慮金屬熱慣性帶來的溫度滯后效應[2]。過渡態葉片溫度的滯后會引起熱變形的滯后，隨之對葉尖間隙產生影響，進一步引起渦輪效率的下降[3]，嚴重時會導致葉尖發生刮磨，因此須要全面地預測過渡態下金屬葉片的溫度場。由于發動機整個飛行包線涉及眾多過渡過程，快速準確的溫度預估方法就顯得尤為重要。

作為目前燃氣渦輪發動機設計的關鍵技術，流熱耦合數值模擬在渦輪換熱領域得到了廣泛的應用。Takahashi[4]應用商業化軟件Fluent 5對燃氣渦輪發動機起動和停車過程中第一級渦輪轉子葉片進行了瞬態流熱耦合計算，總結了來流變化條件下渦輪葉片的換熱規律：前緣和尾緣溫度變化最劇烈，葉根處溫度梯度增加較快，但對葉片內部的溫度并沒有過多的關注。Rossette[5]在Takahashi的工作基礎上，研究了發動機起動過程中帶熱障涂層的高溫合金轉子葉片的對流換熱和導熱，分析了葉片中截面上熱障涂層和葉片基體的溫度梯度隨時間的變化規律，此外，將非定常項引入固體域的流熱耦合模型可以用來評判影響固體熱疲勞機理的熱慣性。為實現有限元分析(Finite Element Analysis，FEA)和計算流體力學(Computational Fluid Dynamics，CFD)的高效耦合，Sun[6]將“SC89”插件程序應用于發動機四級低壓渦輪盤腔的流熱耦合計算中，數值模擬采用英國羅-羅公司開發的軟件Hydra，使用程序調用一組典型時刻CFD算例的熱邊界條件，為有限元分析軟件SC03提供輸入，在確保FEA和CFD計算之間溫度和熱流連續的前提下，實現二者的快速耦合。該過程相當于采取了一種“準穩態處理”的方法，對非穩態中典型時刻的流場進行一組穩態數值模擬，利用中間程序調用各工況下的流動參數，再進行固體的有限元分析，達到縮短計算時間的目的。因為流場響應的時間尺度遠小于固體熱傳導的時間尺度，所以相對固體的熱慣性，流動的瞬態效應可以忽略。

國內在渦輪葉片瞬態溫度場預測方面也開展了相關研究，錢惠華[7]針對某型渦扇發動機高壓渦輪導向葉片的裂紋故障，采用自編程序計算了葉片各部分的換熱系數，結合有限元分析軟件ANSYS模擬了起動過程50 s內的瞬態溫度場，發現導向器在各瞬時的溫度最高點均在葉片尾緣中部，該區域溫度在9 s時已基本穩定，認為這是尾緣的換熱系數比較大的緣故。徐磊[8]應用一維解析方法和數值方法分別對航空發動機整個飛行循環內薄壁類和軸對稱緣盤類熱端部件進行了熱狀態計算分析，并依次研究了各零部件(葉片、輪盤及機匣)的熱響應，為工程設計提供了一種過渡態葉尖間隙的分析方法。關鵬[9]對渦輪導向葉片進行瞬態流熱耦合計算，研究了來流溫度階躍條件下的葉片外表面及內部冷卻通道的溫度場，然后將溫度結果導入ANSYS求解器，計算了葉片的熱應力分布與模態振型變化。屠秋野[10]以雙軸混排渦扇發動機為研究對象，通過一種非定常熱交換模型分析了高壓轉子各部件的吸熱和放熱特點，得到葉片的時間常數[11]最小的結論。王伯鑫[12]提出一種葉片對流冷卻二維模型，并結合CFX對美國通用電氣公司E3發動機高壓渦輪第一級導葉進行數值模擬，將葉片表面不同位置的平均換熱系數和近壁面燃氣溫度輸入該冷卻模型，實現了葉片壁面溫度分布變化的快速求解。

流熱耦合的計算方法雖然可以通過設置流固交界面的方式將燃氣側的對流換熱、金屬固體內的導熱一并考慮，且對傳熱過程的模擬與實際情形較為接近，但就起動過程渦輪葉片換熱計算而言，模擬全程變化的速度場和溫度場需要耗費大量的計算資源，且計算周期長，不易實現葉片瞬態溫度場的快速求解。國內目前對渦輪葉片熱分析中的邊界條件通常根據經驗公式或者約化模型預估，因此得到的結果與燃氣渦輪發動機的真實情況存在一些偏差，且在燃氣渦輪發動機過渡態葉片換熱問題中，對固體溫度滯后效應的關注還比較少。本文將借鑒上文中提到的準穩態處理的思想，以期對葉片的瞬態溫度場進行快速求解，掌握其熱響應規律，為渦輪葉片的熱分析提供支撐，進而縮短發動機熱端部件的研制周期。

在準穩態處理過程中，為準確計算渦輪葉片的溫度場，必須確定燃氣與葉片之間的換熱關系，這就涉及渦輪葉型外部邊界層的流動與換熱問題。目前邊界層流動換熱計算的方法和程序有很多種，但從葉片型面邊界層計算效果看，還是STAN5[13-14]比較適用。20世紀90年代，西北工業大學的劉松齡、朱惠人將STAN5程序引入國內，并在Hylton[14]的工作基礎上，對計算程序添加了駐點模型、轉捩模型等內容，提高了STAN5程序在渦輪葉片流動換熱計算中的準確性。近年來，劉松齡對STAN5程序作了進一步的改進，包括對駐點附近計算模型的改進、加入Mayle[15]建立的轉捩模型、針對逆壓梯度和不合理速度分布造成計算失敗的現象制定科學的修型應對措施等，大大改善了程序對多種流動狀態下葉片換熱系數求解的適用性，改進后的版本命名為“NPUSTAN7Z”。

本文針對燃氣渦輪發動機起動過程渦輪葉片瞬態換熱問題，以某型燃氣渦輪發動機第一級渦輪葉片為研究對象，結合商業化軟件ANSYS CFX 2019 R2和外換熱程序NPUSTAN7Z，探索了一種對過渡態進行準穩態處理[6]的方法。準穩態處理采用非定常計算的方法求解三維Navier-Stokes(N-S)方程，提取有關流動參數輸入外換熱程序，獲得求解固體域的第三類邊界條件，最后進行非穩態導熱計算，得到過渡態下葉片的瞬態溫度場。同時還對過渡態渦輪葉片進行了瞬態流熱耦合計算，即在流體域和固體域中同時求解能量方程，流固交界面施加通量守恒條件，僅在流體域中求解輸運方程，對起動過程進行了總時間3 s的非穩態計算，最終也得到過渡態下葉片的瞬態溫度場。分析了葉片熱慣性帶來的溫度滯后效應，通過對比準穩態處理和瞬態耦合計算兩種方法得到的渦輪葉片對來流熱響應的差異，驗證了準穩態處理的方法對過渡態渦輪葉片換熱特性預測的準確性，最終實現對葉片瞬態溫度場的快速求解。

1 計算模型及方法

1.1 計算模型及邊界條件

以某型燃氣渦輪發動機在地面全程加速過程中五級低壓渦輪第一級葉片為研究對象，約化后建立渦輪計算模型，包括3個導葉、5個動葉， 如圖1所示。本文主要研究轉子葉片的瞬態溫度場，其幾何參數如表1所示，葉頂間隙取葉身高度的0.7%。

(a) 約化后模型

(b) 動葉進口速度三角形示意圖1 計算模型

表1 低壓渦輪轉子葉片幾何參數

燃氣工質設置為理想可壓縮氣體，其粘性系數μ和熱導率k隨溫度的變化滿足薩瑟蘭(Sutherland)公式[16]，設置定壓比熱cp為溫度的表達式[17]。葉片材料選擇K417G高溫合金，其導熱系數、比熱容等參數隨溫度變化[18]。

對燃氣渦輪發動機地面全程加速的試驗數據進行模化，模化范圍對應試車過程的中間14 s，以模化后來流條件變化率最大的3 s(6.2～9.2)為研究區間，截取第一級低壓渦輪進口總溫Tt、進口總壓pt、出口靜壓p和轉速值N，以起始時刻第6.2 s為參考對各變量進行無量綱處理，得到隨時間變化的邊界條件相對值C/C6.2(C代表Tt、pt、p、N),如圖2所示。

圖2 起動過程3 s內邊界條件的變化曲線

燃氣經過渦輪過渡段后，氣流角沿徑向非均勻分布，根據測試數據，進口各速度分量沿徑向的變化如圖3所示,進口湍流度取5%。

圖3 主流進口沿徑向的速度分量

1.2 瞬態流熱耦合邊界設置

瞬態流熱耦合計算采用商業化軟件ANSYS CFX 2019 R2進行數值模擬，計算域包含導葉流體域、動葉流體域、導葉固體域和動葉固體域四個計算域，而非定常計算域僅包含導葉流體域和動葉流體域兩部分，相對流熱耦合計算，非定常計算不含固體域，圖4(a)、4(b)分別為瞬態流熱耦合計算域和非定常計算域。

流熱耦合計算中進出口邊界條件按照圖2與圖3給出的數據設置，以起始時刻的邊界條件進行非定常流熱耦合計算，將計算結果作為瞬態流熱耦合和準穩態計算的初場。轉靜交界面設置為Transient Rotor Stator混合模型，葉片表面(含葉尖)設置為通量守恒的流固交界面，除交界面外所有壁面均滿足無滑移條件。考慮葉根榫槽內充滿冷卻氣體，3 s加速過程葉根設置為800 K恒定的溫度，如圖4(c)所示，輪轂、機匣(上、下面)為絕熱條件，兩側(左、右面)設置周期交界面。

(a) 耦合計算域

(b) 非定常計算域

(c) 各邊界位置

1.3 準穩態處理計算過程

準穩態處理的計算方法結合了商業化軟件ANSYS CFX 2019 R2和外換熱程序NPUSTAN7Z，分以下幾個步驟來開展：CFX求解流體域、數據預處理、表面換熱系數計算、三維節點插值等，其流程如圖5所示。

圖5 準穩態處理流程圖

1.3.1 非定常計算邊界設置

對過渡態進行準穩態處理，選取起動過程研究區間3 s內的7個典型工況，各工況點間隔約0.5 s，分別對應的時刻是A:第6.2 s、B:第6.699 5 s、C: 第7.200 5 s、D:第7.7 s、E:第8.199 5 s、F:第8.700 5 s、G:第9.2 s，以上述時刻的邊界條件進行非定常數值模擬。非定常計算域如上文描述，設置所有壁面為無滑移絕熱條件，轉靜交界面設置為Transient Rotor Stator混合模型，左右兩側設置周期交界面。

1.3.2 計算表面換熱系數

表面換熱系數的計算是準穩態處理中的關鍵步驟。采用非定常的方法求解時均流場后，提取相關流動參數輸入到外換熱程序中，程序采用流函數坐標變換法，求解邊界層微分方程，計算葉片表面換熱系數，作為邊界條件，進一步再進行葉片瞬態溫度場的求解。

本文提取轉子葉片表面10%、20%、…、90%共9個展向高度上的型線數據輸入NPUSTAN7Z程序求解邊界層方程，如圖6所示為不同時刻外換熱程序計算的部分展向高度上沿型面分布的表面換熱系數，橫軸表示無量綱的型面弧長，0表示前緣，區間(-1,0)表示壓力面，區間(0,1)表示吸力面。由于前緣受到主流燃氣的沖擊，邊界層非常薄，因此具有較高的換熱系數。沿壓力面和吸力面型線至尾緣邊界層逐漸增厚，相應地換熱系數較低。隨著發動機起動過程中來流條件的改變，表面換熱系數也逐漸增加。

1.3.3 三維節點插值

進行渦輪葉片導熱計算時須要在網格節點上加載邊界條件，而葉片表面的網格節點與外換熱程序計算的9條已知數據型線的坐標并不一致，不能滿足三維溫度場計算的需要[19]，進一步還應將9條型線上已知的換熱系數插值到固體域的網格節點上。

(a) 70%展向高度

(b) 50%展向高度

(c) 30%展向高度圖6 隨時間變化的表面換熱系數

本文采用MATLAB編寫了三維彎扭葉片表面節點數據的插值程序，采用K-近鄰算法[20]，按照歐式距離，尋找距葉片表面網格節點最近的已知型線上的3～4個坐標點，以坐標點間的距離確定權重，將相應坐標點上的換熱系數插值到葉片表面節點上。尋找近鄰點時，兩型線中間的節點分別尋找上、下型線上距離最近的2點，葉根及葉頂區域尋找展向10%、展向90%型線上的3點，這樣可以避免已知型線附近的網格節點只尋找了該條型線上的數據，造成局部換熱系數偏差過大。如圖7為初始時刻9條型線上的換熱系數插值出的三維葉片表面所有節點的換熱系數分布圖。

圖7 初始時刻插值后的葉片表面換熱系數分布

1.3.4 求解葉片瞬態溫度場

依次對準穩態處理的工況A至G進行求解流體域、計算表面換熱系數、進行節點插值，得到7個時刻下非穩態導熱計算離散點的第三類邊界條件，接著對節點上相鄰兩個時刻的表面換熱系數和絕熱壁溫取平均，獲得中間6個時刻的邊界條件，計算步驟如表2所示。3 s過渡過程的瞬態溫度場計算分12個區間進行，時間步長設置為7.5×10-4s，初次計算時與耦合計算設置相同的初場，之后每次計算均以上次的結果為初場，轉子葉根處給定800 K均溫，最終得到過渡態葉片內部的瞬態溫度場。

表2 葉片溫度場計算步驟

1.4 網格劃分

所有計算域均采用ICEM CFD劃分結構化網格，為減少插值帶來的誤差，流體域和固體域的交界面上設置了相互匹配(相同數目及分布規律)的網格，在固體域周邊布置O型網格，保證了沿葉片表面網格的正交性。

劃分了350 萬、709 萬、1 001 萬三種不同數目的網格進行網格無關性驗證，各網格劃分方案如表3所示。按照初始時刻第6.2 s的邊界條件做非定常流熱耦合計算，通過對比三種方案得到的瞬態初場的溫度、壓力值，綜合考慮計算精度和所需成本，網格劃分選擇方案2，流熱耦合與非定常計算采用同樣的流體域網格。

表3 網格劃分方案

渦輪動葉的網格細節如圖8所示,壁面第一層網格的尺度為0.001 mm，使得壁面無量綱單位y+小于1，滿足湍流模型的要求。為使流固交界面上網格具有較好的匹配性，固體域也采用O型網格，壁面以內第一層網格尺度設置為0.008 mm，沿葉片型線的網格數目和分布規律相同，即交界面兩側網格密度相同。葉頂間隙內沿徑向布置25層網格，兩端網格尺度設置為0.001 mm，輪轂壁面采取同樣的處理方式。

圖8 動葉流體域與固體域的網格細節

2 數值方法

2.1 控制方程

2.1.1 CFX

數值模擬采用商業化軟件ANSYS CFX 2019 R2，其核心求解器CFX Solver中，固定的笛卡爾坐標系(i，j，k)下，微分形式的連續方程、動量方程和能量方程：

(1)

(2)

(3)

式中：SM為動量源項；SE為能量源項；h*為比總焓。且有

(4)

因為在固體區域內沒有流動存在，可以將對流項和擴散項取消，并且在固體域內只有熱傳導傳熱方式存在，所以在固體區域中能量守恒方程簡化為：

(5)

CFX軟件對流體域和固體域均采用有限體積法求解，在流體域中求解三維可壓縮雷諾時均N-S方程，在固體域中僅求解能量方程。瞬態流熱耦合計算中，固體域和流體域滿足能量方程的同時，在流固交界面上施加通量守恒的條件。

耦合計算和非定常計算中，湍流模型均選擇基于雷諾平均方法的二方程模型——剪切應力輸運(Shear Stress Transport, SST)格式的k-ω模型，求解過程中，控制方程對流項均采用高精度差分格式(High Resolution Scheme)，時間項采用隱式歐拉二階向后差分格式(Second Order Backward Euler Scheme)。

2.1.2 NPUSTAN7Z

外換熱程序NPUSTAN7Z是一個邊界層有限差分的計算程序，通過求解壁面邊界層微分方程，得到渦輪葉片表面的換熱系數。

邊界層內連續方程：

(6)

x方向的動量方程：

(7)

式中：τ為邊界層內流體的粘性應力。

y方向的動量方程：

(8)

總焓形式的能量方程：

(9)

式中：q為單位面積的熱流量；X為單位質量流體的徹體力，對于氣體常可以略去；I為流體的總焓。

外換熱程序中采用的湍流模型有混合長度模型、單方程湍流模型，同時考慮了湍流度修正。其中較為常用的是普朗特(Prandtl)雙方程混合長度模型。

2.2 計算結果驗證

非定常計算的時間步長設置為1×10-5s，此時各典型時刻動葉流體域單個扇區(6.372°)通過周期內包含49～73個時間步，能夠保證計算結果的準確性。如圖9為7個工況計算的動葉出口的溫度監控曲線，受到上游尾跡影響，當動葉出口的溫度表現出明顯的非定常特性時，證明計算已收斂。每次計算收斂后，對200步的流場結果取時均值，用于外換熱程序的輸入。

圖9 非定常計算動葉出口監控點的溫度

圖10 耦合計算葉片前緣監控點的溫度

瞬態流熱耦合計算總時間較長，需要消耗大量的計算資源，選擇了三種時間步長1.5×10-2s、1.5×10-3s、1.5×10-4s進行無關性驗證，對起動過程進行總時間0.5 s的瞬態流熱耦合計算，考察轉子葉片中截面附近前緣一點的溫度值，如圖10所示，可以看到隨著時間步長的增加，監控點的溫度隨時間的推移逐漸起伏變化，這與圖2中來流條件的平順變化規律不符，當時間步長為1.5×10-3s 和1.5×10-4s時，溫度變化趨勢更接近真實情況，但綜合考慮耦合計算的經濟性，最終確定1.5×10-3s的時間步長。

3 結果及分析

3.1 過渡態葉片表面瞬態溫度分布

通過準穩態計算與瞬態耦合計算兩種方法，得到轉子葉片的瞬態溫度場，如圖11為壓力面和吸力

(a) 壓力面瞬態溫度場

(b) 吸力面瞬態溫度場圖11 耦合計算和準穩態計算的瞬態溫度場

面溫度隨時間的變化云圖。在發動機加速3 s時間內，隨著來流溫度、壓力和轉速的增加，轉子葉片表面溫度也隨之上升，并且呈現出前緣和尾緣溫度升高快于葉中、由葉尖至葉根溫度逐漸降低、葉根附近存在較大的溫度梯度等特點。

(a) 70%徑向高度

(b) 50%徑向高度

(c) 30%徑向高度圖12 不同徑向高度型面上的溫度變化規律

對型面弧長坐標進行無量綱化，得到葉片表面不同徑向高度型線上溫度的變化曲線，如圖12所示。隨著時間的推移，葉片表面不同徑向高度上的溫度逐漸升高，前緣和尾緣相對中弦位置溫度升高更快。由于耦合計算的結果受到流場的瞬態效應(尾跡、二次流等)影響，沿弧長的溫度分布相比準穩態的計算結果波動較大。兩種計算方法得到的表面溫度分布存在一定差異，末尾時刻溫度差較為顯著，在70%、50%、30%徑向高度上，二者計算的型線平均溫差為6.7 K、10.1 K、8.5 K，前緣附近的差異尤為顯著，末尾時刻準穩態的結果高于耦合計算14.4 K、18.4 K、16.3 K。這源于兩種計算方法耦合方式的不同，準穩態通過求解邊界層方程得到表面換熱系數進而計算溫度場，流熱耦合計算采用設置交界面上通量守恒的方式得到壁面溫度。但總體上二者溫度變化趨勢接近，說明準穩態處理方法具有一定的準確性。

3.2 葉片內部瞬態溫度場

3.2.1 不同區域的瞬態溫度

為詳細了解葉片內部的溫度變化趨勢，查看葉片30%、50%、70%無量綱徑向高度截面上的溫度值，圖13為7個典型時刻兩種計算方法求解的葉片內部的溫度場云圖。由圖可知，葉片內溫度隨著時間的推移而增加，非穩態導熱對葉片內溫度場產生一定影響。隨著來流條件的改變，熱擾動由壁面逐漸向內傳遞，葉片內部的熱響應慢于葉片表面。此外，轉子葉片榫槽內被冷卻空氣環繞，相對葉身溫度較低，因此葉片內部存在由葉身向葉根的導熱，由葉尖至葉根溫度逐漸降低，進而造成葉片內部沿徑向的熱響應不一致，葉片中下部對燃氣的熱響應相對較慢。

圖13 葉片不同徑向截面溫度云圖

葉片內部的等值線方向接近于橫向，表明葉片內部存在縱向熱傳導，熱量由前緣、尾緣向中弦位置傳遞，沿中弧線方向存在由中段指向前緣、尾緣的溫度梯度，最終在葉片中段形成低溫區。準穩態和耦合計算得到的葉片內部溫度變化趨勢仍然是接近的，其中準穩態處理預測的溫度場略高。

圖14為葉片內部9條特征弧線(如圖15中從壓力面穿過葉片至吸力面的圓弧段)上的溫度隨時間的變化曲線，反映了葉片從壓力側到吸力側的橫向溫度變化規律，從左至右，從上至下依次為10%、50%、90%無量綱軸向位置和70%、50%、30%無量綱徑向高度上的溫度值，其中橫軸表示無量綱的周向弧長，-1表示壓力側，(-1，1)表示葉片內部，1表示吸力側。10%、50%、90%軸向位置可以用來表征葉片前緣、葉片中弦位置、葉片尾緣的換熱特性。

在第6.699 5 s時刻，各個位置上的溫度均已偏離初始值，因此經過0.5 s后熱擾動已穿透整個葉片。各特征位置上的熱響應逐漸加快，如圖中70%徑向高度上前緣附近在開始的0.5 s內，耦合計算(圖14(a)中實線橢圓圈出)的溫度平均增加了1.5 K，末尾的0.5 s內溫度平均增加了10.7 K。又如50%徑向高度上中弦位置在開始的0.5 s內，耦合計算(圖14(b)中虛線橢圓圈出)的溫度平均增加了0.1 K，末尾的0.5 s內溫度平均增加了7.9 K。這是由于起動過程中來流雷諾數的增加帶來葉片表面換熱系數的增加，同時受到來流溫度升高的影響，導致越來越多的熱量傳遞到葉片內部，由此熱響應的也越來越快。

(a) 70%徑向高度

(b) 50%徑向高度

(c) 30%徑向高度圖14 葉片內部不同特征弧線上的溫度值

圖15 葉片無量綱相對位置示意圖

觀察葉片各處的溫度變化規律可以發現，葉片的溫度從初場開始升高，之后約0.5 s時間內，各軸向位置幾乎不存在橫向溫差，隨著時間的推移，葉片壓力側、吸力側溫度差開始增加。如70%徑向高度上前緣附近在第6.2 s和第6.699 5 s時刻的橫向溫差為0 K，在第8.700 5 s和第9.2 s時刻，耦合計算的橫向溫差變為2.0 K和3.0 K。又如50%徑向高度上中弦位置在第6.2 s和第6.699 5 s時刻的橫向溫差也是0 K，在第8.700 5 s和第9.2 s時刻，耦合計算的橫向溫差變為1.7 K和1.9 K，其余位置也具有同樣的現象。起動過程中，由于渦輪葉片吸力面和壓力面流動狀態的不一致，使得兩側的換熱特性不同，導致葉片內部存在一定的橫向溫度梯度，且溫度梯度逐漸增加。

在50%軸向位置上，接近末尾時刻可以觀察到溫度呈中間低兩端高的現象，說明中弧線附近區域對燃氣溫度的熱響應較慢，由于前緣、尾緣的橫向尺寸較小，熱擾動極短的時間內就能將其穿透，因此沒有形成明顯的低溫區。對比不同徑向高度前緣、中弦及尾緣溫度隨時間的變化曲線，發現三者中中弦位置溫度隨時間變化最慢，相對來說，局部熱容量較小的前緣、尾緣對燃氣溫度變化的熱響應比葉片中弦位置更快。準穩態計算的溫度略高，但在溫度變化趨勢上與流熱耦合的結果比較一致，對得比較好。

3.2.2 葉片內部平均溫度及其滯后效應

為定量分析葉片內部的溫度變化趨勢，對30%、50%、70%徑向截面的溫度取平均值，得到不同截面的溫度隨時間的變化曲線，如圖16所示。圖中給出了過渡態加速3 s過程內來流截面(位于葉片上游軸向弦長30%處，如圖7所示)平均的溫度變化，耦合計算的來流溫度升高了93.0 K(非定常計算初、末兩個算例來流溫升93.4 K)，30%、50%、70%徑向高度上溫度分別升高了27.3 K、27.2 K、33.9 K，準穩態計算的溫升分別是35.5 K、37.1 K、38.7 K，二者計算的溫升存在一定差異，但計算結果共同反映了葉片中上部相對中下部熱響應更快的現象，這是由葉片中下部熱傳導帶走更多的熱量導致的。由于金屬葉片本身具有熱慣性，葉片內部不同截面上的平均溫度隨時間的變化率低于來流溫度隨時間的變化率，說明葉片的溫升慢于來流的溫升，因此葉片內部存在溫度滯后效應。

圖16 不同徑向截面平均溫度的變化曲線

準穩態計算與瞬態流熱耦合計算的結果總體上吻合較好，需要說明的是，來流溫度初始值的不同是因為二者的計算域不同，準穩態初溫取來流截面時均溫度的平均值，流熱耦合初溫取來流截面瞬態溫度的平均值，二者計算域及邊界條件不一致，所以初始時刻存在偏差。

以初始時刻的溫度為參考值，定義燃氣溫升百分比和金屬溫升百分比：

(10)

(11)

式中：tg0為燃氣初始時刻溫度；tg為燃氣溫度；Δtg為燃氣溫度相對初始時刻溫度之差；ts0為葉片初始時刻溫度，葉片內部各處初始溫度不盡相同；ts為葉片溫度；Δts為葉片溫度相對初始時刻溫度之差。

圖17為不同徑向高度上面平均的溫升隨時間的變化曲線，圖中給出了過渡態加速3 s過程內來流截面平均的溫升曲線，可以看出，來流溫度升高了11.1%(非定常計算初、末兩個算例來流溫升11.3%)，耦合計算和準穩態計算的30%、50%、70%徑向高度面平均的溫升分別是2.7%、2.7%、3.3%和3.5%、3.6%、3.8%，二者計算的溫升百分比偏差最大處位于葉片中部和下部，分別為0.8%、0.9%，偏差較小，因此認為準穩態處理的方法能夠準確預測過渡態渦輪葉片的瞬態溫度場。

圖17 不同徑向截面的平均溫升

葉片中下部、中部、中上部溫度相對來流溫度分別滯后了8.4%、8.4%、7.8%(耦合)，所以中上部的熱響應相對中下部較快，葉片中下部區域滯后效應較強，考慮是因為葉根溫度相對較低，葉片中下部溫度場受熱傳導的影響更大(觀察準穩態計算的結果，可以得到相同的結論)。

圖18為葉片內部不同位置上的線平均溫度滯后百分比的柱狀圖，受葉片熱慣性的影響，葉片內部不同區域的溫度均會存在滯后，且溫度滯后程度不一樣。聯系上文分析，葉片較薄位置(前緣、尾緣)相對葉片較厚位置(中弦區域)熱響應更快，葉片中下部的溫度滯后效應較強。對比圖18(a)中前緣、中弦和尾緣附近的溫度滯后，可知中弦滯后>前緣滯后>尾緣滯后。對比圖18(a)中耦合計算不同徑向高度上的溫度滯后，可知葉片中下部滯后>葉片上部滯后。

葉片內局部溫度滯后最大的位置在50%軸向弦長處的中下部區域，最大滯后達到9.6%，因為該區域靠近葉根，同時也是葉片中較厚的位置，所以該區域相對來流的熱響應最慢。

準穩態計算結果的溫升普遍偏高，對比圖18(a)與圖18(b)不同位置上準穩態與耦合計算的滯后百分數偏差，前緣附近70%、50%、30%徑向高度上的偏差分別為0.7%、1.2%、0.8%，中弦位置分別為0、0.9%、0.8%，尾緣附近分別為0.1%、0.2%、0.4%，準穩態處理的計算結果雖然在各局部位置上均低估了金屬葉片的溫度滯后，但相對耦合計算的溫度場偏差不大。

(a) 耦合

(b) 準穩態

4 結語

本文分別采用瞬態流熱耦合計算和準穩態處理的方法求解過渡態渦輪葉片的瞬態溫度場，通過比較二者的計算結果，發現準穩態處理與瞬態耦合計算的結果吻合較好，僅略微低估了金屬葉片的溫度滯后，因此認為準穩態處理可以實現過渡態渦輪葉片換熱特性的快速求解。同時開發了三維彎扭葉片的節點數據插值程序，解決了工程實際中渦輪葉片瞬態溫度場計算時的邊界條件銜接問題。分析葉片的瞬態溫度場得到如下結論：

(1) 非穩態導熱會對葉片內溫度場產生一定影響，葉片內部的熱響應慢于葉片表面。不同軸向位置溫度滯后程度的不一致，形成了葉片內部的縱向溫度梯度，葉片中段存在低溫區。吸力面和壓力面流動狀態的不一致導致葉片內部存在一定的橫向溫度梯度。

(2) 熱擾動在較短的時間內已穿透整個葉片，起動過程中，隨著時間的推移，來流溫度與表面換熱系數的共同作用使得葉片內的熱響應逐漸加快。

(3) 葉片本身的熱慣性導致其溫升慢于來流的溫升，在發動機起動過程的3 s內，葉片中截面平均溫度相對來流滯后了8.4%。

(4) 葉片不同軸向位置的厚度不一致，帶來了不同區域局部熱容量的差異，反映到轉子葉片溫度場上的就是中弦滯后>前緣滯后>尾緣滯后。

(5) 由于葉根是低溫區，離葉根越近熱傳導帶走的熱量越多，所以葉片中下部滯后>葉片上部滯后，葉根附近存在較大的溫度梯度。