李憲開，張志雨，何淼生，朱斌鑌，柳 軍

(1. 國防科技大學 空天科學學院， 湖南 長沙 410073; 2. 沈陽飛機設計研究所揚州協(xié)同創(chuàng)新研究院有限公司， 江蘇 揚州 225000; 3. 上海交通大學 航空航天學院, 上海 200240)

臨近空間高超聲速飛行器氣動設計需要準確預測壁面熱流，工程經(jīng)驗及理論公式在近連續(xù)流區(qū)內逐漸失效[1-3]，數(shù)值模擬方法仍然是壁面熱流預測的主要工具。然而，當來流克努森數(shù)(Kn)大于0.01時，即便是簡單的高超聲速平板繞流問題仍然因方法的不同而存在較大差異[3-4]。

對于稀薄高超聲速流場，直接模擬Monte Carlo(direct simulation Monte Carlo, DSMC)方法被認為是最可能準確預測壁面熱流的方法[5]。該方法雖物理上滿足Boltzmann方程并對其進行直接模擬[6]，但并不能直接給出稀薄流場熱流輸運的數(shù)學方程。目前，基于宏觀方程的數(shù)值方法還無法在全流域范圍內獲得與DSMC結果一致的熱流結果，這說明高超聲速流動駐點熱流輸運中稀薄效應的作用機制還不夠清楚。

從工程應用的角度來看，納維-斯托克斯(Navier-Stokes，N-S)方程的應用是最廣泛的，其對熱流模擬的基礎是Fourier熱傳導定律。但是，將其直接應用于稀薄高超聲速流場會獲得偏高的熱流結果[2-3]，最直觀的原因是壁面速度滑移或溫度跳躍現(xiàn)象的出現(xiàn)，故而從20世紀80年代至今發(fā)展了諸多滑移或跳躍邊界條件來拓展N-S方程在稀薄流域的應用[7-9]。Maccormack[7,9]等提出的通用型滑移或跳躍邊條可使N-S方程在Kn=0.25的高超聲速圓柱繞流計算中獲得與DSMC誤差不大于10%的駐點熱流結果，但二者的流場結果具有極大的差異[10-11]。近期研究表明[12-14]，對于簡單的高超聲速圓柱繞流和平板繞流問題，即使采用更趨近于真實物理的高階滑移或跳躍邊界條件、非平衡邊界條件也無法獲得與DSMC一致的壁面附近流場信息。

從氣體物理的角度來看，隨著Kn的增大，局部流場的分子速度分布不再滿足平衡態(tài)Maxwell分布，也不再滿足線性本構關系和Fourier熱傳導關系。近年來，不少學者在拓展流體動力學方程上做了大量的研究，以其準確模擬稀薄流域流動問題，包括Eu-type方程[15-16]、Burnett-type方程[17]等。例如，Eu-type方程可以在Kn=0.25、Ma=10的圓柱繞流中獲得與DSMC誤差不大于10%的駐點熱流結果[15]。有趣的是，非線性駐點熱流輸運的相關研究[1-2]表明：基于Burnett方程的二階熱流展開項對非線性熱流的貢獻為正，即稀薄效應下的高階熱流表達應為Fourier熱流加上二階熱流項，然而卻發(fā)現(xiàn)DSMC結果低于N-S方程Fourier熱流，這一矛盾被歸因于展開項系數(shù)的不確定性。

此外，對于稀薄高超聲速流場的解析還可以采用基于速度分布函數(shù)計算的氣體動理學方法，包括氣體動理論統(tǒng)一算法[18]、統(tǒng)一氣體動理學格式[19]和多尺度粒子方法[20]等，這一類多尺度方法的研究也是近十幾年來稀薄計算的熱點。值得一提的是，這一類方法的發(fā)展通常采用DSMC方法結果作為驗證的標準[5]，本文內容屬機理討論，因此采用最可能準確的DSMC方法開展研究。

