陳 杰,張家騏,歐吉輝

(天津大學(xué) 機(jī)械工程學(xué)院 高速空氣動(dòng)力學(xué)研究室,天津 300072)

0 引 言

臨近空間高超聲速飛行器的研發(fā)在中國及世界范圍內(nèi)已受到越來越多的關(guān)注,但同時(shí)也面臨著巨大的技術(shù)挑戰(zhàn)和新的科學(xué)問題[1-3]。周恒和張涵信在他們所寫的《空氣動(dòng)力學(xué)的新問題》一文[1]中分析了在情況下現(xiàn)有的空氣動(dòng)力學(xué)的不足,提出了若干需要考慮的新問題,其中之一就是如何考慮流場中可能出現(xiàn)的局部稀薄氣體效應(yīng)。他們分析了臨近空間高超聲速飛行器的邊界層,指出:由于邊界層中溫度很高,再加上臨近空間飛行器的飛行高度很高,因此邊界層中的氣體密度可能很低,即氣體分子自由程較大。另一方面,邊界層中的剪切很強(qiáng),在這兩個(gè)因素的影響下,邊界層中法向相隔一個(gè)分子自由程的兩層流體的速度差和當(dāng)?shù)氐姆肿舆\(yùn)動(dòng)的平均速度相比不是一個(gè)小量,這會(huì)導(dǎo)致分子運(yùn)動(dòng)的速度分布函數(shù)顯著偏離Maxwell分布,使得現(xiàn)有的各種基于實(shí)際分子自由運(yùn)動(dòng)速度的分布函數(shù)對(duì)Maxwell分布僅有小的偏離假設(shè)下提出的分析計(jì)算方法不再可靠,需要尋找新的方法。根據(jù)這一思想,陳杰和趙磊在文獻(xiàn)[2]中提出了對(duì)應(yīng)的參數(shù)Zh,用以刻畫氣體稀薄效應(yīng)對(duì)剪應(yīng)力的影響。

參數(shù)Zh定義為:相隔為一個(gè)分子自由程的兩層流體的宏觀速度差和分子熱運(yùn)動(dòng)速度之比。這是一無量綱參數(shù),即:

在臨近空間高超聲速飛行器的設(shè)計(jì)中,不僅要考慮氣動(dòng)力,還要考慮氣動(dòng)熱[3,6-10]。對(duì)于長時(shí)間在近空間中以高超聲速機(jī)動(dòng)飛行的飛行器,精確的熱環(huán)境預(yù)測對(duì)低冗度熱防護(hù)設(shè)計(jì)十分重要[6]。臨近空間高超聲速飛行器氣動(dòng)加熱受稀薄氣體效應(yīng)、高溫真實(shí)氣體效應(yīng)、飛行器復(fù)雜外形帶來的復(fù)雜流場結(jié)構(gòu)等多種因素的影響[7-10]。工程上往往針對(duì)駐點(diǎn)熱流基于半經(jīng)驗(yàn)參數(shù)構(gòu)建橋函數(shù)等方式考慮稀薄非平衡計(jì)算[11]。在數(shù)值計(jì)算方面也發(fā)展了多種考慮稀薄效應(yīng)的宏觀方程模型[11-15]。

本文將在之前研究強(qiáng)剪切流中稀薄效應(yīng)對(duì)應(yīng)力計(jì)算工作的基礎(chǔ)上[2,4-5],進(jìn)一步考慮氣體稀薄效應(yīng)對(duì)熱流計(jì)算的影響,將研究以下問題:1)存在強(qiáng)溫度梯度時(shí),氣體的稀薄效應(yīng)如何分析和處理,分子速度分布函數(shù)是否還是Maxwell分布,如果不是,刻畫偏離程度的無量綱參數(shù)是什么;2)如果該參數(shù)的值大了,會(huì)如何影響熱流計(jì)算?傅里葉導(dǎo)熱定律是否成立?是否能定量地找到等效導(dǎo)熱系數(shù)和此參數(shù)的依賴關(guān)系。研究的路線和文獻(xiàn)[2]相同,即以DSMC為手段尋找氣體有稀薄效應(yīng)下的熱傳導(dǎo)規(guī)律。3)初步驗(yàn)證將該規(guī)律納入CFD計(jì)算中是否可以有效改善受稀薄效應(yīng)的熱流計(jì)算?

