基于非結構網格的渦輪氣熱耦合數值方法研究

湯超，葛寧，SHENG Chun-hua

基于非結構網格的渦輪氣熱耦合數值方法研究

湯超1，葛寧1，SHENG Chun-hua2

(1.南京航空航天大學能源與動力學院，江蘇南京210016；2.College of Engineering，University of Toledo，OH 43606，USA)

利用課題組自主開發的三維非結構隱式N-S計算軟件CU_Turbo，采用氣熱耦合計算方法，對MarkⅡ內冷徑向渦輪導向葉片、帶氣膜冷卻渦輪導葉MT1的流場和溫度場進行了數值模擬。計算過程中，隱式時間推進中Jacobi?ans矩陣采用對Roe通量的一種近似方法求解。結果表明，計算值與試驗值吻合良好，驗證了氣熱耦合計算方法的實用性和有效性，為渦輪工程設計提供了一種新的計算分析方法；渦輪葉片通道內附面層的不同流動狀態及氣膜冷卻，對當地換熱都有很大的影響。

渦輪葉片；氣熱耦合；氣膜冷卻；隱式時間推進；非結構網格

1 引言

渦輪是發動機中熱負荷和動力負荷最大的部件，采用有效的冷卻措施是發動機安全可靠工作的有力保證，也是降低材料成本的有效措施。近年來，渦輪前燃氣溫度的逐年提高與渦輪的有效冷卻，特別是葉片內部空氣冷卻技術的迅速提高分不開。有研究表明，渦輪葉片設計中，如果預測的葉片溫度超過實際工作溫度28℃，葉片壽命將減半[1]。因此準確估算葉片的傳熱系數和溫度，有助于預防熱腐蝕和設計高效冷卻系統。

渦輪葉片的數值模擬中，常規方法是，用給定的固體表面溫度作為邊界求出葉片表面的換熱系數分布，再以此作為固體域熱傳導計算的邊界條件。該方法中，壁面溫度的給定更多的是依靠經驗，從而導致計算的不確定性增加，加大了預測熱負荷的難度，尤其是對于復雜的幾何體。而實際中，渦輪部件的溫度場與流場相耦合，單方面研究無法為發動機設計提供準確依據，故提出了氣熱耦合方法[1]。該方法不需要用壁面的溫度分布和換熱系數分布作為邊界條件，而是在整場計算中獲得，很好地解決了常規方法中的不確定性，且可適用于任意復雜的幾何體和工作狀態。同時，由于復雜的冷卻系統導致渦輪葉片造型很不規律，相對于結構網格來說，非結構網格求解器對其有更好的適應性。

本文使用課題組自主開發的基于非結構網格的三維計算軟件CU_Turbo，通過對葉型MarkⅡ和MT1進行數值模擬，驗證了氣熱耦合計算方法的實用性和有效性，為渦輪工程設計提供了一種新的計算分析方法。

2 數值計算方法

2.1控制方程和數值方法

笛卡兒坐標系下，守恒變量形式的三維可壓縮流雷諾平均N-S方程的形式為：

式中：Q=[ρ ρu ρv ρw ρet]為守恒變量，F為通量矢量。

CU_Turbo軟件流體域的計算采用格點格式的有限體積法，無粘通量的離散采用Roe格式，Ven?katakrishnan限制器[2]抑制數值振蕩，時間推進采用LU-SGS[3]隱式方法并加入當地時間步長和隱式殘值光順來加速收斂，高階精度由高斯-格林和最小二乘梯度重構方法實現，邊界條件依據特征波原理[4]，加入Spalart改進的S-A湍流模型把握近壁區域流動。固體域采用熱傳導控制方程，計算方法與流體域保持一致。

2.2非結構網格的氣熱耦合方法

在流固交界面上，必須滿足交界面溫度的連續性和交界面熱通量的連續性兩個物理條件[5]：

對于如圖1所示的非結構網格，由式(2)可得：式中：kf、ks分別為流體域和固體域的熱傳導系數，Tf、Ts分別為流體域和固體域邊界相鄰邊對應點的靜溫，Δnf、Δns分別為流體域和固體域相鄰邊對應點到邊界的距離，m、n分別為交界面上點在流體域和固體域的相鄰邊數，i、j分別指交接面上點相應在流體域和固體域的相鄰邊，wi、wj分別為相應邊的權重且表達式為

由式(3)可求得壁面溫度Tw，再以此作為流體域和固體域的定溫邊界條件分別進行計算。

2.3隱式推進數值方法

