葉片外換熱計(jì)算程序STAN5收斂性改進(jìn)

婁德倉(cāng)

(中國(guó)燃?xì)鉁u輪研究院，四川 成都 610500)

1 引言

渦輪葉片燃?xì)膺叺膿Q熱系數(shù)，可通過(guò)數(shù)值求解流過(guò)型面的粘性流控制微分方程得到。這種解法的優(yōu)點(diǎn)是：可得到壁面附近氣流的速度和溫度分布，且在此基礎(chǔ)上求得的壁面摩阻系數(shù)和換熱系數(shù)較為準(zhǔn)確；能較方便地應(yīng)用近年來(lái)發(fā)展的湍流模型；計(jì)算中用到的經(jīng)驗(yàn)公式少，適應(yīng)性較廣。由于求解邊界層方程的計(jì)算量較求解N-S方程的小得多，因而在工程計(jì)算中得到了較廣泛的應(yīng)用。Crawford等人于1975年采用Patankar-Spalding方法編制了邊界層微分方程計(jì)算程序STAN5[1]，隨后又對(duì)該程序進(jìn)行了諸多改進(jìn)。目前使用的TAXSTAN程序是其最新版本[2]，能考慮葉片表面氣膜出流對(duì)外換熱的影響。

上世紀(jì)80年代，STAN5程序被西工大劉松齡教授等引入國(guó)內(nèi)，經(jīng)改進(jìn)后用于渦輪葉片外換熱計(jì)算，結(jié)合氣膜修正等經(jīng)驗(yàn)公式，能用于帶氣膜冷卻渦輪葉片的熱分析。鄧化愚等人[3]還對(duì)程序進(jìn)行了改進(jìn)，主要是加入了湍流度間隙因子，使程序能考慮來(lái)流湍流度對(duì)外換熱的影響。蔡毅等人[4]通過(guò)微分法、積分法等的對(duì)比，發(fā)現(xiàn)STAN5程序在計(jì)算方法及精度上更具優(yōu)勢(shì)。但在應(yīng)用STAN5程序進(jìn)行外換熱計(jì)算時(shí)，會(huì)出現(xiàn)不收斂等問(wèn)題，這主要是由外邊界的逆壓梯度分布、計(jì)算網(wǎng)格及邊界層卷吸量的控制等因素引起。因此，計(jì)算中經(jīng)常要對(duì)渦輪氣動(dòng)數(shù)據(jù)進(jìn)行人工干預(yù)，這會(huì)造成計(jì)算結(jié)果的人為誤差。為此，本文對(duì)STAN5程序計(jì)算中存在的收斂性等問(wèn)題進(jìn)行分析，并采取相應(yīng)措施進(jìn)行改進(jìn)。

2 邊界層微分方程

描述軸對(duì)稱可壓縮流的邊界層流動(dòng)方程，包括連續(xù)性方程、動(dòng)量方程和能量方程。其中連續(xù)性方程可表示為[5]：

式中：坐標(biāo)x、y的方向分別與物體的型線相切和垂直，u和v分別為x、y方向的分速度，ρ為氣流密度，r為徑向坐標(biāo)。

動(dòng)量方程為：

式中：p為壓力，μ為動(dòng)力粘性系數(shù)，ρ----u′v′為湍流剪切應(yīng)力，u′和v′分別為x、y方向的湍流脈動(dòng)速度，z為單位質(zhì)量流體上的徹體力。

總焓形式的能量方程為：

式中：I為氣體總焓，λ為氣體導(dǎo)熱系數(shù)，cp為氣體定壓比熱，i為氣體靜焓；I′為總焓脈動(dòng)值，S為單位質(zhì)量流體中的源項(xiàng)。

根據(jù)以上方程，用已知的邊界條件和初始條件，便可解出ρ、u、v、I在邊界層內(nèi)的分布。對(duì)于湍流邊界層，還須補(bǔ)充求湍流剪切力的關(guān)系式。

3 STAN5程序求解