1 計(jì)算模型

1.1 一維導(dǎo)熱問題

為單純研究存在強(qiáng)溫度梯度時(shí)氣體稀薄效應(yīng)對(duì)熱流計(jì)算的影響,選取無宏觀流動(dòng)的導(dǎo)熱問題進(jìn)行分析。我們這里想要獲得的是對(duì)于具有一定稀薄程度的氣體在存在強(qiáng)溫度梯度時(shí),傅里葉定律所計(jì)算的熱流誤差,而不考慮壁面引起的效應(yīng)。如圖1(a)所示,該研究對(duì)象為距壁面較遠(yuǎn)且存在強(qiáng)溫度梯度的流體,即在厚度L的薄層內(nèi)存在較大的溫度差Th－Tc。經(jīng)典的兩平板間導(dǎo)熱問題如圖1(b)所示,相距L的兩板分別保持不同的溫度Tc、Th,壁面邊界條件通常采用Maxwell完全漫反射條件。對(duì)于圖1(b)所示問題,可基于兩板間距定義傳統(tǒng)Knudsen數(shù)KnL＝λ/L。如果KnL較小,如KnL＜0.1,兩板間中心區(qū)域可認(rèn)為不受固壁邊界條件影響,等價(jià)于圖1(a)中所要研究的薄層。但隨著KnL增大,壁面效應(yīng)越明顯,中心區(qū)域的熱流計(jì)算將既受到強(qiáng)溫度梯度的影響,也受到壁面影響。雖然可以通過增大兩板間的距離減弱壁面效應(yīng),但是當(dāng)兩板距離L增大后,如在中心區(qū)域?qū)崿F(xiàn)相同的溫度梯度d T/d x,兩板的壁面溫度之差就要相應(yīng)地增大,很難實(shí)現(xiàn)。如保持中心溫度不變,一端溫度甚至可能小于0°,顯然不合理。如要求一端保持不低于0°,則另一端溫度就非常高。即便這樣也不容易實(shí)現(xiàn)較大ZhT時(shí)的情況。因此本文將通過在虛線處施加合適的邊界實(shí)現(xiàn)圖1(a)中所示問題的模擬。

圖1 一維導(dǎo)熱問題Fig.1 1D conduction

1.2 DSMC方法

本文采用基于Bird提出的標(biāo)準(zhǔn)DSMC算法[16]建立的自編程程序。作為第一步,為避免高溫下可能引起的分子振動(dòng)能被激發(fā),甚至導(dǎo)致分子離解等復(fù)雜問題掩蓋氣體稀薄效應(yīng)的本質(zhì)問題,本文考慮單原子氣體氬氣作為模擬氣體,采用硬球模型。氬氣分子質(zhì)量為6.63×10-26kg,直徑為3.659×10-10m。

DSMC方法使用大量的模擬分子模擬真實(shí)的氣體,用一個(gè)模擬分子代表一定數(shù)量的真實(shí)氣體分子,該方法的本質(zhì)是在時(shí)間步長Δt內(nèi),對(duì)分子的運(yùn)動(dòng)和分子間的碰撞進(jìn)行解耦,要求時(shí)間步長小于平均碰撞時(shí)間,網(wǎng)格尺度小于平均自由程。因此,為保證計(jì)算有效導(dǎo)熱系數(shù)的足夠精確,采用了文獻(xiàn)[17]中建議的網(wǎng)格大小和時(shí)間步長。為獲得統(tǒng)計(jì)噪聲較小的熱流以及更準(zhǔn)確的等效導(dǎo)熱系數(shù),模擬使用了大量的模擬分子,初始每個(gè)網(wǎng)格內(nèi)分子數(shù)為1×100個(gè),對(duì)導(dǎo)熱穩(wěn)定后的106個(gè)時(shí)間步進(jìn)行統(tǒng)計(jì)平均。

1.3 稀薄效應(yīng)刻畫參數(shù)

