李心語,劉火星
(北京航空航天大學能源與動力工程學院,北京 100083)
渦輪內部的流場與溫度場是高度耦合的,一方面渦輪內部流場存在著如激波、二次渦系等復雜流動結構,會影響到部件溫度場分布及其換熱特性,另一方面,固體域的溫度分布及冷卻射流與主流的相互作用也將反過來影響流場的流動特征[1]。因此,理解和掌握渦輪內部流場溫度場相互作用機制,對于提高冷卻效率和準確預測熱負荷有著重要意義[2]。
隨著CFD技術的發展和整體計算能力的提升,流熱耦合(Conjugate Heat Transfer,CHT)計算方法被越來越多地應用于渦輪內部流動換熱機理研究及優化設計[3]。NASA的Heidmann等[4]開發了一套三維流熱耦合求解程序,通過流體模塊Glenn-HT求解流體域,使用邊界單元法求解固體域,該程序在真實氣冷渦輪算例上實現了多區耦合模擬。Bohn等[5]使用CHT-flow共軛傳熱程序(基于有限體積格式),對復雜氣冷葉片熱負荷進行了預測與耦合傳熱分析,結果可用于渦輪葉片冷卻系統的優化設計。國內李宇[6]、蘇欣榮[7]等也在有關方面做出了大量工作。隨著流熱耦合計算技術成熟度的提高,相關商用軟件也獲得了廣泛應用。Eifel等[8]使用CFX軟件對復雜冷卻渦輪開展了流熱耦合數值模擬,基于結果對幾何結構進行優化。Mazur等[9]使用STAR CD軟件對氣冷高壓渦輪導葉開展了數值模擬,結果可用于分析葉片強度損傷及壽命預測。此外,哈爾濱工業大學的董平[10]也針對C3X葉片前緣氣冷結構為研究對象開展流熱耦合數值模擬,較為精確地模擬出了冷卻出流過程。
為更深入了解渦輪葉片流熱耦合作用機制,探索流場溫度場匹配的流動換熱機理,本文以2種氣冷渦輪導葉為研究對象,開展了流熱耦合數值模擬。內冷渦輪導葉算例中,對影響耦合計算精度的因素進行定性定量的分析討論,利用修正后的計算模型對該導葉多場特性及耦合機理展開研究。在此基礎上,以帶有氣膜冷卻孔及內冷通道的氣冷渦輪導葉為研究對象,開展葉柵內部流動換熱特性分析的同時,重點圍繞冷卻射流與主流的相互作用,討論近壁邊界層中流熱耦合關系及氣冷效率影響因素等相關問題。
本文選取研究對象為MARK Ⅱ高壓渦輪導葉,NASA于1983年對其流動換熱特性進行了大量實驗研究[11]。MARKⅡ整體為直列葉柵結構,葉身配有10個徑向冷卻孔,尾緣呈鈍角,數值模擬建模參數與實驗文獻中的幾何參數均保持一致。葉片氣動模型及內部完整冷卻結構如圖1所示。

圖1 MARKⅡ幾何氣動外形及內冷結構編號示意圖Fig.1 MARKⅡgeometry and schematic diagram of internal cooling structure numbering
網格采用ANSYS/ICEM 軟件劃分,為保證正交性,本算例使用混合網格。其中結構較復雜的葉片固體域采用非結構化網格,網格數為67萬。流體域采用“H-O-H”結構化網格,設置葉片表面第1層網格高度為1×10-6m以確保壁面y+<1,最終網格數為102萬。MARKⅡ葉片三維網格示意圖如圖2所示。

圖2 MARKⅡ葉片三維網格劃分情況Fig.2 Generation of 3D mesh of MarkⅡblade
本文采用實驗4411號工況進行數值模擬,該工況進氣方向為軸向,主流入口總溫為784 K,進口總壓為342255 Pa,湍流度設置為6.5%,出口背壓為204560 Pa。10個徑向冷卻通道根據經驗公式換算實驗點完成第三類邊界條件(對流換熱系數及參考溫度)設置。由經驗公式可推導出內冷通道對流換熱系數h表達式為

式中:h為對流換熱系數,利用努賽爾數Nu及流體熱傳導系數κ的積除以冷卻通道直徑D得出,Nu可由修正系數Cr、普朗特數Pr及雷諾數ReD(與冷卻通道直徑D及冷卻氣體出口速度有關)的關系推導。
計算中,主流工質設為理想可壓燃氣,分子黏性和熱傳導系數均采用Sutherland公式擬合為與溫度相關的函數。葉片材料設為ASTM 標準310不銹鋼,其密度ρ為8030 kg/m3,比熱Cp為502 J/(kg·K)。固體熱傳導系數κ與溫度T相關的表達式[12]為