總結上述兩個角度的調研，我們不難發(fā)現(xiàn)熱流輸運中的稀薄效應表現(xiàn)在兩個方面：一是壁面速度滑移和溫度跳躍現(xiàn)象，使得真實熱流低于Fourier熱流[3,7-8]；二是非線性本構關系和非線性熱傳導，雖基于Burnett方程展開項分析獲得較Fourier熱流偏高的結果[2-3]，但諸多計算數(shù)據(jù)卻獲得非線性熱流低于Fourier熱流的觀點[2-3,10]。前述兩點稀薄效應的描述均基于宏觀上的理解和認識，在微觀層面上的理解較少，然而深入的理解對進一步發(fā)展準確的駐點熱流預測方法是至關重要的。因此，本文從Fourier熱傳導定律失效的觀點出發(fā)，基于DSMC方法對Kn為0.01～0.1、Ma為5～10的高超聲速圓柱繞流駐點熱流中的稀薄效應進行研究，旨在從微觀視角上給出連續(xù)方法高估駐點熱流的新理解。

1 問題描述

研究對象為半徑R=0.152 4 m的二維高超聲速圓柱繞流，以文獻[3]中Kn=0.05的算例為參照，通過改變來流壓力實現(xiàn)Kn為0.01～0.1的變化范圍。研究的物理問題如圖1所示，為規(guī)避熱力學非平衡帶來的影響，來流氣體選為氬氣(Ar)，相關的氣體參數(shù)如表1所示。

圖1 物理問題示意Fig.1 Schematic diagram of physical problem

表1 氬氣相關氣體參數(shù)

本文計算的算例分為兩組，第一組為不同來流Kn的算例，該組中Case-1至Case-10的來流Kn從0.01等差增長至0.1，公差為0.01；第二組為不同來流Ma的算例，該組中Case-1至Case-11的來流Ma從5等差增長至10，公差為0.5。第一組算例中來流Ma固定為10，第二組算例中來流Kn固定為0.05。

2 方法介紹

采用DSMC方法開展微觀視角下的數(shù)值研究工作，該方法是以唯象論為物理基礎，并借助以往稀薄氣體流動的模擬經(jīng)驗和數(shù)理統(tǒng)計知識而發(fā)展起來的一種計算算法，由Bird教授于20世紀70年代提出[6]。DSMC方法并不是直接對Boltzmann方程進行求解，而是從物理上模擬了一個與Boltzmann方程描述一致的氣體流動過程。20世紀90年代，Wagner[21]證明了DSMC方法的控制方程實為Boltzmann方程，自此之后DSMC方法才在高超聲速領域發(fā)揮重要作用。

DSMC方法的具體介紹本文不再贅述，詳細內容可參考Bird的專業(yè)書籍[6]。本文中采用的分子模型為可變徑硬球(variable hard sphere, VHS)模型，氣-壁相互作用過程通過Maxwell模型實現(xiàn)完全漫反射的固定壁溫邊界條件。對于VHS模型，可根據(jù)Chapman-Enskog理論和分子動力學理論獲得以下幾個平衡態(tài)參數(shù)的定義：

1)平衡態(tài)分子平均自由程為：

(1)

2)參考溫度下黏性系數(shù)可表述為：

(2)

3)黏性系數(shù)冪次關系式為：

(3)

4)熱傳導系數(shù)定義為：

(4)

式(1)～(4)中各參數(shù)的含義可參考Bird專業(yè)書籍[6]，這里不再贅述。

針對駐點熱流展開討論，采用三種不同的熱流計算方式。第一種是DSMC壁面熱流微觀統(tǒng)計結果，下文中采用“DSMC_data”來表示。DSMC的壁面熱流需要統(tǒng)計所有撞擊壁面的模擬分子的能量，如式(5)所示。

(5)

式中：qDSMC_data為壁面熱流密度；treal為真實物理時間；A為單位面積；e是單子氣體分子能量，包括平動能、轉動能和振動能，上標in和re分別代表了入射壁面和從壁面反射離開的模擬分子。

此外，為對比連續(xù)方法和微觀粒子算法在熱流結果上的差異，還采用了Fay-Riddell關系式作為第二種熱流表達方式，下文用“F-R correlation”表示，該關系式可寫成：

(6)

式中，各參數(shù)含義可參考文獻[22]，本節(jié)不再贅述。

