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

湯超，葛寧，SHENG Chun-hua

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

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

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

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

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

1 引言

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

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

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

2 數(shù)值計(jì)算方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

向前推進(jìn)

向后推進(jìn)

求得的修正項(xiàng)Δq更新本時(shí)間步：

2.4通量Jacobians矩陣的計(jì)算方法

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

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

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

3 算例驗(yàn)證及分析

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

3.1MarkII

MarkII葉片是有10個(gè)冷卻孔的渦輪直葉片，葉高76.20 mm，弦長(zhǎng)136.22 mm。圖2給出了葉片及計(jì)算域幾何尺寸，其中左、右表面分別為進(jìn)口和出口，上、下兩邊界面為周期性邊界面，中間為帶冷卻孔葉片，頂端和底端面為對(duì)稱邊界。Hylton已在接近工程運(yùn)行條件下做過(guò)大量試驗(yàn)[10]，本文模擬的是Hyl?ton試驗(yàn)中的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的方法劃分計(jì)算域網(wǎng)格，在貼近葉片表面的壁面處使用邊界層網(wǎng)格，使得第一層y+≤1以保證精確模擬邊界層內(nèi)流動(dòng)。采用商用軟件Gambit生成網(wǎng)格，通過(guò)CU_Turbo接口軟件GMT[9]實(shí)現(xiàn)數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換。

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

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

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

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

圖5 給出了葉片表面中徑處的性能曲線。從圖5(a)中可知，在吸力面，從滯止點(diǎn)開(kāi)始?jí)毫ρ杆俳档?，在x cx=0.45處壓力減小到最小值，激波與壓力分布的位置相對(duì)應(yīng)。

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

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

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

3.2MT1

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

圖6所示為葉片計(jì)算域及網(wǎng)格，左側(cè)為進(jìn)口邊界，右側(cè)為出口邊界，計(jì)算域上、下端設(shè)為對(duì)稱邊界，前、后為周期邊界，所有壁面均設(shè)為定溫邊界。各邊界具體參數(shù)為：主流來(lái)流總壓P*=4.60×105Pa，主流來(lái)流總溫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葉片計(jì)算域及網(wǎng)格Fig.6 Computation domain and grids of MT1 vane

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

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

圖7(b)中，努賽爾數(shù)的定義式為Nu=hd λ，d為特征長(zhǎng)度，λ為流體域熱傳導(dǎo)系數(shù)。對(duì)比有無(wú)氣膜孔時(shí)的努賽爾數(shù)分布，可清楚看出氣膜孔對(duì)葉片表面換熱的影響，在每組氣膜孔后換熱強(qiáng)度均有較大幅度的降低，有效減小了葉片表面的熱負(fù)荷。對(duì)比帶氣膜的計(jì)算值與試驗(yàn)值，壓力面由于計(jì)算得出的溫度附面層偏厚導(dǎo)致努賽爾數(shù)的計(jì)算值稍偏低，吸力面峰值處則略高于試驗(yàn)值，總體趨勢(shì)吻合良好。同時(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)氣膜孔與主流氣摻混的過(guò)程，冷卻氣進(jìn)入主流通道后在葉片表面形成氣膜，有效地把高溫氣體與葉片隔離開(kāi)，達(dá)到了保護(hù)葉片的目的。

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

4 結(jié)論

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

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

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

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

[1]馮國(guó)泰，黃家驊，李海濱，等.渦輪發(fā)動(dòng)機(jī)三維多場(chǎng)耦合數(shù)值仿真試驗(yàn)臺(tái)的數(shù)學(xué)模型研究[C]//.中國(guó)工程熱物理學(xué)會(huì)熱機(jī)氣動(dòng)熱力學(xué)學(xué)術(shù)會(huì)議.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]史萬(wàn)里，葛寧，Sheng Chunhua.粘性計(jì)算中預(yù)處理的隱式算法[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27(5)：565—571.

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

[9]趙滋陽(yáng).非結(jié)構(gòu)網(wǎng)格下氣動(dò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ù)驗(yàn)證計(jì)劃項(xiàng)目

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