不同燃氣輸運特性計算方法對渦輪平面葉柵對流換熱的影響

薛 鈺,劉景源

(南昌航空大學飛行器工程學院,南昌 330063)

以燃料與氧化劑化學放熱產生的熱能轉變為機械功為動力的推進系統一直是研究的熱點。研究性能高、質量輕、尺寸小的航空發動機一直是世界各國追求的目標。提高渦輪前燃氣溫度是達到上述目標的重要因素。而提高渦輪前燃氣溫度的貢獻中,約70%是來自換熱和冷卻技術的發展[1]。因此航空發動機推進系統內的對流換熱問題是研究的重點之一。

文獻[1]針對燃氣渦輪換熱,系統地論述了換熱分析及冷卻等技術及研究進展。文獻[2]對小型航空發動機的換熱問題進行了綜述,指出了燃燒室、渦輪葉片內外表面等獨特的換熱問題。文獻[3]則對一種自激Helmholtz 脈沖燃燒室的彎尾管的換熱問題進行了實驗及數值模擬研究,為此類發動機燃燒室設計提供技術支撐。文獻[4]應用實驗及理論分析研究了燃燒室內在流動方向加裝波紋管以提高渦輪發動機燃燒室在高溫燃氣熱沖擊下的使用期限。由于在波紋管壁面開孔導致了燃氣和冷卻空氣氣膜摻混,因此文獻[4]研究摻混狀態下冷卻效果。關于渦輪葉片對流換熱問題,氣膜冷卻[5-6]及多場耦合[7]是研究的熱點。

上述研究的換熱高溫氣體絕大多數采用高溫空氣代替高溫燃氣,但是燃燒室內產生的高溫氣體為多組分的燃氣,無論高溫燃氣與航空發動機部件壁面的對流換熱還是冷卻空氣和高溫燃氣的摻混,均涉及多組分氣體的輸運過程。多組分的氣體動量輸運導致了黏性的產生,而溫度的差別導致換熱。高溫燃氣和高溫空氣輸運特性不同以及輸運特性計算方法的不同均引起換熱的差別。

為了計及多組分燃氣輸運特性的影響,需要計算燃氣的輸運系數。而燃氣的輸運系數由多組分燃氣中每種氣體的質量分數和每種氣體的輸運特性計算方法決定。每種組分的質量分數是燃燒室燃料成分、燃燒室出口溫度等的函數。并且組分的輸運特性為燃燒室出口溫度的函數。另外,組分的輸運特性的計算方法也有多種。因此研究燃氣的輸運特性計算方法對渦輪換熱的影響,需要考慮燃燒室出口溫度不同導致的燃氣組分質量分數、每種組分不同輸運特性計算方法等對對流換熱的影響。發動機不同工作狀態下,燃燒室噴入的燃油量不同，燃燒室出口溫度、燃氣質量分數亦不同。綜上,上述不同因素的組合數巨大,應提出一種快速計算不同燃氣輸運特性計算方法對渦輪換熱影響的方法,為工程應用提供參考。

隨著數值模擬及計算機技術的不斷進步,應進一步深入研究渦輪葉片對流換熱問題。因此，在計及燃氣多組分換熱的基礎上,提出一種不同燃氣組分輸運特性計算方法對對流換熱影響的快速算法,分析了燃氣各組分的不同輸運特性計算方法之間的不同組合對多組分燃氣對流換熱強度的影響,并以平板及渦輪平面葉柵為研究對象,應用數值模擬及理論分析方法,研究了高溫燃氣動量輸運及熱量輸運計算方法之間不同組合的差異對對流換熱強度的影響范圍。提出的多組分燃氣輸運特性快速計算方法可為燃氣渦輪的對流換熱問題提供支撐。

1 燃氣輸運特性計算

所用燃料的化學方程式為肖保國等[8]對RP-3航空燃油的化學動力學模型研究簡化后的C9.74H20.52。假定燃料在燃燒內與空氣完全燃燒,并且燃燒后燃氣組分僅為O2、N2、CO2及H2O等4種。以下僅對上述4種組分進行計算及分析。