周恒和張涵信[1]指出:相隔一個(gè)分子自由程的兩層間的分子碰撞是產(chǎn)生剪切力和熱傳導(dǎo)的主要機(jī)制。在剪切流中,刻畫稀薄效應(yīng)的無量綱參數(shù)Zh為一個(gè)分子自由程上宏觀速度的差與微觀分子最可幾速度之比。該參數(shù)有效表征了剪切流中的非平衡。根據(jù)同樣的思想,這里將該參數(shù)類推至溫度,流體的溫度是分子平均熱速度的統(tǒng)計(jì)。如果一個(gè)分子自由程上分子熱運(yùn)動(dòng)速度的改變量與熱運(yùn)動(dòng)平均速度本身相比不再是小量,那分子速度分布也不再會(huì)是Maxwell分布。

因此,定義ZhT為相隔一個(gè)分子自由程的兩層流體的特征分子熱速度梯度和當(dāng)?shù)氐钠骄鶡崴俣戎葹閆hT:

需要指出的是ZhT的形式雖然與文獻(xiàn)[18]中的KnGLL-T只差一個(gè)常數(shù),但是意義不同。這里是兩個(gè)速度或者溫度之比,而文獻(xiàn)里KnGLL-T的意義為兩個(gè)長度之比,即分子平均自由程與溫度的標(biāo)尺長度之比。ZhT參數(shù)的物理意義將在2.1節(jié)中通過對(duì)速度分布函數(shù)的分析進(jìn)一步詳細(xì)解釋。

1.4 邊界條件

這一節(jié)探討邊界條件的實(shí)施。如1.1節(jié)所示,我們需要采用合適的邊界條件模擬距壁面較遠(yuǎn)且存在強(qiáng)溫度梯度的一層流體。在稀薄強(qiáng)溫度梯度條件下,氣體顯然不可能處于平衡態(tài),即在DSMC模擬的邊界處不能通過平衡態(tài)Maxwell分布向計(jì)算區(qū)域加入粒子。為此,我們提出了新的邊界條件以實(shí)現(xiàn)邊界處所設(shè)定的溫度以及粒子所處的非平衡特性。

該邊界條件的主要思想是對(duì)于不受壁面影響的導(dǎo)熱問題,模擬區(qū)域邊界處與中心點(diǎn)具有相似的非平衡性,即分布函數(shù)形狀以及偏離Maxwell分布的程度相同,因此通過中心點(diǎn)的粒子信息在邊界處進(jìn)行重構(gòu),使其與中心點(diǎn)具有相似非平衡態(tài)分布函數(shù),但兩組粒子的統(tǒng)計(jì)溫度不同。在下文將此邊界稱為“非平衡克隆邊界”。其具體實(shí)現(xiàn)方法分以下幾個(gè)步驟:在每一次時(shí)間步推進(jìn)后,1)統(tǒng)計(jì)邊界網(wǎng)格中存在的粒子數(shù),用N表示;2)從中間網(wǎng)格隨機(jī)抽取N個(gè)分子復(fù)制到邊界網(wǎng)格。對(duì)低溫端,由于對(duì)應(yīng)的溫度低,所以一般N要大于中間網(wǎng)格的粒子數(shù),需要另外隨機(jī)抽取若干粒子,補(bǔ)足差數(shù);3)統(tǒng)計(jì)抽取的這N個(gè)粒子重新組成的粒子團(tuán)的平均速度。因?yàn)殡S機(jī)性,重新組成的粒子團(tuán)平均速度可能不再是0,有很微小的偏差。因此需要將這個(gè)小的偏差速度減掉,以免造成邊界處不正確的動(dòng)量輸入;4)統(tǒng)計(jì)N個(gè)粒子的溫度T1,根據(jù)壁面設(shè)定的溫度Tw和這個(gè)統(tǒng)計(jì)溫度比重整粒子的速度,即粒子速度乘以。這樣,壁面網(wǎng)格內(nèi)的分子平均速度為0,溫度為Tw,分子的速度分布函數(shù)與中間網(wǎng)格類似。

圖2 概率密度分布函數(shù)Fig.2 Probability density function