3.1 初始條件

求解邊界層微分方程過(guò)程中，常將駐點(diǎn)區(qū)域?qū)恿鬟吔鐚幼鳛橛?jì)算初始條件。在STAN5程序中，根據(jù)Pohlhausen公式獲得層流邊界層內(nèi)的速度分布：

式中：δ為邊界層厚度，Λ為Pohlhausen參數(shù)，ue為邊界層外氣流速度。

邊界層起始點(diǎn)溫度分布(T)可通過(guò)溫度分布與速度分布間的線性關(guān)系得到，關(guān)系式為：

式中：Tw為壁面溫度，Te為邊界層外氣流溫度。

3.2 湍流模型

現(xiàn)有的STAN5程序中采用的湍流模型有混合長(zhǎng)度模型、單方程湍流模型，同時(shí)考慮了湍流度修正。其中較為常用的是Prandtl混合長(zhǎng)度模型。

3.3 邊界層網(wǎng)格變換

邊界層微分方程求解方法很多，STAN5程序采用的是Patankar和Spalding提出的流函數(shù)坐標(biāo)變換法[6]，首先從x、y坐標(biāo)變換為x、ω(x,y)坐標(biāo)，對(duì)邊界層微分方程進(jìn)行求解，然后再將結(jié)果返回到x、y坐標(biāo)下。ω(x,y)及其與y的關(guān)系可用下式表示：

STAN5程序中的網(wǎng)格，其徑向分布采用基于式(7)建立的速度分布與網(wǎng)格密度關(guān)系。為保證邊界層底層網(wǎng)格有足夠的密度，采用了滑移點(diǎn)控制方法。

式中：c1為常數(shù)，β為滑移點(diǎn)參數(shù)。

4 STAN5在計(jì)算中存在的主要問(wèn)題及改進(jìn)

以某渦輪葉片外換熱計(jì)算過(guò)程為例。STAN5程序在計(jì)算葉片外部換熱時(shí)，首先將葉片表面的等熵馬赫數(shù)分布作為邊界層計(jì)算中的外部速度邊界。通過(guò)求解邊界層微分方程，獲得葉片表面換熱系數(shù)分布。典型的渦輪葉片表面等熵馬赫數(shù)分布如圖1所示。從圖中看，在渦輪葉片葉盆側(cè)前緣擴(kuò)張段和尾緣激波區(qū)存在明顯的逆壓梯度，而這些位置通常會(huì)造成外換熱計(jì)算不收斂。因此在外換熱計(jì)算過(guò)程中，通常要對(duì)外部勢(shì)流速度進(jìn)行人為調(diào)整，避免出現(xiàn)逆壓梯度(如圖1中虛框內(nèi)所示)。由于邊界層微分方程為橢圓型微分方程(即下游計(jì)算結(jié)果完全受上游的影響)，因此該調(diào)整會(huì)導(dǎo)致計(jì)算結(jié)果人為誤差較大，影響外換熱計(jì)算結(jié)果的精確性。

圖1 渦輪葉片表面等熵馬赫數(shù)分布Fig.1 Isentropic Mach number distribution of blade surface

4.1 收斂性問(wèn)題及分析

根據(jù)邊界層內(nèi)部流動(dòng)方程的量級(jí)分析，邊界層內(nèi)部的壓力分布可簡(jiǎn)化為式(8)。其中沿流向的壓力梯度決定于外部主流的速度分布。因此，在外部勢(shì)流存在較強(qiáng)逆壓梯度且邊界層底部速度較小時(shí)，計(jì)算中就會(huì)出現(xiàn)速度負(fù)值。而邊界層網(wǎng)格又與速度分布相關(guān)，最終導(dǎo)致計(jì)算發(fā)散。

通常，邊界層內(nèi)部出現(xiàn)速度負(fù)值時(shí)，STAN5程序自身會(huì)對(duì)速度負(fù)值進(jìn)行調(diào)整(如式(9))，同時(shí)減少邊界層外部的卷吸流量，重新求解方程。