發展高可信度數值技術有賴于對影響流動和換熱預測精度的因素進行系統細致的分析[13]。算例使用ANSYS/CFX軟件,使用RANS方程求解流體域,能量方程求解固體域,計算差分格式均采用二階迎風格式。為提高整體流熱耦合數值精度,本節對3種影響精度的因素進行定性定量探討。
圖3中,橫坐標為軸向長度x與葉片弦長L之比,縱坐標為葉表溫度T與總溫T*之比。可知,采取不同壁溫設置的算例在溫度預測方面表現差異較大。其中,絕熱壁面的邊界條件無法體現冷卻作用,故結果偏離實驗值較遠。而基于流熱耦合的計算結果最為貼近真實溫度分布,從而證明了流熱耦合計算方法的必要性。

圖3 不同壁溫條件下葉表溫度分布曲線Fig.3 Surface temperature distribution curves of blade under different wall temperature conditions
不同湍流轉捩模型的求解算例對比如圖4所示。從溫度分布曲線可以看出,由于全湍流模型高估了局部換熱能力,各模型計算結果相差逾60 K。其中,SST-γReθ模型最大誤差小于3%,平均誤差小于1%。綜合來看,SST-γReθ模型較其余模型更符合真實值,因此將采用該模型進行后續計算。

圖4 不同計算模型求解的葉表溫度分布曲線Fig.4 Surface temperature distribution curves of bladesolved by different calculation models
為進一步提高精度,圖5比較了湍流普朗特數PrT為定值和Kays-Crawford表達式[14]求解PrT的多個算例。Kays-Crawford公式列舉如下:

圖5 不同Pr T 條件求解的葉表溫度分布曲線Fig.5 Surface temperature distribution curves of blade solved under different Pr T conditions

式中:μT和μL分別為湍流相關和層流相關的動力黏性系數。可以發現,各能量封閉方法預測的流動換熱數值結果差異較小。結合文獻[15]可得出結論,盡管湍流普朗特數在黏性底層數值變化較大,但底層的能量輸運作用主要受導熱系數和動力黏度影響,PrT對其影響較小。由此推知,選取能量封閉方法時,將湍流普朗特數設為合理常數即可。
結合以上提高精度有關結論,本文的數值模擬算例將均使用流熱耦合方法進行計算,選取SST-γReθ湍流轉捩模型求解流體域方程,湍流普朗特數設為0.7。
采用1.2節所述修正模型對4411號工況開展了流熱耦合計算。結合圖6給出的4411號工況下壓力分布(縱坐標為壓力P與總壓P*之比)及圖7中的馬赫數云圖,分析葉柵內部流動特性。吸力面上,隨著流道面積減小,燃氣開始加速,軸向弦長比0% ~40%位置處馬赫數Ma迅速上升,壓力迅速降低;而至軸向弦長比41%處附近,燃氣達到Ma為1.5的超聲速引發激波,使氣流大幅降速并導致壓力陡然上升。結合云圖來看,該激波與邊界層進行相互作用從而誘發轉捩,流場出現分離泡,流線分布參見圖7中局部放大圖。至軸向弦長比58%后,分離泡消失,近壁氣流再附后吸力面流場已處于完全發展的全湍流階段。另外,在MARKⅡ葉表壓力面處,燃氣自前緣駐點開始緩慢加速,葉柵型面壓力分布呈下降趨勢。需要注意,在壓力面軸向弦長比為95%處會形成一道馬赫數Ma為1.1的弱激波,該激波處壓力梯度達到最大,壓力分布處于最低峰值。結合折線圖發現,經過弱激波后壓力出現回升。

圖6 葉表壓力分布曲線Fig.6 Surface pressure distribution curve of blade