為驗(yàn)證此邊界,我們分別模擬了圖1所示兩種情況。對(duì)于圖1(a)所示模型,Tc＝700 K,Th＝1200 K,L為中心點(diǎn)平均自由程的5倍,兩邊均用克隆邊界。對(duì)于圖1(b)所示模型,Tc＝400 K,Th＝2000 K,中心點(diǎn)附近溫度為1200 K,為避免壁面影響,L為中心點(diǎn)平均自由程的10倍。圖2給出了兩個(gè)計(jì)算模型相同溫度1200 K處的概率密度分布函數(shù)(Probability Density Function,PDF),并與1200 K的Maxwell分布對(duì)比?？梢钥吹綀D1(a)模型中克隆邊界條件所實(shí)現(xiàn)的分布函數(shù)與圖1(b)模型中心點(diǎn)的溫度函數(shù)重合,且相同程度地偏離了平衡態(tài),驗(yàn)證了克隆邊界可以重構(gòu)不受壁面影響的非平衡態(tài)分布函數(shù)。峰值微弱的差別來自于溫度的差,克隆邊界處溫度為1200 K,中心點(diǎn)附近所取的點(diǎn)溫度為1197 K。

1.5 等效導(dǎo)熱系數(shù)

本文最終要獲得的是ZhT對(duì)熱流計(jì)算的影響,求得等效的導(dǎo)熱系數(shù),其定義為:

即DSMC獲得的熱流比上當(dāng)?shù)氐臏囟忍荻?可理解為氣體在線性本構(gòu)關(guān)系假設(shè)下所表現(xiàn)的等效導(dǎo)熱系數(shù),它是傳統(tǒng)的導(dǎo)熱系數(shù)乘以修正系數(shù)Ak(ZhT)。

2 結(jié)果分析

2.1 速度分布函數(shù)

我們計(jì)算了不同ZhT的參數(shù)下計(jì)算區(qū)域中心點(diǎn)的概率密度分布函數(shù)(PDF)。隨著ZhT的增加,分布函數(shù)逐漸偏離Maxwell分布。不同于剪切流中分布函數(shù)呈現(xiàn)單峰、雙峰、三峰等多種形態(tài)[2],導(dǎo)熱問題中不同ZhT的參數(shù)下分布函數(shù)偏離的形態(tài)是一致的,如圖3所示,只是偏離程度不同。圖3給出了ZhT＝0.03和ZhT＝0.1時(shí)的概率密度分布函數(shù)。ZhT＝0.03時(shí)分布函數(shù)的峰值已經(jīng)開始向右傾斜,ZhT＝0.1則明顯偏離。

從分布函數(shù)與Maxwell分布對(duì)比來看,粒子運(yùn)動(dòng)的最可幾速度不再是零而變?yōu)檎Ｆ湓蚩煞治鋈缦?設(shè)圖4中的A、B線分別為一個(gè)網(wǎng)格的兩個(gè)邊界,溫度梯度方向?yàn)閺腁向B,即B處的溫度高于A處的溫度,用NA、NB分別表示在網(wǎng)格線上進(jìn)出的粒子數(shù)目。由于氣體宏觀靜止,因此在網(wǎng)格線處,單位時(shí)間從左往右和從右往左穿過的粒子數(shù)應(yīng)相等,即圖中所示的＝,＝。設(shè)c－為分子熱運(yùn)動(dòng)平均速度,則可認(rèn)為在一個(gè)時(shí)間步長d t內(nèi),位于B或A兩側(cè)一個(gè)寬為的范圍內(nèi)的正向或反向運(yùn)動(dòng)的粒子在d t時(shí)間內(nèi)均可分別向左或向右穿過網(wǎng)格線,對(duì)初始的Maxwell分布來說,正向和反向的粒子數(shù)應(yīng)相等。而位于該范圍的粒子數(shù)則和密度與成正比。由于密度和溫度成反比,而則與成正比,因此在d t時(shí)間內(nèi)穿過網(wǎng)格的粒子數(shù)應(yīng)該和1/成正比。由于B處的溫度比A處的溫度高,所以B處穿過網(wǎng)格的粒子數(shù)小于A處穿過網(wǎng)格的粒子數(shù),也就是圖中所標(biāo)示的,即有凈流入的正向粒子和凈流出的反向粒子,使得分布函數(shù)偏于正的一側(cè)。

圖3 概率密度分布函數(shù)Fig.3 Probability density function with different Zh T parameters

圖4 粒子輸運(yùn)示意圖Fig.4 Schematic diagram of particle transport

根據(jù)以上分析,單位時(shí)間內(nèi)一個(gè)網(wǎng)格的正向粒子凈流入量和通過中心界面的正向粒子流量之比,即:

可以刻畫分布函數(shù)偏離Maxwell分布的程度。而這一參數(shù)和公式(3)的定義一致。

我們還可從DSMC的計(jì)算過程中做統(tǒng)計(jì)以證實(shí)這一趨勢。在計(jì)算中心點(diǎn)網(wǎng)格溫度為1000 K,網(wǎng)格內(nèi)模擬粒子數(shù)在100個(gè)左右,等導(dǎo)熱穩(wěn)定后我們對(duì)中心點(diǎn)網(wǎng)格低溫側(cè)和高溫側(cè)進(jìn)出的粒子在1×106時(shí)間步上進(jìn)行統(tǒng)計(jì)平均。如圖3(a)所示的ZhT＝0.03的情況,中心點(diǎn)網(wǎng)格低溫側(cè)進(jìn)出的粒子均為5.4088個(gè),而高溫側(cè)進(jìn)出的均為5.3893,即網(wǎng)格內(nèi)總的進(jìn)出粒子相等,但低溫側(cè)進(jìn)入的比右側(cè)多0.36%。圖3(b)所示的ZhT＝0.1的情況,低溫側(cè)進(jìn)入的比高溫側(cè)多0.42%,當(dāng)ZhT＝0.36時(shí)低溫側(cè)進(jìn)入的比高溫側(cè)多0.52%。因此,ZhT越大,網(wǎng)格中來自低溫側(cè)的粒子比來自高溫側(cè)的粒子越多。不過DSMC的結(jié)果是最終平衡態(tài)的結(jié)果,當(dāng)數(shù)量更多的正向粒子進(jìn)入網(wǎng)格后,會(huì)和反向運(yùn)動(dòng)的粒子發(fā)生碰撞,使得正反向的粒子數(shù)趨于相等(正向粒子數(shù)雖然較多,但平均速度較小,網(wǎng)格內(nèi)粒子之間相互碰撞的總的效果是使得粒子數(shù)和大小趨于相等,和輸入粒子的作用相反)。這是一個(gè)十分復(fù)雜的過程,很難用簡單的物理過程說清楚。

2.2 壁面效應(yīng)分析

本節(jié)主要分析壁面效應(yīng)的影響,以及驗(yàn)證克隆邊界是否能夠有效消除壁面影響。對(duì)同一問題,如果使用完全漫反射條件,即所模擬的問題為受限于間隔距離很小的兩平板間的導(dǎo)熱問題,壁面附近氣體的溫度與給定的壁面溫度不相等,即所謂的溫度跳躍,壁面附近努森層中的溫度分布也呈現(xiàn)非線性變化,如圖5黑線所示。如果采用克隆邊界,溫度分布則如黑點(diǎn)所示,可基本實(shí)現(xiàn)所施加的溫度梯度,且成近似線性分布。在靠近壁面的一個(gè)點(diǎn)處有一個(gè)小的間斷,但這并不影響我們?cè)谥行膮^(qū)域的分析。產(chǎn)生該小的間斷的原因是,在計(jì)算中,網(wǎng)格尺度是當(dāng)?shù)刈杂沙痰?/5以下,按理一個(gè)網(wǎng)格兩側(cè)5個(gè)網(wǎng)格內(nèi)的粒子都有可能不經(jīng)過粒子間碰撞直接到達(dá)該網(wǎng)格對(duì)熱量交換做出貢獻(xiàn)。但邊界網(wǎng)格內(nèi)的那個(gè)網(wǎng)格在邊界方向只受到臨近一個(gè)網(wǎng)格的影響,而損失了一部分粒子對(duì)它的影響。比如低溫端右側(cè)的第二個(gè)網(wǎng)格,它左側(cè)只受到第一個(gè)網(wǎng)格內(nèi)進(jìn)來的粒子的影響,而實(shí)際情況應(yīng)該是有從更遠(yuǎn)處的網(wǎng)格中更低溫度的粒子參與能量交換。失去了這些更低溫度粒子的參與,自然要求相鄰網(wǎng)格的更大溫差來彌補(bǔ)。高溫端的情況類似。

圖5 兩種邊界條件的溫度分布對(duì)比(Kn L＝0.2)Fig.5 Comparison of temperature distribution for two different boundary conditions(Kn L＝0.2)