時間推進采用LU-SGS隱式算法，把式(1)有限體積法離散并采用向后差分可得：

式中：Vi為控制體i的體積，Rn+1i為時間第n+1層上殘值且代表控制體i的第ij個控制面上的無粘通量，Fv代表該控制面上的粘性通量，n?代表控制面外法矢，Ωij代表控制面面積。

式(5)為一非線性系統，直接求解非常麻煩，通過泰勒展開將其線化可得：

對于上式中每個時間步的系數矩陣，可分解為對角陣D、上三角陣U和下三角陣L：

其中D為純對角線矩陣，L、U分別對應編號小和編號大的鄰居單元：

將式(8)帶入式(7)可進一步改寫為：

向前推進

向后推進

求得的修正項Δq更新本時間步：

2.4通量Jacobians矩陣的計算方法

對流通量Jacobians矩陣的求法是隱式求解的關鍵，區別于Steger-Warming矢通量分裂的解析方法[6]求解Jacobians矩陣，本文采用對Roe通量的一種近似方法求解。該方法簡單可靠且易于編程實現。本課題組史萬里等對以上兩種方法做過比較研究，評估了Jacobians矩陣計算的可靠性和實用性[7]。控制面ij的Roe通值表達式為：

式中：帶上標⌒的項采用修正的Roe平均變量計算，認為該項為當地值并假設為常數。求微分后可得任一控制面對于界面左右兩側的對流通量Jacobians矩陣的表達式：

粘性通量Jacobians矩陣的求解采用薄層假設，求解過程參見文獻[8]。

3 算例驗證及分析

CU_Turbo軟件已經過湍流平板氣動驗證、跨聲流凸包氣動驗證、層流平板氣動熱力耦合驗證[9]，且計算值與試驗值吻合良好。本文對徑向渦輪導葉MarkII和帶有6排圓柱形氣膜孔的渦輪導葉MT1進行詳細的計算和分析，來驗證該軟件的實用性和有效性。

3.1MarkII

MarkII葉片是有10個冷卻孔的渦輪直葉片，葉高76.20 mm，弦長136.22 mm。圖2給出了葉片及計算域幾何尺寸，其中左、右表面分別為進口和出口，上、下兩邊界面為周期性邊界面，中間為帶冷卻孔葉片，頂端和底端面為對稱邊界。Hylton已在接近工程運行條件下做過大量試驗[10]，本文模擬的是Hyl?ton試驗中的5411工況，進出口參數為：總壓332 000 Pa，總溫811 K，背壓171 282 Pa，進口馬赫數0.19，出口馬赫數1.04。

圖2 MarkⅡ葉片幾何尺寸及工況Fig.2 Geometric configuration and boundary conditions of MarkII vane

采用cooper的方法劃分計算域網格，在貼近葉片表面的壁面處使用邊界層網格，使得第一層y+≤1以保證精確模擬邊界層內流動。采用商用軟件Gambit生成網格，通過CU_Turbo接口軟件GMT[9]實現數據結構轉換。

圖3為馬赫數等值線分布云圖。從圖中看，由于高進、出口壓比和吸力面收縮通道的影響，在吸力面前緣附近流動有一個強加速過程，在x cx=0.45 (x cx為橫坐標與軸向弦長之比)位置出現一道強激波，這也將影響到該區域的熱負荷；激波后邊界層出現分離，之后流動減速到亞聲速范圍，在更遠的下游流體再次加速，在尾緣附近出現一道較弱激波。而在壓力面一側，馬赫數增加較緩慢。

圖3馬赫數等值線分布Fig.3 Mach number contours

圖4 所示為溫度等值線分布云圖。可見，溫度最低點發生在第2和第3冷卻孔之間，由于這一區域冷卻孔位置比較集中，冷卻效果要比其它位置的好；吸力面一側，氣流溫度在激波前降低到局部最小值，在激波后又迅速抬升到當地最大值，溫度最高的位置在葉片尾緣。原因為尾緣部分冷卻孔數量沒有前緣部分多，冷氣流量小，加之葉片尾緣的直角形狀，氣體流動復雜，各種渦系相互摻混，所以溫度在這一區域達到最高。

圖4溫度等值線分布Fig.4 Temperature contours

圖5 給出了葉片表面中徑處的性能曲線。從圖5(a)中可知，在吸力面，從滯止點開始壓力迅速降低，在x cx=0.45處壓力減小到最小值，激波與壓力分布的位置相對應。