圖7 MARKⅡ葉表馬赫數云圖Fig.7 Contour of Mach number at surface of MARKⅡblade
內冷葉柵流場與溫度場存在一定的相互作用,因此,結合上述流場分布開展換熱特性分析。經流熱耦合計算得到葉表溫度分布結果,如圖8所示。葉片吸力面前緣氣流開始加速,氣體動能增加內能減少,從而吸力面前緣氣體溫度呈下降趨勢;至軸向弦長比41%位置處出現激波,與邊界層相互作用產生的分離泡使動能損失而溫度場內能增加,軸向弦長比41% ~58%位置上溫度分布呈上升趨勢。結合換熱系數分布(見圖9)分析可知,在軸向弦長比58%位置附近,葉片的冷卻結構——內冷通道起到有效作用,換熱系數h與總換熱系數h0在該處溫度大幅下降。至吸力面軸向弦長比80%處,流場云圖顯示此時湍流已達充分發展階段,強化了葉片型面的對流換熱效果,溫度分布呈回升趨勢。

圖8 葉表溫度分布曲線Fig.8 Surface temperature distribution curve of blade

圖9 葉表換熱系數分布曲線Fig.9 Heat transfer coefficient distribution curve of blade surface
結合溫度分布云圖(見圖10)對葉表壓力面一側進行分析,高溫駐點后流場氣流流速較低,內能轉化為動能的效率較低。而打在前緣處較為集中的3個徑向冷卻孔起到有效作用,軸向弦長比0% ~20%處溫度呈下降趨勢。在下游葉表換熱系數與溫度分布出現波動,隨著軸向弦長比的增大與冷卻孔的尺寸減小,氣冷換熱效果隨著冷氣流量的減小而變弱,至尾緣葉表溫度與來流總溫之比已達0.78。由此推知,葉片內冷通道的對流換熱冷卻效果與其分布位置、尺寸大小、冷氣流量分配有較大關聯,導葉設計時應著重前緣尾緣處的冷卻效果。

圖10 流體域及葉片內部溫度分布云圖Fig.10 Temperature distribution contour in fluid domain and blade interior
利用流熱耦合計算結果,對4411號工況下MARKⅡ葉片開展有限元分析。圖11和圖12給出了4411號工況下葉片的等效應力分布和給定機匣輪轂約束后的總變形量。可見,葉片的前緣尾緣熱負荷均偏高,引起變形較為劇烈。為確保渦輪葉片強度及使用壽命,設計冷卻結構時應予以重點關注。

圖11 葉片內部等效熱應力分布Fig.11 Distribution of equivalent thermal stress inside blade

圖12 給定約束后葉片總變形量Fig.12 Total variation quantity of blade with given constraints
前緣氣膜冷卻在渦輪導葉的熱防護中起到了重要作用,而氣膜孔冷卻射流與高溫主流間的摻混過程較為復雜,涉及到多場相互影響及邊界層對流換熱特性。基于第1節內冷渦輪流熱耦合計算方法及結論,將對氣冷渦輪導葉開展多場耦合計算。
本節所研究的氣冷渦輪葉片改型自MARKⅡ高壓渦輪導葉,葉片氣動造型與第1節保持一致。為節省計算資源,葉片展向高度變為原高度76.2 mm的四分之一。為減少氣動損失,氣膜冷卻孔打于葉片壓力面弦長12%處(前緣附近),葉片展向高度的50%處。由于本節涉及渦輪葉片上氣膜冷卻孔射流的基礎性機理研究,設計冷卻射流角度時,主要考慮不同吹風比下冷卻射流與橫向高溫主流相互作用、孔附近流場結構的變化及其對流熱耦合機理的影響,而非尋求全場冷卻最優解。因此,從機理性研究角度出發,為探尋冷卻射流與主流間更直接的相互影響,不同于工程常用的射流傾角,本節圓形氣膜單孔設為垂直出流。角度參考文獻[16],如圖13所示,定義α為展向夾角,β為流向夾角,該氣膜孔直徑為2 mm,出氣方向α和β角均為90°。此外在保留其余9個原始徑向冷卻通道前提下,將2號徑向冷卻通道改為矩形對該氣膜冷卻孔供氣,整體葉片造型如圖14所示。

圖13 氣膜冷卻孔出氣角度示意圖Fig.13 Schematic diagram of outlet angle of film-cooling hole

圖14 氣冷葉片幾何造型Fig.14 Geometry model of air-cooled blade
網格劃分思路同第1節,利用混合網格技術劃分改型MARK Ⅱ葉片。葉片燃氣域采用分塊結構化H-O-H網格,壁面第1層網格高度符合計算模型要求。氣冷流體部分劃分蝶形網格,根據圖13定義方向對其進行加密,詳見圖15。由于固體域并不涉及冷卻射流作用及邊界層擾動,其網格精度相對流體域網格要求較低,固體域采用非結構化網格,組裝后形成最終計算網格,組裝示意圖如圖16所示。整體網格數為85萬,計算時將3塊網格區域導入CFX中組裝計算。