如圖2所示，修正后的速度分布與邊界層內(nèi)真實(shí)的速度分布有很大差別。若繼續(xù)計(jì)算，仍會(huì)出現(xiàn)速度負(fù)值，導(dǎo)致計(jì)算發(fā)散。有時(shí)雖能得到收斂結(jié)果，但所采用的方法已改變了邊界層外部的氣流速度分布規(guī)律(圖3)，并且這種收斂方式會(huì)對(duì)換熱系數(shù)(α)產(chǎn)生較大影響(圖4)。

4.2 收斂性的兩點(diǎn)改進(jìn)

通常在邊界層方程計(jì)算中，若出現(xiàn)速度負(fù)值，表明邊界層已破壞，這時(shí)邊界層方程不再適用。程序自身的修正實(shí)質(zhì)上是假設(shè)該區(qū)域內(nèi)仍然為邊界層型流動(dòng)。在渦輪流道中，大部分區(qū)域?qū)儆陧槈禾荻鹊倪吔鐚恿鲃?dòng)，因此該局部區(qū)域內(nèi)采用此假設(shè)仍然合理。本文在此基礎(chǔ)上對(duì)逆壓梯度下的速度分布u(y)及其網(wǎng)格密度進(jìn)行修正。

4.2.1 逆壓梯度下的速度分布

針對(duì)葉片外部速度場(chǎng)急速變化時(shí)產(chǎn)生的較強(qiáng)逆壓梯度，可采取兩種方法進(jìn)行速度修正。

一是在當(dāng)前計(jì)算站速度梯度出現(xiàn)負(fù)值時(shí)，對(duì)速度分布進(jìn)行修正，如圖5(a)所示。當(dāng)前計(jì)算點(diǎn)xi的速度分布形式ui(y)，與前一計(jì)算點(diǎn)xi-1的速度分布形式ui-1(y)保持相似，即：

圖5 邊界層內(nèi)部速度分布修正Fig.5 Modification of the velocity profile inside the boundary layer

結(jié)果表明該方法不能完全避免速度負(fù)值問(wèn)題。

二是在當(dāng)前計(jì)算站存在逆壓梯度且速度二次梯度出現(xiàn)負(fù)值時(shí)，對(duì)速度分布進(jìn)行修正，如圖5(b)所示。當(dāng)前計(jì)算點(diǎn)xi的速度分布形式ui(y)，與前一計(jì)算點(diǎn)xi-1的速度分布形式ui-1(y)保持相似，并保持外邊界的速度值uei不變，即：

結(jié)果表明該方法能明顯改善計(jì)算不收斂。

4.2.2 基于壓力梯度變化的自適應(yīng)網(wǎng)格

葉盆側(cè)尾緣激波區(qū)馬赫數(shù)變化劇烈，壓力梯度迅速由順壓梯度變?yōu)槟鎵禾荻?，因此該區(qū)域?qū)α飨蚓W(wǎng)格變化較敏感。雖采取了上述速度分布修正，但在有些馬赫數(shù)分布下仍會(huì)出現(xiàn)計(jì)算不收斂。通常邊界層在尾緣區(qū)域較厚，沿流向的空間步長(zhǎng)較大(Δx=0.3δ)，這時(shí)會(huì)出現(xiàn)兩種情況：一是當(dāng)前節(jié)點(diǎn)位于最大(或較大)逆壓梯度點(diǎn)，導(dǎo)致計(jì)算發(fā)散；二是當(dāng)前節(jié)點(diǎn)跨過(guò)最大(或較大)逆壓梯度點(diǎn)，計(jì)算穩(wěn)定收斂。

為使邊界層內(nèi)部沿流向網(wǎng)格步長(zhǎng)與主流壓力梯度相聯(lián)系，在逆壓梯度下，采用式(12)計(jì)算網(wǎng)格步長(zhǎng)：