第三種熱流表達方式耦合了DSMC流場計算結果和Fourier熱傳導定律，即基于DSMC流場溫度的Fourier熱流，下文用“DSMC_Fourier”來表示，具體表達式為：

(7)

式中：κ為熱傳導系數(shù)，可由式(4)求得；η為壁面法向方向。

最后，為使不同來流條件下的熱流具有可比性，采用來流參數(shù)對其進行無量綱化：

(8)

3 結果與討論

3.1 不同稀薄效應下駐點熱流規(guī)律

圖2所示為不同來流Kn下駐點熱流變化規(guī)律，橫坐標為來流Kn的1/2次方，縱坐標為熱流系數(shù)cq。其中：“F-R correlation”通過式(6)求得；“DSMC_data”為DSMC壁面熱流結果，由式(5)確定；“DSMC_Fourier”為基于DSMC流場溫度的Fourier熱流結果，由式(7)求得。此外，圖2中采用藍色十字圖標標出了文獻[3]中的無滑移N-S方程的結果，與紅色實線的匹配結果可驗證本文F-R關系式的正確性。

圖2 不同表達方式駐點熱流隨來流Kn的變化規(guī)律Fig.2 Trend of stagnation point heat fluxes with different Kn

根據(jù)式(6)的推導，容易得到F-R關系式熱流結果與Kn1/2成正比，圖中“F-R correlation”曲線所示亦如此。黑色虛線所示的DSMC壁面熱流結果從Kn1/2>0.1(Kn>0.01)開始逐漸偏離F-R關系式預測結果，偏離程度隨Kn增大而增大，并且不再保持與Kn1/2的線性關系。前述觀察到的現(xiàn)象具有共識性，基于對稀薄效應的宏觀認識(壁面滑移或跳躍現(xiàn)象)可以較好地理解。

有趣的是，黑色實線所描述的基于DSMC流場溫度的Fourier熱流結果“DSMC_Fourier”與DSMC壁面熱流結果存在差異。這意味著即便N-S方程等連續(xù)方法能夠借助于壁面滑移或跳躍邊界條件獲得與DSMC結果一致的宏觀流場，也不能得到準確的壁面熱流結果。更詳細地，不同Kn下“DSMC_Fourier”駐點熱流較DSMC微觀統(tǒng)計結果高，且相對誤差隨著Kn的增加逐漸從20%減小至10%左右，這隱喻著Fourier熱傳導定律的失效。

通常，稀薄效應由來流密度變小或研究對象變小所引起，可以理解為流場狀態(tài)的稀薄程度，如圖2所示的Kn變化。來流Ma增大，會使激波及壁面附近的梯度增大，進而局部Kn隨之增大，可以理解為流動的稀薄程度。圖3給出了不同來流Ma時，三種不同熱流表達方式下的駐點熱流變化規(guī)律?？梢钥吹?，由來流Ma引起的稀薄效應增強同樣使得駐點熱流結果增大，但三種不同表達方式下的熱流增長具有差異性。黑色虛線所示的DSMC壁面熱流結果均小于F-R關系式預測結果，且二者的差異隨來流Ma的增加而增大。此外，“DSMC_Fourier”所示的Fourier結果在當前Ma范圍內仍然高估駐點熱流，且相對誤差隨著Ma的增加逐漸從24%減小至12%左右。

圖3 不同表達方式駐點熱流隨來流Ma的變化規(guī)律Fig.3 Trend of stagnation point heat fluxes with different Ma

總結圖2與圖3所示的熱流差異規(guī)律，可以做如下理解：

1) F-R關系式等連續(xù)方法與DSMC微觀統(tǒng)計的駐點熱流之間的差異主要由溫度跳躍現(xiàn)象引起，這是稀薄效應的第一種體現(xiàn)；

2) 基于DSMC宏觀溫度的Fourier熱流和DSMC微觀統(tǒng)計的駐點熱流之間的差異表明，F(xiàn)ourier熱傳導定律不再適用于稀薄高超聲速駐點熱流的預測，這是稀薄效應的第二種體現(xiàn)。