圖5(b)中，滯止點溫度達到局部最大值。沿壓力面一側，由于第2和第3冷卻孔間冷卻效果較好，以及葉片輪廓由凸轉凹、凸面流動換熱量不斷減小使得溫度降低[11]，壓力面溫度從滯止點開始逐漸降到局部最小值。沿吸力面一側，由于與壓力面相同的原因，溫度先降低到局部最小值，激波后邊界層發生分離誘導轉捩，導致換熱量急劇上升，溫度也隨之迅速增加，接近尾緣時由于葉片變薄和冷卻不足致使溫度進一步上升。溫度計算值最大誤差小于3%，壁面平均誤差小于1%。

圖5(c)中，換熱系數h定義為h=qw(T0-Tw)，qw為葉片表面熱流密度，T0為氣體來流溫度。可見，溫度影響換熱，換熱系數與表面溫度分布呈現出很強的相關性，計算值與試驗值趨勢基本吻合。駐點后壓力面和吸力面的熱傳導系數急劇下降，吸力面激波處誘導分離轉捩，導致熱傳導系數開始迅速上升。由于計算的溫度邊界層比實際的厚，使得壓力面換熱系數的計算值比試驗值偏低，靠近尾緣區域的誤差，還可能是由于當地附面層流動受相鄰葉片壓力面強激波影響所致。

圖5 MarkII葉片表面中徑處性能曲線Fig.5 Performance profiles of MarkII vane

3.2MT1

MT1葉片的冷卻系統由2個圓柱形冷卻腔和6排圓柱形氣膜孔組成，冷卻氣自冷卻腔上端進入、經氣膜孔流出，與主流摻混在葉片表面形成氣膜，冷卻腔下端封閉，兩腔在吸力面均只有1排氣膜孔，在壓力面有叉狀排列的2排氣膜孔[12]。

圖6所示為葉片計算域及網格，左側為進口邊界，右側為出口邊界，計算域上、下端設為對稱邊界，前、后為周期邊界，所有壁面均設為定溫邊界。各邊界具體參數為：主流來流總壓P*=4.60×105Pa，主流來流總溫T*=444 K，前冷卻腔進口總壓P*1=6.28×105Pa，后冷卻腔進口總壓P*2=4.88×105Pa，冷卻腔進口總溫T=286.5 K，Tw=288 K。

圖6 MT1葉片計算域及網格Fig.6 Computation domain and grids of MT1 vane

在貼近葉片表面的壁面處使用邊界層網格，使得第一層y+在1的量級以保證精確模擬邊界層內流動。氣膜孔進、出口處也采用了與壁面邊界層一樣的加密方式，以保證氣膜孔與主流通道交界面上網格過渡平滑。流體域網格170萬。

圖7示出了葉片表面等熵馬赫數和葉片表面中徑處努賽爾數的分布情況。從圖7(a)中可看出，計算值與試驗值吻合較好，吸力面后端(x=30 mm)有一道弱激波，氣膜孔對吸力面等熵馬赫數分布的影響較大，波動明顯強于壓力面。

圖7(b)中，努賽爾數的定義式為Nu=hd λ，d為特征長度，λ為流體域熱傳導系數。對比有無氣膜孔時的努賽爾數分布，可清楚看出氣膜孔對葉片表面換熱的影響，在每組氣膜孔后換熱強度均有較大幅度的降低，有效減小了葉片表面的熱負荷。對比帶氣膜的計算值與試驗值，壓力面由于計算得出的溫度附面層偏厚導致努賽爾數的計算值稍偏低，吸力面峰值處則略高于試驗值，總體趨勢吻合良好。同時，可觀察到在吸力面約x cx=0.30處換熱量急劇增加，這是由于此處發生邊界層轉捩所致。

圖7 MT1葉片表面等熵馬赫數及努賽爾數分布Fig.7 Isentropic Mach number and Nu of MT1 vane

圖8 示出了冷卻氣經氣膜孔與主流氣摻混的過程，冷卻氣進入主流通道后在葉片表面形成氣膜，有效地把高溫氣體與葉片隔離開，達到了保護葉片的目的。

圖8 冷卻氣與主流氣摻混跡線Fig.8 Mixing stream trace of cooling air and gas

4 結論

CU_Turbo是一個基于非結構網格的三維計算軟件，由于其對復雜內冷結構渦輪葉片的適應能力及氣熱耦合方法的實用性，使得軟件非常適用于渦輪葉片的流動分析和換熱機理研究，從而優化渦輪葉片的冷卻流路設計，有較大的工程應用價值。本文選擇對帶內冷的渦輪導葉MarKⅡ和帶6排氣膜孔的渦輪導葉MT1，進行流場和溫度場耦合計算，得出如下結論：