1.1 燃氣輸運特性計算方法

單組分氣體分子黏性系數和熱導率計算公式有多種,但常用的有簡單的冪次律及Sutherland公式[9]。

分子黏性系數μ的冪次律及Sutherland公式為

(2)

熱導率k的冪次律及Sutherland公式為

(4)

式中:μ0、k0分別為常溫度T0下的分子黏性系數及熱導率；S、n為常數,皆取自文獻[9]。

由于燃氣為多組分混合氣體,可應用Wilke[10]建議的公式,計算燃氣的分子黏性系數μmix和熱導率kmix。

(6)

式中：μi、ki分別為各組分的分子黏性系數和熱導率,由式(1)～式(4)計算可得；xi為燃氣各組分的摩爾分數，可由文獻[11]中的方法計算給出；φij為配分函數；ns為混合氣體組分數目。

配分函數的表達式[10]為

(7)

式(7)中：μi、Mi分別為第i個組分的分子黏性系數及摩爾質量。

只需把式(7)中μi的換成ki即為熱導率k的配分函數。

所用的第3種燃氣不同組分的分子黏性系數及熱導率的計算方法由文獻[12]給出。

最后一種計算方法由國際工質物性計算軟件REFPROP得到[13]。

1.2 燃氣輸運特性計算結果

圖1所示為不同計算方法時燃氣4種組分的分子黏性系數μ隨溫度變化曲線。其中μN為文獻[12]給出的值,μP為式(1) 的計算值,μS為式(2) 的計算值,μR是軟件REFPROP的計算值(其中只給出了2 000 K的H2O的值)[13]。由圖1可見,燃氣4種組分的μ皆隨溫度的升高單調增加,且不同方法給出的各組分μ的差別也隨溫度的升高而逐漸增大。

圖1 不同計算方法時燃氣各組分μ隨溫度的變化Fig.1 Different methods of calculating viscosity coefficient of gas components at different outlet temperatures

圖2所示為不同方法給出的燃氣4種組分的熱導率k隨溫度變化曲線。其中各熱導率k的含義與圖1類似。圖中燃氣4種組分的k皆隨溫度的升高單調增加,且不同方法給出的k的差別也隨溫度的升高而逐漸變大。

圖2 不同方法時燃氣各組分k隨溫度的變化Fig.2 Different methods for calculating thermal conductivity of gas components at different temperatures

綜上,不同計算方法得到的μ和k不同,且之間的差別隨來流溫度的升高而升高。k不同計算方法之間的差別比μ的大。

1.3 燃氣各組分的普朗特數

對于燃氣多組分的對流換熱問題,動量的運輸及熱量的傳遞對對流換熱的影響可以用無量綱的Reynolds數及普朗特數表征。

混合燃氣的普朗特數為

Pr=μmixcPmix/kmix

(8)

式(8)中：燃氣定壓比熱cPmix由燃氣各組分的cPi和質量分數xi加權計算；各組分的cPi則用溫度的四次多項式給出。

圖3所示為不同計算方法給出的燃氣各組分的Pr隨來流溫度變化曲線。其中PrN為文獻[12]中各組分的Pr；PrR為軟件REFPROP計算的值(只給到2 000 K的H2O的值)；而μP-kP為式(1)和式(3)計算μ及k后得到的(下標P及S分別為用冪次律及Sutherland公式給出的μ及k)。由于Pr數是μ和k的函數,而不同計算方法得到的μ和k不同,因此不同計算方法組合后得出的Pr數有一定的差別。

圖3 不同計算方法燃氣各組分Pr隨來流溫度的變化Fig.3 Different methods for calculating Pr of gas components at different temperatures

1.4 不同燃氣輸運特性計算方法對對流換熱影響的快速算法

雖然燃氣中每種組分的輸運特性計算有多種組合,由于燃氣當地努賽爾數主要由其雷諾數及普朗特數決定(組分擴散系數影響很小),因而給出不同組合下燃氣普朗特數隨溫度變化所有曲線中的上下邊界曲線后,即可分析不同燃氣輸運特性計算方法的不同組合對對流換熱的影響程度。

