高壓渦輪導葉傳熱數值計算與分析

王建華，葛寧

(南京航空航天大學 能源與動力學院，江蘇 南京 210016)

0 引言

在現代燃氣渦輪發動機中，高壓渦輪導葉與燃氣熱交換率最高。因渦輪進口溫度遠超過金屬葉片的熔點溫度，故必需大量使用冷卻氣體，這就要求在設計過程中準確地預測渦輪葉片表面溫度分布。丁林[1]搭建了渦輪葉片測溫平臺，利用光學掃描和輻射測溫方法對葉片表面溫度分布進行測量。近年來，研究人員對氣熱耦合方法進行了廣泛研究。LUO J等[2]研究了渦輪導葉在3種典型工況下的氣動和傳熱特性，采用不同的湍流模型計算分析了葉中截面處溫度分布和外換熱系數分布。褚云會等[3]對跨音速渦輪數值模擬，分析了網格量和湍流模型對計算結果的影響。CROCE G[4]將開發的程序和商業軟件進行了耦合，即流體域和固體域使用不同的網格，通過插值程序將流固交接面上的信息進行傳遞。郭兆元等[5]采用耦合和非耦合方式對單級渦輪機葉片數值計算，表明兩種計算方法對流場結構影響不大，只有耦合計算能準確地預測溫度分布。雖然耦合計算正確地設置了流固交接面上邊界條件，但是，在計算復雜的內部冷卻結構時，這種方法在工程應用上受到限制。因此，在研究冷卻問題時，有學者提出用外換熱系數定量描述對流換熱現象，這也是工程上初步設計冷卻系統的常用方法。

本文首先在典型工況下對C3X計算域進行數值敏感性驗證。然后采用非耦合方法計算分析了壁面溫度對外換熱系數(HTC)的影響。采用耦合方法對比了不同湍流模型下HTC分布。最后對改進的換熱量預測方法進行可行性分析并對比不同計算方法下的HTC。

1 求解方法及驗證

采用課題組自主研發的NUAA-TURBO三維數值求解器和ANSYS-CFX求解器求解雷諾平均Navier-Stokes方程(RANS)。通過Numeca對C3X葉型建立計算域網格，如圖1所示，其中進、出口為H型網格，葉片表面建立O型網格。由于熱邊界層受網格影響較大，為了正確模擬傳熱，葉片表面的y+值<1，網格節點數達2×107。計算域邊界條件如下：來流總壓P0=3.218×105Pa，來流總溫T0in=783K，出口馬赫數Maisout=0.89，湍動量Tu=8.3%，湍流長度尺度為0.005m。

圖1 計算域網格及邊界條件

1.1 網格無關性驗證

首先改變沿葉片輪廓周圍邊界層中的O型網格節點，使用NUAA-TURBO求解器來進行網格無關性驗證。在保證離壁距離相同的情況下，比較了3種不同的網格節點。如圖2所示，在溫度比TR=Tw/T0in=0.65時，計算了30、40和50個節點下葉片表面的無量綱換熱量q/qref分布。無量綱換熱量的計算方法見式(1)，換熱量q通過壁面附近流體分子之間導熱計算，即q=-k(?T/?n)。C為葉片軸向弦長，k為進口溫度條件下的導熱系數。從整體上看，3種網格節點下的計算結果差異很小，在之后的計算中選擇網格節點為40。

(1)

圖2 TR=0.65時不同邊界層節點HTC分布

1.2 非耦合計算

通過NUAA-TURBO求解器的定量計算來分析TR對HTC的影響。計算域的邊界條件如上所述，在不同工況下，僅僅改變壁面溫度邊界條件。假設空氣為理想氣體，理想氣體定律和薩瑟蘭黏性定律來求解氣體的狀態和熱物理性質。從圖3中可以清楚地看出葉片壁面溫度對HTC的影響。在葉片尾緣的相同位置處，在TR=0.60情況下的HTC比TR=0.95情況下高30%左右。同時，在葉片前緣的駐點區域，TR對HTC的影響不如尾緣區域明顯。這也表明上游歷史效應對邊界層的影響。在葉片前緣附近，由于邊界層剛形成，對壁面的溫度影響較小，但是，隨著邊界層在葉片表面的發展，相比于邊界層外部流動，溫度不同的流體在邊界層內部逐漸積累并且不斷影響壁面溫度。這表明當地壁面溫度條件影響HTC，同時受冷卻表面的外換熱系數還與上游邊界層發展歷程相關。

圖3 外換熱系數HTC依賴于TR

1.3 耦合計算

在非耦合計算中只是考慮外部氣動影響。正如在數值計算時，都是給定等溫邊界條件，這就改變了葉片表面和流體域之間的真實交界面邊界條件。為了排除這種任意邊界條件的影響，采用CFX求解器中的熱-固耦合模型數值計算。圖4給出了耦合計算域，外部流體域與非耦合計算下的邊界條件一致，葉片內部固體域材料選擇熱導率較低的ASTM310型不銹鋼，在葉片內部有10個冷卻通道，每個冷卻通道進口給定半徑、進口質量流量、進口靜溫和湍流度。這里選擇文獻[6]中4422工況。

圖4 耦合計算域

圖5顯示了不同湍流模型下，葉片表面無量綱靜壓分布，計算結果表明湍流模型對氣動特性的影響較小。在壓力面，從葉片前緣到x/Cax=-0.5附近，壓力幾乎保持不變，這個階段流體加速緩慢，葉片表面沒有分離，邊界層內的靜壓梯度為0，葉片表面的壓力分布受湍流模型的影響很小。從壓力面中部到尾緣附近，壓差增大，流體急劇加速。在吸力面，從前緣到x/Cax=0.4附近，靜壓快速下降到進口壓力的一半。隨后流體開始減速，到尾緣附近，流體又經歷了緩慢地加速的過程。