我們還進(jìn)一步分析了壁面效應(yīng)對(duì)于等效導(dǎo)熱系數(shù)計(jì)算的影響。為此,在固定ZhT的同時(shí),改變兩壁面間的間距,看大到什么程度后其等效導(dǎo)熱系數(shù)不再改變(大小以間距L相當(dāng)于多少個(gè)中心位置的自由分子程數(shù)λ0計(jì),即KnL＝λ0/L)。對(duì)初始ZhT分別為0.04和0.08的情況采用漫反射邊界和克隆邊界做了兩組算例。圖6給出了所得的導(dǎo)熱系數(shù)k的修正值A(chǔ)k隨ZhT的變化。圖中空心圓和空心三角為完全漫反射邊界計(jì)算的結(jié)果,并標(biāo)出了相應(yīng)的KnL值。紅色實(shí)心圓和實(shí)心三角為克隆邊界計(jì)算的結(jié)果。對(duì)比發(fā)現(xiàn),漫反射邊界計(jì)算所得的Ak隨兩板間距減小而急劇減小,而隆邊界條件計(jì)算的Ak的值則幾乎不受兩板間距的影響。對(duì)我們所關(guān)心的情況,用克隆邊界條件時(shí),計(jì)算域最多只要2.5倍的中心位置處的分子自由程就夠了,而對(duì)兩端為固壁的情況,則即使ZhT僅為0.04,計(jì)算域也要達(dá)到10個(gè)分子自由程以上才可以。而由于溫度有不能小于零的限制,對(duì)兩端固壁情況,有時(shí)不可能滿足要求。

圖6 兩種邊界條件計(jì)算的導(dǎo)熱系數(shù)修正值對(duì)比Fig.6 Correctional coefficients of thermal conductivity calculated by two boundary conditions with parameter Zh T

這部分結(jié)果說明:采用壁面漫反射條件在計(jì)算域不夠大時(shí)不能獲得無量綱參數(shù)ZhT和修正系數(shù)Ak之間的規(guī)律性變化;而克隆邊界可有效消除固壁影響,只需要較小的計(jì)算域就能消除壁面對(duì)熱流計(jì)算產(chǎn)生的影響。對(duì)我們研究的情況,只要兩個(gè)邊界間距離大于2.5倍中心處的分子自由程,即可獲得相對(duì)可靠的結(jié)果。

2.3 導(dǎo)熱系數(shù)修正曲線A k(Zh T)

我們計(jì)算了不同中心點(diǎn)溫度、不同溫度梯度的50個(gè)工況,ZhT在[0,0.4]范圍內(nèi)變化,所得的導(dǎo)熱系數(shù)修正曲線如圖7所示。圖中不同顏色符號(hào)代表了不同的中心點(diǎn)溫度,從1000 K到2000 K變化,我們可以看到,該曲線不依賴于當(dāng)?shù)氐臏囟?只要ZhT相同,所獲得的導(dǎo)熱系數(shù)修正值就相同。所有工況的導(dǎo)熱系數(shù)修正值A(chǔ)k隨著ZhT遵循統(tǒng)一的規(guī)律,呈單調(diào)遞減變化,也就是等效導(dǎo)熱系數(shù)比傳統(tǒng)的導(dǎo)熱系數(shù)偏小,傅里葉定律計(jì)算的熱流值偏大。

圖7 導(dǎo)熱系數(shù)修正值隨無量綱參數(shù)Zh T的變化Fig.7 Correctional coefficients of thermal conductivity with parameter Zh T

3 圓柱繞流

為初步驗(yàn)證上述規(guī)律是否可以有效提高對(duì)駐點(diǎn)熱流的預(yù)測精度,分別采用DSMC和考慮導(dǎo)熱系數(shù)修正的CFD模擬了單原子氣體氬氣繞半圓柱的流動(dòng)。CFD模型的計(jì)算框架及計(jì)算格式詳見文獻(xiàn)[4,5]。模擬算例為文獻(xiàn)[19]中Kn＝0.05的工況,圓柱半徑為0.1524 m,壁面溫度Tw＝500 K,來流溫度T∞＝200 K,來流氣體的數(shù)密度為8.494×1019/m3,來流馬赫數(shù)分別為Ma＝10和Ma＝25。