上述兩點理解通過兩種宏觀現(xiàn)象闡述，還需要基于微觀視角進行進一步的解釋，主要回答兩個基礎問題：一是溫度跳躍現(xiàn)象的微觀理解是什么？二是稀薄效應下Fourier熱傳導定律的失效機制是什么？后續(xù)兩個小節(jié)將分別針對前述兩個基礎問題展開。

3.2 駐點溫度跳躍現(xiàn)象的微觀理解

通常認為，F(xiàn)-R關系式及N-S方程等連續(xù)方法在稀薄條件下高估熱流的主要原因是未考慮壁面滑移或跳躍現(xiàn)象。如同設計溫度跳躍邊界條件一樣，這是基于宏觀的一種認識，通常借助局部克努森數(shù)來設計。局部克努森數(shù)可定義[8]為：

(9)

式中，Q是溫度、速度和壓強等流場宏觀參數(shù)，η為壁面法向方向。

此外，一種簡單且常用的溫度跳躍邊界條件可以是Smoluchowski跳躍邊界條件[7,9]：

(10)

若采用溫度計算KnGLL并代入跳躍邊界條件，可以得到：

(11)

由于壁面附近溫度變化不算太大，式(11)中γ和Pr可作常數(shù)處理。因此可認為在當前計算范圍內，式(11)所示的比值基本保持不變。

式(10)或式(11)為溫度跳躍的數(shù)學或宏觀表達式，為印證其合理性，統(tǒng)計了不同來流Kn和來流Ma下式(11)左側所示的比值變化規(guī)律，如圖4中帶三角符號的紅色直線所示?？梢钥吹剑摫戎低琄n和Ma基本保持無關性，這表明在當前考慮的稀薄效應程度下，采用KnGLL衡量連續(xù)介質假設失效下的溫度跳躍現(xiàn)象是合適的，文獻[10]也表明采用式(10)所示的溫度跳躍邊界條件可以在Ma=10、Kn=0.05條件下獲得與DSMC誤差不超過10%的駐點熱流結果。

(a) 來流Kn變化(a) Freestream Kn changes

在微觀層面，當邊界層內、外存在足夠大的溫差時可以認為溫度跳躍現(xiàn)象必然存在。無論第一層網(wǎng)格多小，其內所含模擬分子不可能全部由壁面反射分子組成，第一層網(wǎng)格的物理量變化更多地取決于其他網(wǎng)格和該網(wǎng)格之間的流通量。故而，第一層網(wǎng)格宏觀溫度必然大于壁面溫度，而這種跳躍溫差取決于分子間碰撞頻率的大小。當流動逐漸偏離平衡態(tài)時，分子間碰撞頻率會減小，可能入射的模擬分子群和壁面反射分子群之間的碰撞減少，導致壁面附近流體溫度越發(fā)高于壁面溫度，溫度跳躍現(xiàn)象增強?；谏鲜隼斫猓x一個微觀統(tǒng)計量非平衡度Dn來表征當?shù)亓鲌龅姆瞧胶舛龋?/p>

(12)

式中，νsimu為DSMC碰撞過程中實際的碰撞頻率，νtheo為VHS模型對應的理論平衡態(tài)碰撞頻率。那么，Dn=0為平衡態(tài)，Dn>0或Dn<0為非平衡態(tài)。值得注意的是，采用的氬氣為單原子分子，因此這里的非平衡指的是分子速度分布偏離Maxwell平衡態(tài)分布，也可稱為平動非平衡。

為驗證前述所提出的偏離非平衡態(tài)和溫度跳躍之間的關系，圖4中列出了不確定度Dn、溫度跳躍Tj和來流Kn、來流Ma之間的關系。為保證Dn和Tj在同一坐標中具有可比性，圖中藍色虛線所示為溫度跳躍的萬分之一倍(Tj/104)。從圖4中不難發(fā)現(xiàn)，Dn和Tj均與Kn和Ma呈一次線性正相關關系。此外，圖中Dn與Tj曲線的斜率是十分接近的，可認為非平衡度在一定程度上可以用于描述溫度跳躍的大小，這也表明了溫度跳躍現(xiàn)象在微觀層面上是由流動偏離平衡態(tài)所引起。