由式(5)～式(7)易證,μmix隨μi的增大單調遞增,kmix隨ki的增大單調遞增。另外,從物理上也易得,當每種燃氣組分的質量分數(不同來流溫度下組分的質量分數可由燃燒室熱力計算給出,下同)及燃氣溫度均給定下,μmix隨μi的增大單調遞增,kmix隨ki的增大單調遞增。

由于μmix隨μi的增大單調遞增,kmix隨ki的增大單調遞增,因此,由式(8)燃氣Pr的上邊界由各組分μ的計算方法中最大值和各組分k的計算方法中最小值代入式(5)和式(6)計算得出；燃氣Pr的下邊界由各組分μ的計算方法中最小值和各組分k的計算方法中最大值代入式(5)和式(6)計算得出。圖4為燃氣Pr的上下邊界隨來流溫度變化曲線。由圖可見,燃氣Pr上下邊界之間的差別隨著燃氣溫度的升高逐漸增大。

圖4 燃氣Pr的上下邊界隨來流溫度變化曲線Fig.4 Maximum and minimum values of Pr of the gas at different temperatures

2 平板對流換熱數值模擬

2.1 幾何模型及邊界條件

平板幾何模型如圖5所示。整個計算區域300 mm×60 mm,其中板長250 mm,前緣點到來流邊界為50 mm。

圖5 平板幾何模型及邊界條件Fig.5 The geometry and boundary conditions of flat plate

數值計算采用不可壓N-S方程組及燃氣多組分輸運方程組。平板繞流的雷諾數為25 000,數值模擬假設層流流動。采用速度入口、壓強出口及壁面溫度為1 000 K的等溫、無滑移壁邊界條件。來流為多組分燃氣,溫度分別為1 200、2 200 K(不同來流溫度下組分的質量分數可由燃燒室熱力計算給出,下同)；平板前緣到來流邊界為對稱條件。

2.2 計算結果及分析

2.2.1 網格無關性驗證

采用網格點數為9×104、17×104、34×104的3套網格,對來流溫度為1 200 K時平板對流換熱進行網格無關性驗證。

圖6(a)和圖6(b)所示為不同網格數時出口溫度分布曲線和當地努塞爾數隨Rex的變化曲線。其中圖6(b)中Nux為式(9)給出的值[14]。

圖6 出口溫度分布及當地Nux變化曲線Fig.6 Temperature and Nux distribution at the outlet

(9)

從圖6可見，3套網格的出口溫度分布曲線和當地努塞爾數隨Rex變化曲線幾乎重合,表明當網格數大于9×104時,計算結果不受網格數影響。因此以下均采用17×104網格進行計算。圖6(b)中當地努塞爾數和式(9)吻合良好,證明采用的數值模擬方案正確。

2.2.2 不同來流溫度下燃氣Pr數的上下邊界對對流換熱的影響

圖7給出了來流溫度分別為1 200、2 200 K時燃氣Pr數分別取上下邊界時平板出口溫度分布曲線。由圖7可見,隨著平板法向距離的增加,不同來流溫度下的出口溫度均單調升高。這是由于流體與平板距離越小,受到對流換熱的影響越大。圖7中燃氣Pr的上下邊界出口溫度分布曲線差別較小,上邊界溫度曲線略大于下邊界,入口溫度為2 200 K時兩者之間差別略大于入口溫度1 200 K的。在固定來流雷諾數不變下,燃氣Pr上下邊界出口溫度曲線差別隨溫度升高逐漸增大,來流2 200 K時兩者最大差別達到2.6%。

圖7 來流溫度不同時燃氣Pr上下邊界出口溫度分布Fig.7 Outlet temperature distributions of the gas at different inlet temperatures under the upper and lower boundaries of the Pr