(1)CU_Turbo軟件能準確預測帶內冷的渦輪葉片MarkⅡ和MT1的流場及溫度場。MarkⅡ在超聲速工況下計算得到的表面溫度與試驗值符合良好，壁面溫度計算值最大誤差小于3%；MT1的換熱分布同樣與試驗值吻合較好。

(2)該軟件采用隱式時間推進，隱式通量Jaco?bians矩陣采用對Roe通量的一種近似方法求解，進一步在渦輪葉片計算中驗證了將Jacobians矩陣與Roe格式結合這種方法的實用性和可靠性。

(3)附面層內不同流動狀態(層流、湍流和轉捩)對當地的傳熱過程有很大影響，氣膜冷卻可顯著減少氣膜孔下游葉片與高溫燃氣的熱交換，有效保護葉片。

[1]馮國泰，黃家驊，李海濱，等.渦輪發動機三維多場耦合數值仿真試驗臺的數學模型研究[C]//.中國工程熱物理學會熱機氣動熱力學學術會議.2000.

[2]Venkatakrishman V.Convergence to Steady State Solu?tions of the Euler Equations on Unstructured Grids with Limiters[J].Journal of Computational Physics，1995，118 (1)：120—130.

[3]Luo Hong，Baum J D，L?hner R.A Fast,Matrix-Free Im?plicit Method for Computing Low Mach Number Flows on Unstructured Grids[R].AIAA 99-3315，1999.

[4]Sheng C，Wang X.Characteristic Variable Boundary Con?ditions for Arbitrary Mach Number Algorithm in Rotating Frame[R].AIAA 2003-3976，2003.

[5]Sheng Chun-hua，Xue Qing-luan.Aerothermal Analysis of Turbine Blades Using an Unstructured Flow Solv?er-U2NCLE[R].AIAA 2005-4683，2005.

[6]Wang X.A Preconditioning Algorithm for Turbomachinery Viscous Flow Simulation[D].MS：Mississippi State Univer?sity，2005.

[7]史萬里，葛寧，Sheng Chunhua.粘性計算中預處理的隱式算法[J].空氣動力學學報，2009，27(5)：565—571.

[8]王靖宇.渦輪內非結構網格氣動熱力耦合數值計算方法研究[D].南京：南京航空航天大學，2012.

[9]趙滋陽.非結構網格下氣動熱力耦合數值方法研究[D].南京：南京航空航天大學，2010.

[10]Hylton L D，Mihelc M S，Turner E R，et al.Analytical and Experimental Evaluation of the Heat Transfer Distribution overtheSurfacesofTurbineVanes[R].NASA CR-168015，1983.

[11]Bohn D，Heuer T.Conjugate Flow and Heat Transfer Cal?culation of a High Pressure Turbine Nozzle Guide Vane [R].AIAA 2001-3304，2001.

[12]Adami P，Martelli F，Chana K S，et al.Numerical Predic?tionsofFilmCooledNGVBlades[R].ASME GT2003-38861，2003.

Study of Turbine Conjugate Heat Transfer Simulation Based on Unstructured Grids

TANG Chao1，GE Ning1，SHENG Chun-hua2
(1.College of Energy and Power Engineering，NUAA，Nanjing 210016，China；2.College of Engineering，University of Toledo，OH 43606，USA)

This paper describes an attempt to predict the aerodynamic pressure loading and heat transfer of two cases with an unstructured Navier-Stokes flow solver,which include a NASA turbine vanes-MarkII and NGV MT1.A three-dimensional Reynolds-averaged Navier-Stokes unstructured flow solver,referred to as CU_Turbo,has been modified to solve the governing equations in the flow field and a thermal diffusion equa?tion inside the solid region in a coupled manner.During implicit time iterative,an approximate method of Roe flux is used for finding the Jacobians matrix.The results showed that the numerical simulation results were agreement well with the experimental values.The results verified the feasibility and effectiveness of the code.So the code provides a new method of calculation and analysis for the engineering design of tur?bine;flow condition in boundary layer and film cooling have great impacts on the local heat transfer in tur?bine vane channels.

turbine vane；conjugate heat transfer；film cooling；implicit time iterative；unstructured grids

V231.1

A

1672-2620(2013)02-0033-05

2012-08-06；

2012-12-22

航空推進技術驗證計劃項目

湯超(1988-)，男，安徽巢湖人，碩士研究生，主要從事葉輪機械氣動熱力學研究。