3.3 Fourier熱傳導定律的失效機制

回顧圖2和圖3中所示的熱流規(guī)律可知，基于DSMC流場溫度的Fourier熱流結果與DSMC壁面熱流仍然存在差異，這表明稀薄效應增強時，采用Fourier定律描述沿駐點線的熱流輸運是不合適的。通常認為該問題由線性本構關系失效導致，進而發(fā)展了諸多高階熱流關系式來解決該問題。然而，類如Wang、Singh等的研究工作表明二階熱流項的貢獻為正[1-2]，體現(xiàn)在壁面熱流上應該獲得較Fourier熱流更大的結果，但數(shù)值結果獲得相反的結論。

本節(jié)通過對比DSMC和基于DSMC流場溫度的Fourier表達沿駐點線的熱流結果，嘗試給出線性Fourier熱傳導定律的微觀失效機制。圖5所示為Kn=0.01、Ma=10的高超聲速圓柱繞流駐點線熱流輸運結果，其中“DSMC_Fourier”為基于DSMC流場溫度的Fourier熱流，由式(7)計算求得；“DSMC_data”為DSMC微觀統(tǒng)計結果，通過對網(wǎng)格內模擬分子的能量進行統(tǒng)計而寫成熱流通量的形式[23]。

(13)

式中：c為矢量熱運動速度；int為單個氣體分子的內能，包括轉動能及振動能。當氣體為單原子時，略去式(13)第二項。

圖5 駐點線熱流輸運及非平衡度變化趨勢Fig.5 Variation of the heat transfer and non-equilibrium degree along the stagnation streamline

根據(jù)熱運動速度分解，可寫出x方向上的熱流通量表達式[23]：

(14)

此外，由于速度分布偏離平衡態(tài)分布是線性本構關系不再成立的內在原因，因此圖5中還給出了沿駐點線的非平衡度Dn的變化趨勢，如圖中藍色曲線所示。

圖5中，綠色實線所示的“DSMC_Fourier”線性熱流輸運在激波內和壁面附近與DSMC微觀統(tǒng)計結果存在明顯差異。另外，藍色曲線所示的非平衡度在激波和壁面附近均大于零(分子速度分布偏離平衡態(tài))，可認為局部稀薄誘導的非平衡效應是Fourier熱傳導定律不再適用的本質原因。更詳細地，線性Fourier表達在激波內低估傳熱，在壁面附近高估傳熱。激波內大梯度引起的偏離平衡態(tài)分布和壁面附近大梯度引起的偏離平衡態(tài)分布在分布形式上是不同的，盡管兩種偏離均引起碰撞頻率的減小，但最終引起的線性熱傳導失效是不同的。另外，還可以觀察到一個有趣的現(xiàn)象：非平衡度的峰值并沒有對應著熱流最大的位置，而是比較靠近激波的前邊界，也許這也是諸多可計算氣體動理學方法在激波前緣位置與DSMC結果吻合不太好的緣故。

圖5右上角的紅色實心圓點為DSMC駐點熱流結果，該結果高于壁面第一層網(wǎng)格內的DSMC微觀統(tǒng)計熱流。對比了所有算例的第一層網(wǎng)格熱流和駐點熱流，結果表明：駐點處第一層網(wǎng)格熱流均小于駐點熱流。并且，該差異隨Kn和Ma的增大而增大，在本文計算的最小Kn和最小Ma條件下，二者是十分接近的。第一層網(wǎng)格熱流和壁面熱流之間的差異正是溫度跳躍存在的一種證明。

此外，綠色曲線所示的“DSMC_Fourier”熱流輸運在靠近壁面幾層網(wǎng)格內出現(xiàn)了突增的現(xiàn)象，本文所計算的算例均出現(xiàn)了不同程度的Fourier熱流突增現(xiàn)象，并且這種突增現(xiàn)象均發(fā)生在距離壁面約3倍當?shù)胤肿悠骄杂沙痰目臻g內。正好，根據(jù)Maxwell速度分布可知，在一個平均碰撞時間內95%以上的分子可自由移動的距離均在3倍分子平均自由程內。同理，從壁面反射的分子在無分子間碰撞的條件下可自由移動的距離基本上均小于3倍分子平均自由程，即壁面反射分子直接影響流場的作用范圍約為3倍分子平均自由程。因此，將觀察到的Fourier熱流突增現(xiàn)象解釋為壁面約束效應，這和文獻[24]中觀察到的近壁面分子平均自由程變化規(guī)律在物理上是相通的。