高超聲速繞圓柱流動(dòng)在駐點(diǎn)附近存在強(qiáng)的溫度梯度,流體沿圓柱表面向后運(yùn)動(dòng)時(shí)又存在強(qiáng)的剪切。因此,該問題需要同時(shí)考慮強(qiáng)剪切和強(qiáng)溫度梯度引起的稀薄氣體效應(yīng)的影響。在計(jì)算中,根據(jù)流場當(dāng)?shù)氐膮?shù)Zh和ZhT對(duì)導(dǎo)熱系數(shù)進(jìn)行了修正。2.3節(jié)中獲得了ZhT對(duì)導(dǎo)熱系數(shù)k的修正值A(chǔ)k(ZhT),Zh對(duì)導(dǎo)熱系數(shù)k的修正值A(chǔ)k(Zh)在文獻(xiàn)[2]中獲得。在此算例中,兩者以加權(quán)平均的方式獲得最終的導(dǎo)熱系數(shù)的修正值,即:

需要說明的是這里作為初步嘗試采用了加權(quán)平均的方式考慮參數(shù)Zh和ZhT對(duì)導(dǎo)熱系數(shù)的影響,主要基于流動(dòng)圖像,認(rèn)為在駐點(diǎn)附近應(yīng)是ZhT起主導(dǎo),沿圓柱表面向后運(yùn)動(dòng)的邊界層中應(yīng)為Zh起主導(dǎo)。更加合理的方式將在今后的工作中探討。

圖8和圖9分別給出了Ma＝10和Ma＝25時(shí)圓柱表面熱流系數(shù)的計(jì)算結(jié)果,同時(shí)給出了完全漫反射邊界條件計(jì)算的DSMC結(jié)果,以及使用滑移邊界條件的N-S方程和不使用滑移邊界條件的N-S方程結(jié)果。橫坐標(biāo)?為從駐點(diǎn)開始沿圓柱表面順時(shí)針旋轉(zhuǎn)的角度,縱坐標(biāo)為無量綱的熱流系數(shù)CH＝q/0.。從圖中可以看出,無論是否采用滑移邊界,傳統(tǒng)的N-S方程計(jì)算的熱流均高于DSMC結(jié)果。而通過本文所建立的導(dǎo)熱系數(shù)修正算法對(duì)不同馬赫數(shù)的圓柱表面熱流,尤其駐點(diǎn)附近熱流預(yù)測有了很大的改善。

圖8 高超聲速圓柱繞流表面熱流系數(shù)(Ma＝10)Fig.8 Surface heating coefficient for hypersonic flow over a cylinder(Ma＝10)

圖9 高超聲速圓柱繞流表面熱流系數(shù)(Ma＝25)Fig.9 Surface heating coefficient for hypersonic flow over a cylinder(Ma＝25)

4 結(jié) 論

本文以臨近空間高超聲速飛行器氣動(dòng)熱計(jì)算為背景,考慮具有一定稀薄程度的氣體在存在強(qiáng)溫度梯度時(shí)熱流的計(jì)算。主要獲得以下結(jié)論:

1)提出了表征該問題的稀薄效應(yīng)刻畫參數(shù)ZhT;

2)獲得了依賴參數(shù)ZhT的導(dǎo)熱系數(shù)修正規(guī)律;

3)通過導(dǎo)熱系數(shù)修正可有效提高高超聲速圓柱繞流壁面熱流的預(yù)測精度。

值得說明的是,本文僅考慮了一種簡單的情形,即單原子氣體不存在速度梯度和壓力梯度,僅存在強(qiáng)溫度梯度的情況下對(duì)熱流計(jì)算的影響。而實(shí)際飛行器問題則更為復(fù)雜,駐點(diǎn)附近可看作純的導(dǎo)熱問題,但是過了駐點(diǎn)后,法向還存在很大的速度梯度、壓力梯度,并且伴隨高溫的振動(dòng)能激發(fā),電離離解等效應(yīng),需要在今后的研究中逐步分析,以解決真實(shí)飛行器流動(dòng)中的問題。另外,飛行器的壁面條件應(yīng)如何處理,也是要專門研究的難題。

致謝：感謝周恒院士的指導(dǎo)和通過有啟發(fā)性的討論給予的幫助。