式中：下標(biāo)old表示上一步長(zhǎng)的參數(shù)，Δx0為當(dāng)前計(jì)算站下的初始步長(zhǎng)，Δx1為經(jīng)壓力梯度修正后的計(jì)算步長(zhǎng)。圖6對(duì)比了壓力梯度修正前后葉片尾緣附近的網(wǎng)格，可見(jiàn)修正后的網(wǎng)格明顯變密。

經(jīng)上述修正，STAN5程序計(jì)算中的不收斂現(xiàn)象得到解決。修正過(guò)程中，保持外邊界速度不變，就能保證計(jì)算結(jié)果誤差在合理范圍內(nèi)。

圖6 修正前后的沿流向網(wǎng)格Fig.6 Comparison of the stream-wise mesh before and after modification

5 計(jì)算結(jié)果分析

圖7 外換熱計(jì)算結(jié)果分布Fig.7 External heat transfer coefficient of the blade

根據(jù)圖1中的馬赫數(shù)分布，改進(jìn)后的STAN5程序外換熱計(jì)算結(jié)果與人工修正馬赫數(shù)分布計(jì)算結(jié)果間的對(duì)比如圖7所示。從圖中看出，兩者間的外換熱分布趨勢(shì)一致，但在葉盆側(cè)，兩者的計(jì)算值相差約±4%，最大達(dá)20%。采用新型網(wǎng)格的計(jì)算結(jié)果中，對(duì)尾緣激波區(qū)域(如圖中葉背尾緣)的換熱系數(shù)分布更加詳細(xì)。圖8為邊界層內(nèi)部的速度分布云圖。

圖8 邊界層內(nèi)部的速度分布云圖Fig.8 Contour of velocity inside the boundary layer

6 結(jié)論

本文分析了渦輪葉片外換熱計(jì)算程序STAN5中的求解方法，包括邊界層求解的微分方程、計(jì)算初始條件及網(wǎng)格變換等。針對(duì)程序在外流場(chǎng)存在較強(qiáng)逆壓梯度時(shí)所出現(xiàn)的計(jì)算不收斂及錯(cuò)誤收斂等現(xiàn)象，對(duì)該程序進(jìn)行了兩點(diǎn)改進(jìn)：修正了逆壓梯度下邊界層內(nèi)部速度分布，建立了邊界層內(nèi)部沿流向的網(wǎng)格分布與主流壓力梯度間的聯(lián)系。改進(jìn)后的程序有效地避免了計(jì)算發(fā)散等問(wèn)題。算例分析結(jié)果表明，本文提出的改進(jìn)方法與以前的人工調(diào)整方式獲得的結(jié)果相比，外換熱系數(shù)分布趨勢(shì)一致，局部誤差約為4%，最大誤差達(dá)20%。

[1]Crawford M E，Kays W M.STAN5-A Program for Numeri?cal Computation of Two-Dimensional Internal and Exter?nal Boundary Layer Flows[R].NASA CR-2742，1976.

[2]Kays W M，Crawford M E，Weigand B.Convective Heat and Mass Transfer[M].4thed.McGraw-Hill，2005.

[3]鄧化愚，蔡 毅，劉玉芳.渦輪氣冷葉片局部外換熱系數(shù)計(jì)算方法的改進(jìn)[J].燃?xì)鉁u輪試驗(yàn)與研究，1989，2(2)：1—5.

[4]蔡 毅，劉松齡.渦輪葉片外換熱系數(shù)計(jì)算方法和比較[J].燃?xì)鉁u輪試驗(yàn)與研究，1995，8(1)：1—8.

[5]航空發(fā)動(dòng)機(jī)手冊(cè)總編委會(huì).航空發(fā)動(dòng)機(jī)手冊(cè)：第16冊(cè)——空氣系統(tǒng)及傳熱分析[K].北京：航空工業(yè)出版社，2001.

[6]Patankar S V，Spalding D B.Heat and Mass Transfer in Boundary Layers:A General calculation Procedure[M].London：International Textbook Company Ltd.，1970.