圖8分別給出了來流溫度1 200 K和2 200 K時燃氣Pr上下邊界的出口速度分布。由于平板壁面黏性滯止,出口速度隨平板法向距離的增加而增加。圖8中燃氣Pr上下邊界出口速度分布曲線差別較大,上邊界對應的速度比下邊界大,而來流溫度為2 200 K時兩者之間差別比來流溫度為1 200 K的大很多。這表明在平板來流雷諾數不變時,燃氣Pr上下邊界對應的出口速度差別隨來流溫度的升高而增大,到2 200 K時兩者的最大差別達17.9%。

圖8 來流溫度不同時燃氣Pr上下邊界出口速度分布Fig.8 Outlet velocity distributions of the gas under the upper and lower boundaries of the Pr of at inlet different temperatures

圖9分別給出了來流溫度分別為1 200、2 200 K時燃氣Pr上下邊界的當地Nux隨平板當地雷諾數Rex變化曲線。由式(9),當固定Rex時,Nux隨著Pr的增大而增大,但由于Pr變化相對Rex較小,所以Nux的變化并不大。圖9中平板的Nux隨著Rex的增加單調增加,Pr上邊界對應的Nux高于下邊界,入口溫度為1 200 K時兩者最大差別為2.7%,2 200 K時則達到4.0%。這表明在雷諾數相同時,燃氣Pr上下邊界之間Nux的差值隨著溫度的升高逐漸增大,其變化趨勢與圖4中Pr上下邊界的變化類似。

圖9 來流溫度不同時燃氣Pr上下邊界對應Nux的變化Fig.9 Local Nux distributions of the gas under the upper and lower boundaries of the Pr at different inlet temperatures over the flat plate

3 渦輪葉柵對流換熱數值模擬

以微小型渦輪發動機渦輪靜葉1/2葉高處平面葉柵對流換熱為例,應用數值模擬的方法,研究不同燃氣輸運特性計算方法對其對流換熱的影響。

3.1 幾何模型及邊界條件

渦輪葉柵幾何模型及邊界條件如圖10所示。其中流場入口及出口高度均為12.5 mm,葉柵弦長為14.5 mm。

圖10 渦輪葉柵幾何模型及邊界條件Fig.10 The geometric model and boundary conditions of turbine cascade

數值模擬的控制方程為N-S方程組及燃氣多組分輸運方程組。入口流場為湍流,湍流模型為可實現性的二方程k-ε方程,湍流度為5%。采用壓力基耦合求解器的SIMPLE算法求解數值模擬的控制方程組。葉柵模型采用速度入口、壓力出口及1 000 K等溫壁面邊界條件。入口溫度分別為1 200、1 500 K；入口速度根據雷諾數的定義及本節給定的數值10 000,根據給定的入口溫度,由1.4節給出的不同溫度下的燃氣輸運特性及燃氣的熱力參數求得；出口壓強為1 atm；流場上下邊界為周期性對稱邊界條件。

3.2 計算結果及分析

3.2.1 網格無關性驗證

采用網格數為14×104、18×104、24×104的3套網格驗證入口溫度1 200 K時的網格無關性。

圖11所示為不同網格數時出口溫度分布曲線和當地Nux隨葉柵Rex(原點為葉柵前緣點)的變化曲線。由圖11(a)可見，3套網格曲線給出的數值模擬結果基本重合,而圖11(b)和圖11(c)中3套網格在葉背和葉盆上的Nux曲線也基本重合。由此,以下計算均采用18×104網格進行。

Fig.11 出口溫度分布及當地Nux變化曲線Fig.11 Temperature and Nux distributions of the skin surface of turbine cascade at the cascade outlet

3.2.2 不同入口溫度下燃氣Pr上下邊界對葉柵對流換熱的影響

圖12分別給出了入口溫度為1 200、1 500 K時燃氣Pr上下邊界對應的出口溫度分布曲線。由圖12可見,出口溫度中間低兩邊高。這是因為流體流過呈周期性分布的葉片時,會受到與葉片對流換熱的影響,離中部葉柵越遠,受到的影響就越小。另外,圖12中燃氣Pr上邊界的出口溫度低于下邊界,入口溫度為1 500 K時兩者之間差別大于入口溫度1 200 K的。這表明燃氣Pr上下邊界出口溫度的差別隨溫度的升高而增大。