圖5 不同湍流模型下葉片表面靜壓分布

如圖6所示，不同湍流模型下換熱系數與實驗值對比。從整體上看，邊界層內部的轉捩流動對HTC有較大影響。在吸力面前緣沿葉片表面下游，不同湍流模型對外換熱系數的預測都有不同程度的偏差。高雷諾數湍流模型(k-ε類)主要模擬湍流狀態，使用半經驗公式簡化邊界層流動，計算所得的湍流強度比較高，同時在吸力面前緣附近，激波-附面層相互干擾的現象極為明顯，對計算結果有較大影響。盡管低雷諾數湍流模型(k-ω類)計算所得的湍流強度較低，但是仍然反映邊界層內部的湍動特性。對于吸力面前緣附近的層流區域，無法準確地分辨。如果在k-ωSST湍流模型上使用轉捩模型，預測的HTC與實驗值吻合較好。

圖6 不同湍流模型下HTC對比

2 非線性方法可行性分析

在典型非耦合計算中假設：氣動性能決定換熱性能。對于給定的葉型結構和流動條件，在牛頓冷卻定律q=h(T-Taw)中，有兩個由氣動決定的未知量(h，Taw)。由于這兩個未知量不受壁溫影響，h值能正確地表示為每一點處的曲線斜率，可將其改寫成式(2)的形式。

(2)

當q與Tw呈現明顯的非線性關系時，曲線斜率表現為在當地變化。式(2)改寫成式(3)。此時式(3)中的外換熱系數h可以解釋為微分算子dq/dTw的有限差分近似，對于給定的壁溫Tw，在極小的溫差ΔT范圍內，曲線的局部斜率能正確地表示外換熱系數。

(3)

為了確定式(3)中ΔT的取值大小，在TR=0.70條件下，選取3K、5K和10K進行驗證研究。為了在圖中顯示清晰，去掉了尾緣回流區域內和轉捩附近HTC分布的峰值點。如圖7所示，在所選取的3個溫差下，計算所得外換熱系數基本相同，因此，在下面兩點法的計算中選擇5K。

圖7 3種溫差下葉中截面外換熱系數分布

根據非耦合計算結果中的分析，非線性方法主要考慮了HTC與Tw呈現線性關系，即HTC=h0+h1Tw。不考慮絕熱壁溫隨壁面溫度的變化。具體表達式如下：

q=(h0+h1Tw)(Tw-Taw)

(4)

在3個壁溫下，求解式(4)中的未知變量h0,h1和Taw。這些未知變量確定后，就可以根據式(4)確定未知Tw下的換熱量分布以及對應壁面溫度下的HTC。未知變量也是葉片表面上的每個網格點的常數。這種當地校正方法與相關經驗關系式有明顯地區別。求解葉片表面HTC分布的常規做法是需要不同壁溫下兩個CFD計算解。該方法只需要3個CFD解，計算成本低，而預測精度得到質的提高。

在溫度比為0.95、0.75和0.65條件下，使用NUAA-TURBO求解器得到每個網格點的換熱量來求解式(4)。用非線性方法可以預測溫度比為0.6～0.95范圍內的換熱量。圖8比較了非線性方法的計算結果與數值計算結果。這種方法可確保在整個Tw范圍內，預測結果與CFD數值計算結果的匹配非常好。為了突出HTC獨立于Tw的假設所帶來的誤差，還使用傳統的兩點法預測了溫度比分別為0.7、0.8和0.9下的換熱量。如圖9所示，結果表明兩點法無法體現HTC隨壁面溫度的變化，從而在換熱量的預測中帶來誤差。

圖8 CFD解與非線性方法的換熱量對比

圖9 CFD解與兩點法的換熱量對比

3 外換熱系數分布

圖10對比了葉片表面平均溫度TR=0.84情況下，兩點法、非線性方法以及耦合計算下的HTC分布。在校正過程中，兩點法在轉捩附近出現了尖峰值，與兩點法相比，非線性方法預測的外換熱系數值更加穩定，這種穩定性主要是對不同壁面穩定范圍內的多點曲線擬合，而傳統的兩點法只是基于單個壁面溫度上局部區域內的差分近似。如圖10所示，非線性方法在轉捩位置處與實驗值和耦合計算有些誤差，從整體上看，可以較為準確地預測外換熱系數。

圖10 TR=0.84條件下外換熱系數分布對比

4 結語

本文以典型高壓渦輪平面葉柵C3X為例，采用NUAA-TURBO求解器和CFX求解器主要研究高壓渦輪葉片外換熱系數與壁溫的關系。得出以下結論：

1)非耦合計算時，壁溫對葉片前緣處HTC的影響較小，但是在葉片尾緣有明顯的非線性。特別地，在TR=0.60與TR=0.95兩種情況下，HTC相差達到30%左右。耦合計算時，帶轉捩的k-ωSST湍流模型能較為準確地預測外換熱系數分布。

2)基于牛頓冷卻定律的非線性形式，采用非線性方法來計算不同壁溫下外換熱系數。與傳統方法相比，可以有效地對HTC-Tw依賴性進行校正。非線性方法是在給定壁面溫度的情況下獲得外換熱系數，這種基于當地的而不是全局的校正方法，提高了計算精度。