湯超，葛寧，SHENG Chun-hua

基于非結(jié)構(gòu)網(wǎng)格的渦輪氣熱耦合數(shù)值方法研究

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

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

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

渦輪葉片；氣熱耦合；氣膜冷卻；隱式時間推進(jìn)；非結(jié)構(gòu)網(wǎng)格

1 引言

渦輪是發(fā)動機中熱負(fù)荷和動力負(fù)荷最大的部件，采用有效的冷卻措施是發(fā)動機安全可靠工作的有力保證，也是降低材料成本的有效措施。近年來，渦輪前燃?xì)鉁囟鹊闹鹉晏岣吲c渦輪的有效冷卻，特別是葉片內(nèi)部空氣冷卻技術(shù)的迅速提高分不開。有研究表明，渦輪葉片設(shè)計中，如果預(yù)測的葉片溫度超過實際工作溫度28℃，葉片壽命將減半[1]。因此準(zhǔn)確估算葉片的傳熱系數(shù)和溫度，有助于預(yù)防熱腐蝕和設(shè)計高效冷卻系統(tǒng)。

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

本文使用課題組自主開發(fā)的基于非結(jié)構(gòu)網(wǎng)格的三維計算軟件CU_Turbo，通過對葉型MarkⅡ和MT1進(jìn)行數(shù)值模擬，驗證了氣熱耦合計算方法的實用性和有效性，為渦輪工程設(shè)計提供了一種新的計算分析方法。

2 數(shù)值計算方法

2.1控制方程和數(shù)值方法

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

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

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

2.2非結(jié)構(gòu)網(wǎng)格的氣熱耦合方法

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

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

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

2.3隱式推進(jìn)數(shù)值方法

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

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

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

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

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

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

向前推進(jìn)

向后推進(jìn)

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

2.4通量Jacobians矩陣的計算方法

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

式中：帶上標(biāo)⌒的項采用修正的Roe平均變量計算，認(rèn)為該項為當(dāng)?shù)刂挡⒓僭O(shè)為常數(shù)。求微分后可得任一控制面對于界面左右兩側(cè)的對流通量Jacobians矩陣的表達(dá)式：

粘性通量Jacobians矩陣的求解采用薄層假設(shè)，求解過程參見文獻(xiàn)[8]。

3 算例驗證及分析

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

3.1MarkII

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

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

采用cooper的方法劃分計算域網(wǎng)格，在貼近葉片表面的壁面處使用邊界層網(wǎng)格，使得第一層y+≤1以保證精確模擬邊界層內(nèi)流動。采用商用軟件Gambit生成網(wǎng)格，通過CU_Turbo接口軟件GMT[9]實現(xiàn)數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換。

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

圖3馬赫數(shù)等值線分布Fig.3 Mach number contours

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

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

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

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

圖5(c)中，換熱系數(shù)h定義為h=qw(T0-Tw)，qw為葉片表面熱流密度，T0為氣體來流溫度?？梢?，溫度影響換熱，換熱系數(shù)與表面溫度分布呈現(xiàn)出很強的相關(guān)性，計算值與試驗值趨勢基本吻合。駐點后壓力面和吸力面的熱傳導(dǎo)系數(shù)急劇下降，吸力面激波處誘導(dǎo)分離轉(zhuǎn)捩，導(dǎo)致熱傳導(dǎo)系數(shù)開始迅速上升。由于計算的溫度邊界層比實際的厚，使得壓力面換熱系數(shù)的計算值比試驗值偏低，靠近尾緣區(qū)域的誤差，還可能是由于當(dāng)?shù)馗矫鎸恿鲃邮芟噜徣~片壓力面強激波影響所致。

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

3.2MT1

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

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

圖6 MT1葉片計算域及網(wǎng)格Fig.6 Computation domain and grids of MT1 vane

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

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

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

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

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

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

4 結(jié)論

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

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

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

(3)附面層內(nèi)不同流動狀態(tài)(層流、湍流和轉(zhuǎn)捩)對當(dāng)?shù)氐膫鳠徇^程有很大影響，氣膜冷卻可顯著減少氣膜孔下游葉片與高溫燃?xì)獾臒峤粨Q，有效保護(hù)葉片。

[1]馮國泰，黃家驊，李海濱，等.渦輪發(fā)動機三維多場耦合數(shù)值仿真試驗臺的數(shù)學(xué)模型研究[C]//.中國工程熱物理學(xué)會熱機氣動熱力學(xué)學(xué)術(shù)會議.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.粘性計算中預(yù)處理的隱式算法[J].空氣動力學(xué)學(xué)報，2009，27(5)：565—571.

[8]王靖宇.渦輪內(nèi)非結(jié)構(gòu)網(wǎng)格氣動熱力耦合數(shù)值計算方法研究[D].南京：南京航空航天大學(xué)，2012.

[9]趙滋陽.非結(jié)構(gòu)網(wǎng)格下氣動熱力耦合數(shù)值方法研究[D].南京：南京航空航天大學(xué)，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

航空推進(jìn)技術(shù)驗證計劃項目

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