圖15 流體域網格劃分Fig.15 Mesh generation of fluid domain

圖16 葉片網格組裝示意圖Fig.16 Schematic diagram of assembled mesh of blade
本節算例模擬工況依然為4411號工況,主流域、冷氣域及葉片材料分別選為理想可壓燃氣、理想氣體和ASTM 標準310不銹鋼,各材料相關物理量及邊界條件參見第1節。由于本節著重氣膜冷卻孔附近微小結構的流熱耦合機理,多個算例給定冷氣總溫不變且射速流量跨度較大,以吹風比表示效果可能不甚明顯。因此,為更直觀表現單孔冷氣射流對整體流場溫度場的影響,基本算例將根據動量比及吹風比給定射流進口速度10 m/s,入口溫度300 K。后續算例將同樣根據吹風比等無量綱系數給定射速1~50 m/s的進口邊界條件,總溫均為300 K。
根據第1節結論,使用ANSYS/CFX軟件對改型MARKⅡ葉片開展基于SST-γReθ湍流轉捩模型的流熱耦合數值計算。計算中,每50步保存一次結果以便進行觀察。
分析馬赫數云圖(見圖17)得知,氣膜冷卻孔對整體流場影響不大,吸力面激波強度及分離位置幾乎并未改變,壓力面馬赫數分布也變化不大。而對比溫度分布云圖(見圖18)發現,與第1節結果差異較大,葉表壓力面前緣處高溫燃氣速度梯度較小,溫度呈緩慢下降趨勢。至軸向弦長比為20%處,冷卻氣體開始發揮作用,從圖18中的放大圖可以看到,冷卻出流穿透邊界層后被自由流壓回壁面,在葉表形成了一層低溫冷卻氣膜。而隨著軸向弦長比的增加,葉表溫度分布呈上升趨勢,這是由于冷氣與燃氣的摻混作用越來越劇烈,氣膜保護范圍逐漸減小。

圖17 改型MARKⅡ葉表馬赫數云圖Fig.17 Contour of Mach number at surface of retrofitted MARKⅡblade

圖18 葉表溫度云圖Fig.18 Contour of temperature at surface of blade
由以上分析可推知,氣膜冷卻的冷卻效果與其射流強度有重要關聯。為更深入研究冷卻射流與主流相互作用機制,在同等條件下設置氣膜冷卻孔射速為1~30 m/s的6個算例,探尋不同吹風比下摻混作用對氣膜冷卻效率的影響規律。不同吹風比條件下,算例壓力面溫度分布對比如圖19所示,Vi代表氣膜冷卻孔射速。結合流線放大圖20(a)可知,在冷氣初始位置,冷氣射速在一定范圍內越大,則熱防護效果越好。這是由于冷氣射速越小越難穿透主流邊界層,隨著摻混作用的發生,其冷卻區域變短,氣膜層無法高效保護下游葉片壁面,此時溫度分布曲線顯示冷卻效果與冷氣射速存在正相關關系。而當冷氣射速達到一定值,如圖20(b)所示,冷氣在壓力面主流中擾流加劇,導致燃氣在射流位置下游產生高溫區。隨后在下游自由流重新將冷氣壓回壁面,氣膜層開始發揮作用降低葉表溫度,且冷氣流量的冷卻效率隨著射速增加而變高。增加射速至30 m/s高速時,隨著氣膜層向下游發展,大量冷卻氣體發揮其優勢,從圖20(c)中可以明顯看到,氣膜層對葉片壁面保護較好,溫度隨下游距射流孔距離的增加而降低。但結合圖19還需注意,大吹風比高速射流條件下,冷氣噴射可能出現較輕微的冷氣吹離現象,為達到最高冷卻效率,后續氣膜圓孔的設計中應合理進行兩者的權衡。

圖19 不同吹風比下壓力面溫度分布Fig.19 Temperature distribution of pressure surface under different blowing ratios