圖12 來流溫度不同時燃氣Pr上下邊界出口溫度分布Fig.12 Outlet temperature distributions of the gas under the upper and lower boundaries of the Pr at different inlet temperatures over the turbine cascade

圖13分別給出了入口溫度為1 200、1 500 K時燃氣Pr上下邊界對應的出口速度分布曲線。由圖13可見,出口速度曲線的變化趨勢與圖12的出口溫度分布曲線類似。這是因為燃氣組分的μ會隨著溫度的升高而增加,而來流雷諾數固定,則來流速度會隨μ的增加而變大。另外,圖13中燃氣Pr上邊界的出口速度大于下邊界,入口溫度為1 500 K時上下邊界對應的速度之間差別大于入口溫度1 200 K的。這表明燃氣Pr上下邊界出口速度的差別隨溫度的升高而增大。

Fig.13 來流溫度不同時燃氣Pr上下邊界出口速度分布Fig.13 Outlet velocity distribution of the gas under the upper and lower boundaries of the Pr at different inlet temperatures

圖14分別為入口溫度為1 200、1 500 K時燃氣Pr上下邊界對應的葉柵當地Nux的變化曲線。從圖14中可見,葉背及葉盆的Nux隨葉柵當地Rex的增大而增大,但在葉片尾緣處受附近當地流場變化的影響而變小。圖14中燃氣Pr上邊界對應的Nux大于下邊界。入口溫度為1 200 K時,兩者在葉背及葉盆上最大差別分別為4.74%、4.58%；入口溫度1 500 K時兩者最大差別為9.15%及9.70%。這表明,燃氣Pr不同對對流換熱強度的影響較大,且隨著溫度的升高,燃氣Pr上下邊界對應的Nux的差別逐漸增大。

Fig.14 來流溫度不同時葉背和葉盆上燃氣Pr上下邊界對應的當地Nux變化曲線Fig.14 Local Nux distributions of the gas under the upper and lower boundaries of the Pr over the Turbine cascade at different temperatures

4 結論

提出一種不同燃氣組分輸運特性計算方法對對流換熱影響的快速算法。雖然多組分燃氣的輸運特性計算有多種組合,但由于當地努賽爾數由雷諾數及普朗特數決定,因而僅需得到不同組合下燃氣普朗特數隨溫度變化所有曲線中的上下邊界曲線,并以上下邊界點對應的燃氣輸運特性的值作為輸入條件,由給定的外形進行數值模擬即可評估不同輸運特性計算方法的差別對對流換熱影響程度。而后以平板和渦輪平面葉柵為例研究了不同來流溫度下,不同燃氣輸運特性計算方法組合對對流換熱的影響，得到如下結論。

(1)對平板對流換熱,隨著平板來流溫度的升高(來流雷諾數不變),燃氣普朗特數上下邊界對應的當地努賽爾數的差別變大。當來流溫度從1 200 K升高到2 200 K時,當地努賽爾數差別從2.7%增大到4.0%。

(2)對渦輪平面葉柵對流換熱,燃氣普朗特數上下邊界對應的當地努賽爾數之間的差別隨溫度的升高而增大(來流雷諾數不變)。當來流溫度從1 200 K升高到1 500 K時,當地努賽爾數差別從4.58%增大到9.70%。

(3)在雷諾數相同的條件下,燃氣組分輸運特性計算方法的不同組合對對流換熱強度的影響較大。多組分燃氣在燃燒室出口溫度不同時各組分的質量分數也不同,并且單組分氣體輸運特性有多種計算方法,導致在研究熱機及推進系統的對流換熱問題中的燃氣輸運特性和Pr時,有眾多的計算方法的組合。因此在進行對流換熱計算時,可快速給出不同燃氣輸運特性計算方法所有組合的努賽爾數的變化區間的左右邊界,為工程應用提供參考。