為量化前述壁面約束效應，表2列出了各算例條件下的Fourier熱流突增的強度，表中的數(shù)值為近壁面3倍分子平均自由程范圍內的Fourier熱流最大值和最小值之比。從表2中可看到，量化Fourier熱流壁面約束效應的比值隨Kn的增大而減小，也隨Ma的增大而減小。因此，來流Kn和Ma的增大均使得Fourier熱流突增減弱?；仡檲D5中“DSMC_Fourier”曲線，壁面附近平動非平衡效應使Fourier表達高估熱流，壁面約束效應同樣使Fourier表達高估熱流。圖4和表2表明，壁面附近平動非平衡度隨Kn、Ma增加而增加，壁面約束效應隨Kn、Ma增加而減弱。二者共同作用下，F(xiàn)ourier表達下的駐點熱流高于DSMC結果，且相對偏差隨Kn、Ma增加而減小，從而可認為壁面約束效應對Fourier熱傳導失效的貢獻更大。

表2 壁面約束效應下Fourier熱流突增強度規(guī)律

由于DSMC方法的特點是隨機，雖獲得了流場的非平衡度描述，但無法確定哪些模擬分子代表非平衡態(tài)，哪些模擬分子代表平衡態(tài)。也正由于模擬分子的隨機碰撞和自由移動，可在統(tǒng)計壁面熱流時以非平衡度的大小按概率統(tǒng)計平衡態(tài)熱流和非平衡態(tài)熱流。若非平衡度為0.1，那么與壁面發(fā)生碰撞的模擬分子有10%的能量交換統(tǒng)計為非平衡態(tài)熱流，剩余90%的能量交換統(tǒng)計為平衡態(tài)熱流。

圖6 非平衡態(tài)熱流隨來流Kn、Ma的變化規(guī)律Fig.6 Variation of the nonequilibrium heat flux with respect to the freestream Kn and Ma

基于上述觀點，可獲得如圖6所示的非平衡態(tài)熱流隨Kn和Ma的變化規(guī)律，且容易看到：非平衡態(tài)熱流與Kn和Ma均呈現(xiàn)出線性正相關關系。同時，由于Kn為零時，流動處于平衡態(tài)，此時非平衡態(tài)熱流也應為零，因此可以認為非平衡態(tài)熱流與Kn成正比。

4 結論

針對工程上關心的稀薄高超聲速駐點熱流預測中存在的方法適用性問題，采用Fay-Riddell關系式、DSMC方法和基于DSMC流場溫度的Fourier傳熱三種熱流表達方式，對不同來流Kn和Ma的高超聲速圓柱繞流駐點熱流規(guī)律及差異開展了研究，針對經(jīng)典連續(xù)方法在稀薄流區(qū)高估駐點熱流的現(xiàn)象給出了如下幾點理解：

1) 基于微觀統(tǒng)計的非平衡度參數(shù)Dn和壁面溫度跳躍Tj均與Kn、Ma呈線性正相關關系且斜率接近，Dn所描述的平動非平衡正是溫度跳躍現(xiàn)象在微觀層面上的誘因。溫度跳躍會削弱溫度梯度導致熱流降低，從而無滑移連續(xù)方法會高估駐點熱流。

2) 流場局部稀薄效應使得壁面附近分子速度分布偏離平衡態(tài)，線性本構關系不再保持，F(xiàn)ourier熱傳導定律失效且高估熱流，從而基于DSMC流場溫度的Fourier熱流高估駐點熱流。

3) 受到壁面反射分子的影響，基于DSMC流場溫度的Fourier熱流在壁面附近約3倍分子平均自由程內高估熱流，稱之為壁面約束效應，這也是Fourier熱傳導定律失效且高估駐點熱流的第二種體現(xiàn)。