圖20 不同冷氣射速下氣膜冷卻孔下游流線放大圖Fig.20 Enlarged streamline distribution of downstream hole at different flow velocities of cold air
由此得出結論,想要加強氣膜層的熱防護,不能一味提高冷氣射流強度。射速較低時,冷氣無法穿透燃氣附面層,增加冷氣流量可提高氣膜層冷卻效率。而當射速較高使冷氣可以突破附面層時,主流繞過射流產生高溫區后,氣膜冷卻層開始發揮作用。因此,高速的冷卻氣體在局部高溫區上游冷卻效果較差,而在下游冷卻效率較高。
真實渦輪工作條件下,流場和溫度場相互耦合。一方面,渦輪內部流動特性會對部件熱負荷產生顯著影響,另一方面,氣流在換熱邊界的相互作用也將影響流動圖畫。研究表明,氣流相互作用使下游流場出現復雜渦系[17]。本節在2種高速射流下取冷氣孔下游2 mm、6 mm、10 mm、14 mm位置作為橫向流線提取點,流線分布圖如圖21所示。可以看到,流線在4個位置上捕捉到了流場中渦系結構的生成與消失。從圖21中可以看出,隨著與氣膜冷卻孔距離增加,主流的剪切作用使冷氣方向發生改變,流場中出現成對反向的腎形渦,其強度隨著剪切效應的變化而經歷先遞增后遞減的過程。此外,由于氣膜冷卻孔打在壓力面,冷卻射流受到主流速度分量的影響導致渦核向上方偏移。對比流線分布得知,射速越快,則腎形渦對的強度就越高且形狀更飽滿。需要注意,在射速為50 m/s的大吹風比算例中,冷氣與高溫主流摻混劇烈,繞進射流下游區域的少量主流存在貼近壁面的運動趨勢。這種向下的速度分量受腎形渦對壁面附近氣流的強剪切作用,將導致壁面附近的流場出現次生小型二次渦對,如放大圖22所示。隨著距離的增加,該渦系結構將隨著主渦系強度的降低而消失。由于該次生對轉渦尺寸較小并很快消失,從流線及溫度分布來看,并未對流動結構及換熱特性造成一定影響。

圖21 高速射流條件下不同位置橫向流線發展過程Fig.21 Development of transverse streamlines at different positions under high-speed jet conditions

圖22 高速射流條件下受剪切作用形成的二次渦對Fig.22 Secondary vortex pair induced by shear under high-speed jet conditions
觀察不同氣冷出流速度下壁面極限流線圖(見圖23)得知,隨著冷氣流量的逐漸增大,流場中將出現更明確的射流擾流圖譜。氣膜冷卻孔上游的主流燃氣由于受到冷氣阻礙將形成一個鞍點,該鞍點分離出2條分離線構成馬蹄渦。射速越高,則擾流效應越強,氣流間的相互阻礙作用也將越大,馬蹄渦形狀也更為飽滿。而在氣膜冷卻孔下游,冷卻出流與主流的剪切效應導致了第二個鞍點的出現,其附近流線形成前文所描述的腎形渦對。隨著氣膜邊界層的發展,馬蹄渦分離線及腎形渦分離線在壓力梯度作用下完成再附,流線方向與中截線趨于一致。綜上所述,冷卻出流與主流的相互作用產生了復雜的流動結構,而這些流動結構(如腎形渦、馬蹄渦等)也將影響溫度場的分布,如前文所述氣流的摻混作用會導致氣膜冷卻孔附近形成局部高溫區。

圖23 高速射流條件下局部壁面極限流線Fig.23 Local wall limit streamlines under high-speed jet conditions
1)采用流熱耦合計算方法和合適的湍流轉捩模型有利于提高氣冷渦輪導葉換熱特性的計算精度,湍流普朗特數取值影響不大。
2)經流熱耦合計算方法分析發現,導葉前緣尾緣附近溫度及應力水平較高,應著重上述2個位置進行合理的冷卻結構設計。
3)氣冷渦輪導葉流場中,冷氣與主流相互作用產生的復雜流動結構對溫度場具有較大影響,如會在氣膜冷卻孔附近產生局部高溫區。
4)氣膜冷卻孔冷氣流量與葉片冷卻效率并非簡單正比關系。射速低時,增加冷氣流量可提高氣膜層冷卻效率。增加冷氣流量到一定值時,冷氣流量增加將導致氣膜冷卻孔后上游冷卻效果變差,下游冷卻效果變好。
5)射速較高條件下,冷卻出流因主流的剪切效應形成反向腎形渦對,主流也會受到冷氣擾流影響而在氣膜冷卻孔前產生強度較小的馬蹄渦,對導葉溫度分布產生影響。