高超聲速計(jì)算中的氣體動(dòng)理學(xué)格式

徐昆， 陳松澤

香港科技大學(xué) 數(shù)學(xué)系， 香港

高超聲速計(jì)算中的氣體動(dòng)理學(xué)格式

徐昆*， 陳松澤

香港科技大學(xué) 數(shù)學(xué)系， 香港

回顧了高超聲速連續(xù)流部分的計(jì)算流體力學(xué)(CFD)方法，總結(jié)了近些年興起的氣體動(dòng)理學(xué)格式。闡述了該格式的構(gòu)造機(jī)制，強(qiáng)調(diào)了將物理規(guī)律直接用于構(gòu)造數(shù)值方法的思路。結(jié)合一些應(yīng)用實(shí)例，例如激波相互作用、激波邊界層相互作用以及邊界層分離等高超聲速問題，說明了這種構(gòu)造思路給數(shù)值模擬帶來的優(yōu)點(diǎn)。從高超聲速的發(fā)展歷程來看，氣體動(dòng)理學(xué)格式的構(gòu)造過程包含了更基礎(chǔ)的物理規(guī)律，而且具有多尺度的特性。這些特性有助于研究復(fù)雜的高超聲速問題。介觀或者微觀角度直接構(gòu)造數(shù)值方法的發(fā)展趨勢為高超聲速計(jì)算工具指出了可能的發(fā)展方向。

高超聲速； 氣體動(dòng)理學(xué)格式； 滑移邊界條件； 激波相互作用； 激波邊界層干擾

通常將來流馬赫數(shù)大于5的流動(dòng)稱為高超聲速流動(dòng)。20世紀(jì)40年代，中國著名科學(xué)家錢學(xué)森首先提出了高超聲速的概念[1]。20世紀(jì)60年代，美國進(jìn)行的2次X15飛行試驗(yàn)就是高超聲速流動(dòng)的直觀例子[2]。其中一次試驗(yàn)中飛行馬赫數(shù)達(dá)到了4.9，飛行器沒有出現(xiàn)明顯的燒蝕現(xiàn)象；而另一次試驗(yàn)中飛行馬赫數(shù)達(dá)到了6.9，飛行器燒蝕非常嚴(yán)重。前者可以看作普通的超聲速飛行，而后者反映了高超聲速流動(dòng)中的一個(gè)重要特點(diǎn)——高能高溫。由于來流速度非常大，高速氣體靠近固體減速時(shí)會(huì)將大量的動(dòng)能轉(zhuǎn)化為內(nèi)能，從而使得流場中出現(xiàn)極高的壓力和溫度。

隨著馬赫數(shù)的提高，高溫高壓會(huì)逐漸造成許多非平衡效應(yīng)。例如，氣體分子平動(dòng)和轉(zhuǎn)動(dòng)能量的非平衡、振動(dòng)能量激發(fā)、氣體分子離解、原子電離以及高溫氣體輻射效應(yīng)等。而且由于高超聲速飛行器的飛行環(huán)境往往是稀薄大氣，所有的這些非平衡還需要疊加上稀薄效應(yīng)?？梢哉f高超聲速學(xué)科是一個(gè)內(nèi)容異常豐富的交叉學(xué)科。本文不涉及以上的非平衡問題，而主要關(guān)注高超聲速中連續(xù)流和近連續(xù)流的氣動(dòng)力和氣動(dòng)熱計(jì)算相關(guān)問題。

20世紀(jì)70年代以前由于計(jì)算機(jī)和數(shù)值算法限制，只有很小一部分問題能夠用計(jì)算方法研究。這個(gè)時(shí)期是高超聲速空氣動(dòng)力學(xué)理論研究的黃金時(shí)期，牛頓流模型、高超聲速無黏小擾動(dòng)理論、可壓縮邊界層理論等都得到了長足的發(fā)展。然而，面對復(fù)雜外形或者復(fù)雜流動(dòng)時(shí)，理論分析方法就顯得力不從心。

熱防護(hù)問題是設(shè)計(jì)空天飛行器的一個(gè)重要問題。一般來說，高溫對飛行安全是有害的。必須考慮如何通過優(yōu)化設(shè)計(jì)來減小熱流峰值以及如何布置熱防護(hù)層?，F(xiàn)在的高超聲速飛行器大量使用鈍體外形。如今我們知道駐點(diǎn)熱流率和鈍頭半徑的平方根成反比。然而，在20世紀(jì)40年代到50年代初，人們一直嘗試用帶有尖端的細(xì)長體做高超聲速飛行器[3]。當(dāng)時(shí)人們對氣動(dòng)加熱的影響估計(jì)不足，希望直接沿用普通超聲速飛行器的設(shè)計(jì)經(jīng)驗(yàn)，走了不少彎路。后來，經(jīng)過了大量的理論研究和實(shí)驗(yàn)后，人們才認(rèn)識到鈍體外形在高超聲速飛行器設(shè)計(jì)中的重要地位。這種設(shè)計(jì)上的優(yōu)化往往只能緩解嚴(yán)重的氣動(dòng)加熱問題，而布置熱防護(hù)層則是和氣動(dòng)加熱問題直接對抗。簡單外形的高超聲速繞流大多數(shù)是附著流動(dòng)，流場結(jié)構(gòu)簡單。駐點(diǎn)處的熱流往往是最大的。這種情況下，只要在相應(yīng)的位置布置合適的熱防護(hù)層就能較好地解決氣動(dòng)熱問題。對于較為復(fù)雜的外形，即使是雙錐外形，由于激波相互作用、激波邊界層相互作用、黏性相互作用以及湍流等原因，高超聲速的流場結(jié)構(gòu)變得異常復(fù)雜。不同的流動(dòng)結(jié)構(gòu)會(huì)產(chǎn)生各種嚴(yán)重的局部氣動(dòng)熱問題。熱防護(hù)層的布置就需要預(yù)先了解流場結(jié)構(gòu)和飛行器表面的熱流分布。對于稍復(fù)雜一些的外形，理論方法就難于勝任了。從這個(gè)例子可以看出，由于人們以往經(jīng)驗(yàn)和知識的片面性和局限性，往往很難通過理論方法直接得到正確的結(jié)論。正如Rasmussen[4]在他的書中提到“在現(xiàn)代條件下，這些(分析方法)努力所得的回報(bào)越來越少，還不如直接通過計(jì)算機(jī)求解控制方程”。

高超聲速學(xué)科從誕生開始就和工程應(yīng)用緊緊地聯(lián)系在一起。從美國的X系列飛行器到獵鷹HTV-2(Hypersonic Technology Vehicle 2)滑翔式高超聲速飛行器，工程應(yīng)用在高超聲速學(xué)科發(fā)展中一直扮演著重要角色。這些工程應(yīng)用給高超聲速學(xué)科提供了大量的試驗(yàn)數(shù)據(jù)。高超聲速試驗(yàn)可以分為飛行試驗(yàn)和地面試驗(yàn)。飛行試驗(yàn)自然可以獲得第一手資料，但是其成本高昂且試驗(yàn)難度大。地面試驗(yàn)一般在各種類型的高超聲速風(fēng)洞中進(jìn)行。由于實(shí)際飛行中的高超聲速流動(dòng)條件非常復(fù)雜，控制參數(shù)多，地面試驗(yàn)往往只能滿足部分流動(dòng)條件。例如，低湍流度來流、稀薄條件、高焓條件以及由此帶來的一系列效應(yīng)等，許多地面試驗(yàn)設(shè)備無法很好地模擬。這就造成地面試驗(yàn)得到的流動(dòng)狀態(tài)常常和正式飛行中出現(xiàn)的現(xiàn)象相左。以X15為例，風(fēng)洞試驗(yàn)數(shù)據(jù)顯示流場基本全是層流，而飛行測試表明流動(dòng)基本是湍流[5]。這種差異使得地面試驗(yàn)的價(jià)值大打折扣。事實(shí)上，許多威脅到飛行安全的嚴(yán)重的氣動(dòng)加熱問題都是在實(shí)際高超聲速飛行中才被人們發(fā)現(xiàn)的[2]。

高超聲速的流場復(fù)雜，使得從理論框架來考慮各種因素的影響變得不太現(xiàn)實(shí)。而地面試驗(yàn)無法完全重現(xiàn)飛行條件。這就需要一種工具來將地面試驗(yàn)的結(jié)果推廣到真實(shí)的飛行條件。計(jì)算流體力學(xué)(CFD)在高超聲速工程應(yīng)用中正好扮演了這種角色。人們通過地面試驗(yàn)獲得各種流場數(shù)據(jù)，再將這些數(shù)據(jù)作為標(biāo)準(zhǔn)來檢驗(yàn)計(jì)算模擬程序。然后，通過驗(yàn)證后的計(jì)算程序再用來計(jì)算真實(shí)條件下的流動(dòng)狀態(tài)。以往的經(jīng)驗(yàn)說明從計(jì)算中可以獲得大量并且詳細(xì)的信息。這些信息能夠補(bǔ)充和完善實(shí)驗(yàn)和理論分析的結(jié)果，幫助人們更好地理解基本的物理過程。現(xiàn)代的高超聲速問題要求人們更多地借助于計(jì)算機(jī)模擬等手段進(jìn)行探索和研究。

1 簡述計(jì)算流體力學(xué)的發(fā)展

1.1 捕捉間斷

計(jì)算流體力學(xué)發(fā)展早期，由于計(jì)算機(jī)技術(shù)和數(shù)值方法的限制，只有很小一部分問題能用計(jì)算模擬得到較好的結(jié)果。1966年，Moretti和Abbett[6]針對鈍體繞流問題采用有限差分時(shí)間推進(jìn)格式求解了的非定常的Euler方程，并且得到了漸近的穩(wěn)態(tài)流場。這是計(jì)算流體力學(xué)的重大突破。他們的工作對計(jì)算流體力學(xué)有兩大貢獻(xiàn)[7]。

1)時(shí)間推進(jìn)方式能夠使得非定常Euler方程無論在超聲速還是在亞聲速流動(dòng)區(qū)都保持雙曲性。

2)證明了有限差分方法捕捉能夠捕捉到Rankine-Hugoniot激波間斷條件。

1969年，MacCormack[8]發(fā)表了一篇里程碑式的論文，論文介紹了一種顯式數(shù)值格式。這種數(shù)值格式后來以他的名字命名。這種2階精度算法編程簡潔且穩(wěn)定性好，對非常大的梯度流場問題也適用。MacCormack顯式格式已經(jīng)作為計(jì)算流體力學(xué)的關(guān)鍵研究工具被廣泛使用了許多年。這兩項(xiàng)工作開啟了計(jì)算流體力學(xué)的高超聲速應(yīng)用的大門。在20世紀(jì)70年代，三維高超聲速壓縮拐角問題已經(jīng)可以用MacCormack顯式格式模擬[9]。

20世紀(jì)80年代以前，計(jì)算流體力學(xué)格式以有限差分格式居多。此后，隨著非結(jié)構(gòu)網(wǎng)格出現(xiàn)有限體積格式才真正流行起來[7]。許多復(fù)雜工程問題都能方便地使用有限體積格式求解。這得益于有限體積法中蘊(yùn)藏的守恒律的思想。宏觀守恒律是流體運(yùn)動(dòng)的基本特征。無論是推導(dǎo)宏觀方程，還是構(gòu)造數(shù)值格式，都需要將這種物理規(guī)律體現(xiàn)出來。守恒律的概念實(shí)際上可以定義在任何控制體上?？刂企w內(nèi)的各種因變量，不論是質(zhì)量、動(dòng)量的各個(gè)分量還是內(nèi)能，都可以由通過控制體表面的流量寫出平衡關(guān)系。有限體積法將這種關(guān)系應(yīng)用于每個(gè)網(wǎng)格，嚴(yán)格保證每個(gè)網(wǎng)格單元內(nèi)以及整個(gè)計(jì)算區(qū)域內(nèi)的守恒關(guān)系。相比于有限差分格式，有限體積格式更不容易受到網(wǎng)格幾何奇異性的影響，在復(fù)雜的流動(dòng)應(yīng)用問題中有更有優(yōu)異的表現(xiàn)。當(dāng)然有限差分法也通過發(fā)展守恒型差分格式體現(xiàn)了守恒律的重要性。文獻(xiàn)[7]指出人們不能期望在真實(shí)解是不可微的地方通過求解微分方程從而找到一種合理的解。因此，激波捕捉方法必須建立在積分形式基礎(chǔ)之上而不是微分的守恒律,而有限體積法也正好契合了這種思想。

有限體積格式的基本公式需要重構(gòu)垂直于網(wǎng)格邊界的通量。其大部分的工作是關(guān)于如何在離散的空間中描述間斷的流體現(xiàn)象“黎曼問題”。Godunov[10]發(fā)表了一篇標(biāo)志性的論文，給出了能夠預(yù)測包含激波的多維流場和接觸面問題方法。他假設(shè)在控制體內(nèi)數(shù)據(jù)符合連續(xù)線性分布，在網(wǎng)格邊界處通過求解黎曼問題來處理雙曲系統(tǒng)的間斷問題。通過在整個(gè)計(jì)算區(qū)域上求解一組黎曼問題，這個(gè)方法能夠根據(jù)不同波的傳播方向正確的給出物理量的依賴區(qū)。Steger和Warming[11]提出的通量分裂格式也是基于特征理論的。他們系統(tǒng)地給出了實(shí)特征量和特征向量以及分裂通量的公式。通量分裂格式已經(jīng)被證明能高效捕捉和解析激波，但是在聲速點(diǎn)和駐點(diǎn)會(huì)有一些瑕疵。此后，由van Leer[12]和Roe[13-14]開發(fā)的通量分裂格式廣泛地用于解析激波的計(jì)算流體力學(xué)應(yīng)用中。他們掀起了應(yīng)用近似黎曼求解器改進(jìn)迎風(fēng)格式計(jì)算對流通量的熱潮，出現(xiàn)了一批優(yōu)秀的近似黎曼求解器，如HLLE (Harten-Lax-van Leer-Einfeldt)[15]和HLLC (Harten-Lax-Van Leer-Contact)[16]等。

為了消除或減小激波附近的振蕩現(xiàn)象，Harten[17]提出了TVD (Total Variation Diminishing)格式的概念。這種格式實(shí)現(xiàn)了數(shù)值解的總變差不增，從而消除了非物理的振蕩現(xiàn)象。此后，中國的張涵信院士提出NND格式(Non-oscillatory containing No free parameter and dissipative Difference scheme)[18]。該格式從物理問題滿足熵條件出發(fā)也具有總變差不增的性質(zhì)。NND格式在國內(nèi)獲得了廣泛的應(yīng)用。

可以看到，守恒律、特征理論以及TVD性質(zhì)(熵條件)，都可以直接應(yīng)用在數(shù)值格式中，如有限體積法、Godunov格式和NND格式。同時(shí)也可以通過求解控制方程這種間接的方式來實(shí)現(xiàn)，如守恒型格式、通量分裂格式和TVD格式。也許這些方法細(xì)節(jié)有些許不同，但都試圖將某種物理規(guī)律或者物理考慮應(yīng)用到數(shù)值方法中。從這些方法的發(fā)展來看，最精確和高效的求解問題的數(shù)值方法是那些能最好模擬物理的方法。

1.2 提高精度

以上的工作都是關(guān)于如何求解有間斷物理量的流場問題。在科研人員的努力之下，現(xiàn)在的數(shù)值方法在捕捉間斷方面已經(jīng)比較成熟，能夠處理高超聲速中大多數(shù)的激波問題。雖然有這些顯著的成就，數(shù)值算法的分辨率依然還有許多改進(jìn)的空間。例如，TVD格式在極值點(diǎn)附近會(huì)降階。這就影響了格式在光滑流場區(qū)域的使用。因此，后來發(fā)展出來的ENO(Essentially Non-Oscillatory)格式弱化了TVD格式性質(zhì)，從而在極值點(diǎn)附近保持了格式的精度[19]。ENO格式選擇最光滑的重構(gòu)模板，而WENO(Weighted ENO)格式將ENO中考慮的模板進(jìn)行加權(quán)組合實(shí)現(xiàn)了光滑區(qū)更高階的精度[20]。這種高精度高分辨率的需求主要來自于黏性流模擬，特別是湍流模擬。1992年Lele[21]提出了緊致有限差分格式,這種格式很快在領(lǐng)域內(nèi)流行起來。該格式空間分辨率和譜方法相近。

緊致格式首先被用于不可壓或弱可壓縮流動(dòng)。當(dāng)被應(yīng)用于有激波的問題時(shí)，緊致格式會(huì)在間斷附近產(chǎn)生非物理的振蕩。這個(gè)問題可以通過避免在激波附近使用緊致格式來解決。具體來說可以使用激波裝配法[22-24]或者在間斷附近使用激波捕捉格式如ENO或WENO[25-26]。這些格式常常被用于高超聲速邊界層轉(zhuǎn)捩和可壓縮湍流研究。

無論是激波問題、黏性問題、湍流問題還是它們之間的相互作用，都是高超聲速流動(dòng)中的關(guān)鍵流動(dòng)現(xiàn)象。計(jì)算流體力學(xué)多年的發(fā)展可以說一直都是圍繞著這些問題展開的。

1.3 控制方程的變化

現(xiàn)在數(shù)值模擬中使用的控制方程多為Euler方程和Navier-Stokes方程。但高超聲速問題中黏性作用非常突出。比如，高超聲速邊界層厚度不僅和雷諾數(shù)相關(guān)還與馬赫數(shù)的平方成正比。許多條件下，高超聲速流場的激波層中都不存在無黏區(qū)域。因此，高超聲速模擬中黏性問題就更為重要。高超聲速模擬中求解黏性流方程常常用到以下3種方程，黏性激波層(Viscous Shock-Layer, VSL)方程[27]、拋物化(Parabolized)的Navier-Stokes(PNS)方程[28-29]、完整的Navier-Stokes方程。VSL方程在經(jīng)典的邊界層方程基礎(chǔ)之上，弱化了壓力穿透條件，并考慮了物面法向的梯度。PNS方程比VSL方程更進(jìn)一步，只忽略了流向梯度的黏性效應(yīng)。這2種方程都是Navier-Stokes方程的某種近似。它們在應(yīng)用中有諸多限制。例如，VSL方程不能計(jì)算分離流動(dòng)，PNS方程只能計(jì)算橫向分離流動(dòng)。盡管如此，在高馬赫數(shù)低雷諾數(shù)時(shí)，PNS方程也比無黏方法得到的結(jié)果好。人們選擇使用近似的VSL方程和PNS方程主要是為了減小計(jì)算消耗。在計(jì)算機(jī)高速發(fā)展的今天，計(jì)算資源對計(jì)算模擬的限制越來越小，諸如VSL方程和PNS方程的應(yīng)用空間也相應(yīng)變小了。相比于計(jì)算資源的消耗，人們更加關(guān)心方程是否能精確地反映物理現(xiàn)象，模擬是否精確?？梢钥吹剑S著計(jì)算機(jī)技術(shù)的不斷發(fā)展，計(jì)算模擬使用的控制方程將會(huì)包含更豐富的物理內(nèi)容以便更精確地描述物理現(xiàn)象。

2 氣體動(dòng)理學(xué)格式

氣體動(dòng)理學(xué)格式(Gas Kinetic Scheme, GKS)[30-32]是一種求解Navier-Stokes方程的有限體積格式。這種格式與Godunov格式類似，通過求解網(wǎng)格邊界處的黎曼問題來實(shí)現(xiàn)對物理過程的模擬。但和Godunov格式不同的是GKS使用的物理規(guī)律更加基礎(chǔ)，它使用的是玻爾茲曼(Boltzmann)方程。

2.1 GKS的構(gòu)造思路

2.1.1 更基礎(chǔ)的物理規(guī)律

人們使用的控制方程逐漸向完整的Navier-Stokes方程過渡。這種趨勢延續(xù)下去，自然就需要使用比Navier-Stokes方程更為基礎(chǔ)的氣體控制方程Boltzmann方程。Euler方程和Navier-Stokes方程描述宏觀物理量的演化規(guī)律，而Boltzmann方程描述的是微觀氣體分子的分布函數(shù)的演化規(guī)律。根據(jù)氣體動(dòng)理學(xué)理論可以通過Boltzmann方程的Chapman-Enskog展開推導(dǎo)得到Navier-Stokes方程[33]。這就是Boltzmann方程求解器能夠求解Navier-Stokes方程的基礎(chǔ)。GKS就是通過求解Boltzmann方程，甚至是求解簡化碰撞模型的Boltzmann模型方程(如Bhatnagar-Gross-Krook(BGK)方程)達(dá)到求解Navier-Stokes方程的目的。式(1)給出了BGK方程,其表達(dá)式為

(1)

式中：f為氣體分子分布函數(shù)；g為f對應(yīng)的平衡態(tài)分布函數(shù)；u為分子速度；x為空間位置；t為時(shí)間；τ為系統(tǒng)的弛豫時(shí)間，同分子間的平均碰撞時(shí)間有類似的物理意義。

f和g是分子速度u、空間位置x、內(nèi)部自由度ξ和時(shí)間t的函數(shù)。在氣體動(dòng)理學(xué)的描述下，所有的流動(dòng)宏觀物理量W都可以用氣體分布函數(shù)的矩來表示，表達(dá)式為

(2)

式中：ρ為氣體密度；U為氣體宏觀速度；E為單位質(zhì)量氣體的總能量。

微觀粒子分布函數(shù)能夠包含遠(yuǎn)多于宏觀量的信息，這就為數(shù)值格式提供了更好的擴(kuò)展平臺。

2.1.2 基于微分方程的直接建模方法

同傳統(tǒng)的數(shù)值格式不同，GKS并沒有直接求解BGK方程?？刂品匠滩]有顯式的出現(xiàn)在格式當(dāng)中，而是通過求解網(wǎng)格邊界的BGK方程的黎曼問題再經(jīng)由控制體內(nèi)的守恒律來體現(xiàn)控制方程的作用。這和Godunov格式比較接近。如果將以求解控制方程為目標(biāo)的數(shù)值方法稱為間接方法，以模擬物理問題為目標(biāo)將物理規(guī)律直接應(yīng)用于離散空間的數(shù)值方法就可以稱作直接方法。GKS可以看作一種直接方法，其使用的物理規(guī)律是BGK方程的局部精確解,即

(3)

式中：x0為初始時(shí)刻t=0時(shí)的位置；f0為初始時(shí)刻t=0的氣體分子分布函數(shù)；t′為積分變量。

可以看到這個(gè)積分解包含了空間和時(shí)間的所有信息。對式(3)取矩就能獲得x0位置沿n方向的完整的通量，即

(4)并且給定物理量的空間分布后，時(shí)間方向是精確的。

2.2GKS的方法特色

2.2.1 黏性和無黏部分統(tǒng)一處理

這個(gè)局部精確解是氣體系統(tǒng)的黎曼問題的解。而Godunov格式使用的是Euler方程的黎曼問題的解。Euler方程的黎曼問題只是氣體系統(tǒng)黎曼問題的一個(gè)特例，無法考慮黏性的影響。事實(shí)上，在采用傳統(tǒng)的迎風(fēng)格式得到Navier-Stokes方程解的過程中，黎曼問題的通量或者其他等價(jià)的方法得到的通量只是用來處理無黏項(xiàng)，而黏性和熱流項(xiàng)則需要通過中心格式來處理。理論上，對黏性項(xiàng)和無黏項(xiàng)使用不同的流動(dòng)分布函數(shù)是人工的，當(dāng)黏性和無黏相互作用非常強(qiáng)烈的時(shí)候可能會(huì)引入數(shù)值誤差。例如，迎風(fēng)格式的耗散特性可能會(huì)在大梯度區(qū)引起不穩(wěn)定機(jī)制使得無法得到穩(wěn)定的定常解。即使對黏性項(xiàng)和無黏項(xiàng)分別引入高階離散，這種分開處理的方式帶來的誤差也會(huì)由于不同的初始條件而難于消除，特別是當(dāng)輸運(yùn)過程和耗散加熱過程耦合的時(shí)候。這可能是預(yù)測傳熱的決定性因素。相反，對于GKS，無黏和黏性部分在局部精確解中都有體現(xiàn)，是由同一個(gè)分布函數(shù)(見式(3))來表示的。在連續(xù)流范圍內(nèi)，利用Chapmann-Enskog展開，精確解式(3)為

(5)

其中第2個(gè)式子表示碰撞守恒律。在給定了氣體分布函數(shù)的空間分布后，時(shí)間導(dǎo)數(shù)可以經(jīng)第2個(gè)式子得出。整個(gè)氣體系統(tǒng)的黎曼問題就是可解的。而且系數(shù)τ對應(yīng)的就是黏性作用。對于有間斷的初始條件，請參看文獻(xiàn)[32]。GKS能夠基于同一個(gè)初始條件得到包含黏性和無黏效應(yīng)的分布函數(shù)的演化，進(jìn)而通過氣體分布函數(shù)的積分得到熱流和黏性應(yīng)力。這種做法相比于Navier-Stokes方程通過物理量的梯度得到熱流和黏性應(yīng)力的方式對網(wǎng)格的敏感性更低。因此，氣體動(dòng)理學(xué)能夠?qū)?fù)雜黏性流動(dòng)問題給出更加貼切的描述。Ohwada研究了動(dòng)理學(xué)方程的數(shù)值通量，并闡述了碰撞過程的影響。根據(jù)Ohwada[34-35]給出的理論，這種GKS能夠被證明是對Navier-Stokes方程2階精度的格式。

2.2.2 多尺度的中心格式和迎風(fēng)格式混合

在傳統(tǒng)的動(dòng)理學(xué)通量分裂(KineticFluxVectorSplitting，KFVS)格式中，無碰撞的Boltzmann方程常常用來推導(dǎo)數(shù)值通量。而GKS在求解網(wǎng)格邊界的通量時(shí)，不僅考慮氣體分子的輸運(yùn)過程，還考慮了分子間的碰撞過程。具體來看邊界處的分布函數(shù)(見式(3))。前半部分代表著氣體分子的自由輸運(yùn)過程，一般用迎風(fēng)差分格式重構(gòu)。而后半部分代表著碰撞產(chǎn)生的平衡態(tài)的貢獻(xiàn)，一般用中心差分重構(gòu)。兩者通過系數(shù)e-t/τ及其變換公式加權(quán)得到最終的分布函數(shù)。這種構(gòu)造方式結(jié)合了中心格式和迎風(fēng)格式的優(yōu)點(diǎn)。2001年關(guān)于GKS的綜述文章給出了該格式計(jì)算激波和邊界層的算例[32]。這些數(shù)值結(jié)果顯示，GKS能夠很好地捕捉激波，同時(shí)也能高效地分辨邊界層。同時(shí)具有這2種能力得益于中心和迎風(fēng)特性的結(jié)合。

中心和迎風(fēng)是就離散方式而言的，它們背后的物理過程也值得關(guān)注。迎風(fēng)格式離散的部分f0，對應(yīng)著分子輸運(yùn)過程。這個(gè)過程有著明確的特征速度和特征方向。它描述了每個(gè)分子的獨(dú)立運(yùn)動(dòng)過程，是微觀尺度的描述。而中心格式離散的部分平衡態(tài)g的積分，描述了氣體分子的宏觀行為。這種宏觀尺度的行為是由于氣體分子大量的相互碰撞形成的。它可以由Navier-Stokes方程和Euler方程描述。這種宏觀和微觀尺度過程通過積分解能夠在離散空間中正確地反映出來。對于Navier-Stokes方程，真實(shí)的激波的尺度遠(yuǎn)遠(yuǎn)小于網(wǎng)格的尺度。這時(shí)，通過調(diào)節(jié)弛豫時(shí)間τ使得分子的輸運(yùn)過程f0起主要作用[36]，從而使激波厚度和網(wǎng)格尺度相當(dāng)。而當(dāng)處理光滑區(qū)域時(shí)，平衡態(tài)g的積分起主要作用?；贓uler方程的方法不具備這種多尺度特性，在計(jì)算高超聲速鈍體繞流問題時(shí)，往往伴隨著激波不穩(wěn)定性，出現(xiàn)所謂的“紅斑”現(xiàn)象。從已有的研究來看，這是基于Euler方程的數(shù)值方法所不可避免的。而GKS可以通過其粒子自由輸運(yùn)的部分有效地克服激波不穩(wěn)定性[37]。對高超聲速流，流體的粒子性的貢獻(xiàn)越來越突出?；诹Ｗ蛹僭O(shè)的高超聲速牛頓流理論的成功也從一個(gè)側(cè)面印證了粒子描述的重要性。多尺度特性在計(jì)算高超聲速過渡流問題時(shí)，顯得尤其重要。這是因?yàn)椋^渡流區(qū)域中這種宏觀和微觀的效應(yīng)都對流場有著顯著影響。單獨(dú)用宏觀方程或者只用粒子方法來處理這部分問題都有很大的困難。而通過積分解式(3)能將兩者在不同的網(wǎng)格尺度和時(shí)間尺度的離散空間結(jié)合起來，從而正確地描述這類問題。由于這部分涉及到許多稀薄非平衡問題，在此不作詳述。感興趣的讀者可以參看文獻(xiàn)[38]。

2.2.3 多維性質(zhì)

Euler方程的廣義黎曼問題非常困難，這限制了傳統(tǒng)格式的發(fā)展[13]。多數(shù)的黎曼求解器只考慮界面兩側(cè)物理量的點(diǎn)值，而不考慮斜率。而BGK方程的廣義黎曼問題的解，依然由式(3)表示。這就使得GKS能夠處理網(wǎng)格界面處更復(fù)雜的物理量分布問題。GKS能夠非常方便地應(yīng)用多維的粒子輸運(yùn)機(jī)制來計(jì)算通量。這種方式可以同時(shí)考慮界面的法向和切向梯度。多維的信息能夠改善激波穩(wěn)定性，對非結(jié)構(gòu)網(wǎng)格上的數(shù)值模擬有著重要意義[39]。同時(shí)，由于求解通量使用的是統(tǒng)一的分布函數(shù)，這樣得到的數(shù)值通量中既包含了無黏通量也包含黏性通量。而黏性通量不僅和網(wǎng)格界面的法向?qū)?shù)有關(guān)也和切向?qū)?shù)有關(guān)。因此，如果想要得到精確的黏性作用和傳熱，也需要使用多方向信息來構(gòu)造數(shù)值格式。事實(shí)上，如何精確模擬高超聲速黏性流動(dòng)仍然是當(dāng)下計(jì)算流體力學(xué)的一個(gè)難點(diǎn)，特別是如何精確模擬熱流。

多維GKS已經(jīng)被應(yīng)用于高超聲速層流黏性流問題。例如，馬赫數(shù)為8.03的雙錐高超聲速繞流問題[40]。這個(gè)軸對稱高超聲速雙錐分離流動(dòng)是驗(yàn)證高超聲速算法的一個(gè)經(jīng)典算例[41]。許多不同的高超聲速算法都模擬了這個(gè)問題[42-44]。雙錐外形的第1級錐的半頂角為25°，而第二級錐的半頂角為55°。試驗(yàn)中來流條件為：來流密度為ρ∞=0.654 5×10-3kg/m3，來流速度為U∞=2 664.0 m/s，來流溫度為T∞=185.56 K，壁面溫度為Tw=293.33 K。對應(yīng)的馬赫數(shù)和雷諾數(shù)分別為Ma=9.59，Re=13 090。

這個(gè)問題包含了各種復(fù)雜的流動(dòng)結(jié)構(gòu)，包括激波與激波相互作用、激波和邊界層相互作用以及流動(dòng)分離等。所有的這些問題都對數(shù)值方法的耗散特性非常敏感。多維GKS得到的數(shù)值結(jié)果和實(shí)驗(yàn)非常接近。

在這個(gè)流動(dòng)條件下，第1級錐體產(chǎn)生了附體斜激波，而第2級錐體產(chǎn)生了一個(gè)脫體弓形激波。這2個(gè)激波相互作用，形成一個(gè)透射激波打在靠近拐角的第2級錐的表面。錐體連接處和透射激波產(chǎn)生的逆壓梯度造成了很大的分離區(qū)，并由此產(chǎn)生了分離激波。激波相互作用在透射激波撞擊第2級錐的位置產(chǎn)生非常高的表面壓力和熱流率。圖1為雙錐繞流壓力分布圖[40],其中L1為第1級錐高度。

圖1 雙錐繞流壓力分布圖[40]Fig.1 Pressure distribution around double cone geometry [40]

圖2為沿圓錐表面的壓力pw和熱流qw的實(shí)驗(yàn)和計(jì)算結(jié)果?？梢钥闯?，數(shù)值結(jié)果的分離區(qū)的大小和圓錐表面的熱流都和實(shí)驗(yàn)數(shù)據(jù)吻合得很好，特別是捕捉到了沿著第2級錐表面的復(fù)雜流動(dòng)。

圖2 沿錐表面的壓力和熱流分布[40]Fig.2 Distribution of pressure and heat flux along cone surface [40]

高超聲速層流雙錐繞流計(jì)算結(jié)果捕捉到了激波相互作用、激波邊界層相互作用以及無黏和黏性相互作用等復(fù)雜的流動(dòng)結(jié)構(gòu)。計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)吻合得非常好。特別是物體表面準(zhǔn)確的熱流預(yù)測說明多維GKS是一種可靠的黏性流求解器。

2.2.4 滑移邊界條件

飛行器從低層大氣飛向高層大氣的過程中，大氣密度逐漸降低。在相同的飛行速度下，更小的大氣密度意味著更小的雷諾數(shù)。而雷諾數(shù)與努森數(shù)成反比。努森數(shù)是分子平均自由程和流場特征長度的比值，表征氣體的稀薄程度。一般認(rèn)為雷諾數(shù)較大、努森數(shù)較小時(shí)(小于0.001)能夠使用Navier-Stokes方程描述氣體流動(dòng)。但是由于高超聲速流動(dòng)中，往往存在大梯度區(qū)域，許多全局努森數(shù)很小的流動(dòng)問題在某些局部流場已經(jīng)出現(xiàn)了顯著的稀薄效應(yīng)。

許多研究都使用Navier-Stokes求解器加上滑移邊界條件來模擬這類近連續(xù)流問題。比較著名的滑移邊條件有FK(Fukui-Kaneko)模型[45]。這種邊條件一般先通過氣體動(dòng)理學(xué)理論推導(dǎo)得到滑移速度和壁面物理量梯度之間的關(guān)系，然后再應(yīng)用到宏觀方程上。當(dāng)然這種物理過程也能夠通過數(shù)值格式直接反映出來。在GKS中，所有計(jì)算都是基于粒子分布函數(shù)的，因此滑移邊界條件能夠自然地通過氣體和壁面作用模型加入到算法中。例如，壁面附近的入射分子可以通過流場內(nèi)部的插值得到，而出射分子可以通過壁面模型得到。再加上無穿透條件，就能讓數(shù)值方法自動(dòng)“推導(dǎo)”出滑移邊條件。這樣的做法相對于Navier-Stokes方程增加滑移邊條件的做法有幾個(gè)好處。

1) 使用范圍相對較廣

推導(dǎo)得到滑移邊界條件的宏觀表達(dá)式常常有某些假設(shè)。這些假設(shè)在不同的流場條件下可能會(huì)失效。而數(shù)值方法直接給出的動(dòng)理學(xué)邊界條件(自動(dòng)恢復(fù)滑移邊界條件)假設(shè)相對較少。著名的DSMC(Direct Simulation Monte Carlo)方法就是直接使用的動(dòng)理學(xué)邊條件。動(dòng)理學(xué)邊條件能夠適用于從自由分子流到過渡流各種流場條件。

2) 擴(kuò)展性好

由于免去了復(fù)雜的理論推導(dǎo)和假設(shè)，動(dòng)理學(xué)邊條件可以直接從分子和墻壁作用的微觀過程給出邊條件。

最典型的壁面和氣體分子作用模型有如下2種。

(1)漫反射模型

從墻面反射的氣體分子服從Maxwell分布，即：

(6)

式中：m為分子內(nèi)部自由度的維度；Tw為反射分布函數(shù)的溫度，其為壁面溫度；R為出射或反射的解釋；ρw為密度，由無穿透條件決定:

(7)

式中：n為固體邊界的法向，指向固體內(nèi)部。

(2)鏡面反射模型

反射的氣體分子的切向速度不變，法向速度反號，即：

(8)

式中：fin為入射分布函數(shù)。

考慮到氣體分子漫反射過程需要一定的作用時(shí)間，其反射的分布函數(shù)中的溫度將不再是固體溫度Tw而是某個(gè)弛豫溫度Tr。漫反射的氣體分子分布函數(shù)為

(9)

定義熱適應(yīng)系數(shù)α，取值范圍為0～1時(shí)的表達(dá)式為

α=(Ein-Er)/(Ein-Ew)

(10)

式中：Ein為入射氣體的總能量；Er為以溫度Tr反射的氣體內(nèi)能；Ew為墻壁對應(yīng)的內(nèi)能。它們的表達(dá)式分別為

(11)

式中：Uin為入射氣體速度；Tin為入射氣體溫度；Tr為反射氣體溫度；γ為來流氣體比熱比。

將這幾種反射模型結(jié)合，則固壁附近的氣體分子分布函數(shù)為

(12)

式中：σ為氣體分子漫反射的比例。

σ反映了反射氣體分子的切向動(dòng)量與固體切向動(dòng)量之間的適應(yīng)程度。由于GKS可以給出入射和反射的氣體分布函數(shù)，同時(shí)求解出其對應(yīng)的入射和出射的動(dòng)量和能量，它們的差值有效地給出了通過壁面的熱流和應(yīng)力(包括黏性應(yīng)力)。這種計(jì)算方式比Navier-Stokes方程通過導(dǎo)數(shù)求熱流和黏性力的方法更加符合物理過程。

GKS在計(jì)算空心圓柱-裙組合體高超聲速繞流問題的時(shí)候使用的就是上述動(dòng)理學(xué)邊條件[46]。圖3為空心圓柱-裙組合體的物體外形和網(wǎng)格。

圖3 空心圓柱-裙組合外形及計(jì)算網(wǎng)格示意圖[46]Fig.3 Schematic view of hollow-cylinder flare and computational mesh[46]

自由來流從左向右，定義為x方向，垂直于x方向的縱向方向?yàn)閥方向。來流參數(shù)根據(jù)文獻(xiàn)[47]給出：來流壓力為p∞=6.3Pa，來流馬赫數(shù)為Ma∞=9.91，來流溫度為T∞=51K，壁面溫度為Tw=293K。來流氣體比熱比γ=1.4, 普朗特?cái)?shù)Pr=0.718。

這組條件下整個(gè)計(jì)算區(qū)域都是層流。工作氣體使用N2，其黏性滿足：

(13)

式中：μref=1.663×10-5Pa·s；Tref=273.15 K；Ts=106.7K；

這項(xiàng)研究選取了無滑移邊界(NoSlip)條件和3種不同參數(shù)組合的動(dòng)理學(xué)邊條件：

1)σ=1,α=1，用KB(KineticBoundary)表示。

2)σ=0.5,α=1，用KB1表示。

3)σ=1,α=0.5，用KB2表示。

Navier-Stokes求解器模擬這個(gè)問題一般會(huì)過高地預(yù)測分離區(qū)的長度[47-49]。出現(xiàn)這種情況的原因是前緣的影響不能通過Navier-Stokes求解器精確地模擬，特別是沒考慮稀薄效應(yīng)的時(shí)候。一般而言，引入滑移邊條件會(huì)使得結(jié)果更加接近于實(shí)驗(yàn)和DSMC結(jié)果。圖4給出了x/L=0.3處的密度沿y方向的分布曲線，圖中ρinf=ρ∞。

圖4 x/L=0.3處的密度剖面 [46]Fig.4 Density profile at x/L = 0.3 [46]

由圖可以看到，滑移邊條件得到的曲線整體上比無滑移邊界條件的結(jié)果更接近實(shí)驗(yàn)和DSMC數(shù)據(jù)，而且圖中顯示前緣激波更靠近固體。也就是說，運(yùn)用滑移邊條件不僅能影響壁面附近的流動(dòng)，也會(huì)影響更遠(yuǎn)處的前緣激波。因此，滑移邊界條件高超聲速計(jì)算有顯著的影響。

表1為不同邊條件情況下分離點(diǎn)xsep和再附點(diǎn)xreat的位置。

由表1可知，當(dāng)條件為KB1時(shí)分離點(diǎn)向下游移動(dòng)，同時(shí)再附點(diǎn)向上游移動(dòng)。因此，分離區(qū)的長度明顯減小?？紤]不完全熱適應(yīng)邊條件(KB2)，計(jì)算得到的分離點(diǎn)位置也會(huì)向下游移動(dòng)，并且再附點(diǎn)也會(huì)向下游移動(dòng)。

表1 分離和再附點(diǎn)位置[46]

圖5為沿表面的壓力系數(shù)Cp和熱流(由Stanton數(shù)St表示)的分布曲線。

圖5 壓力系數(shù)和Stanton數(shù)的分布曲線[46]

GKS的計(jì)算結(jié)果給出的變化的趨勢和Markelov等[47]使用相同適應(yīng)系數(shù)的DSMC結(jié)果一致。

這項(xiàng)研究表明GKS能夠通過動(dòng)理學(xué)邊條件方便地實(shí)現(xiàn)各種滑移邊條件，顯著改善高超聲速流的計(jì)算效果。

2.3 未來發(fā)展

GKS取得了很大的進(jìn)步，但還有許多急需發(fā)展的方向。結(jié)合之前介紹的方法特色，未來的發(fā)展方向可以分為兩大類。

1)高階格式的發(fā)展。由于高階格式的實(shí)際表現(xiàn)性能由它的短板決定，這就要求在構(gòu)造和實(shí)現(xiàn)高階格式時(shí)必須全面考慮各種物理和數(shù)值問題，做到面面俱到。GKS的多維、無黏和黏性統(tǒng)一處理以及時(shí)間方向精確求解等特性，為高階格式提供了一個(gè)全面完整的物理描述。這種構(gòu)造格式的思想和高階離散方法已經(jīng)取得了一些初步進(jìn)展，出現(xiàn)了幾種高階格式版本[36,50-52]，也出現(xiàn)了一些湍流模擬應(yīng)用[53]。

但總體而言，GKS的湍流應(yīng)用相對較少，還需要繼續(xù)發(fā)展高分辨率格式并進(jìn)行湍流模型以及DNS研究。特別對于工程領(lǐng)域，有效實(shí)用的三階格式依然是一大難題,還需要更多的研究和探索。

2)多尺度問題模擬。多尺度的物理過程在高超聲速流動(dòng)中占有重要地位。例如湍流問題、稀薄氣體問題、輻射問題以及化學(xué)反應(yīng)問題等都是多尺度問題。過去的許多研究稀薄輻射化學(xué)反應(yīng)的工作都是建立在氣體動(dòng)理學(xué)理論基礎(chǔ)之上的。因此，利用GKS開發(fā)針對多尺度問題的數(shù)值格式也是一種行之有效的方法，GKS的多尺度特性也能發(fā)揮所長，這方面研究也取得了一些進(jìn)展。

針對湍流模擬，已有學(xué)者基于GKS提出了雷諾平均方程—(RANS)模型[54]。該項(xiàng)工作通過適當(dāng)給出系統(tǒng)的弛豫時(shí)間引入湍流黏性。對于稀薄氣體問題，統(tǒng)一氣體動(dòng)理學(xué)格式(Unified Gas-Kinetic Scheme，UGKS)[38]通過具體求解氣體分子分布函數(shù)能夠處理任意努森數(shù)范圍的稀薄氣體問題，但在工程應(yīng)用中，計(jì)算的不同區(qū)域非平衡程度有所不同。例如，在空心圓柱-裙組合體繞流問題中，前緣部分就需要考慮稀薄效應(yīng)，而其他部分則可以用宏觀的Navier-Stokes方法處理。這種基于微觀描述和宏觀方法的多尺度混合模擬方法也成為了一個(gè)熱點(diǎn)。這樣做的好處是既能正確地考慮非平衡效應(yīng)，又能提高計(jì)算效率。針對輻射問題，在統(tǒng)一氣體動(dòng)理學(xué)框架下也實(shí)現(xiàn)了光學(xué)薄到光學(xué)厚不同機(jī)制下的多尺度模擬[55]。GKS為多尺度數(shù)值方法提供了一個(gè)平臺。這也是高超聲速GKS的一個(gè)重要發(fā)展方向。

3 結(jié) 論

1)為了反映真實(shí)的物理過程，計(jì)算使用的控制方程越來越精確。

2)構(gòu)造計(jì)算方法需要考慮控制方程背后的物理過程。

3)基于Boltzmann方程的GKS能更基礎(chǔ)更精確的描述氣體系統(tǒng)，其通量求解器能處理廣義黎曼問題和黏性問題等。

4)GKS的多維多尺度特性為高階格式以及非平衡模擬提供了更好的擴(kuò)展平臺。

[1] Tsien H S. Similarity laws of hypersonic flows[J]. Journal of Mathematical Physics, 1946, 25(3): 247-251.

[2] Bertin J J, Cummings R M. Critical hypersonic aerothermodynamic phenomena[J]. Annual Review of Fluid Mechanics, 2006, 38: 129-157.

[3] Anderson Jr J D. Hypersonic and high-temperature gas dynamics[M]. 2nd. Virginia: AIAA, 2006: 200.

[4] Rasmussen M. Hypersonic flow[M]. New York: John Wiley & Sons, Inc., 1994: 322.

[5] Bushnell D M. Hypersonic flight experimentation-status and shortfalls[C]∥Future Aerospace Technology in the Service of the Alliance. Hampton,VA: NASA Langley Research Center, 1997.

[6] Moretti G, Abbett M. A time-dependent computational method for blunt-body flows[J]. AIAA Journal, 1966, 4(12): 2136-2141.

[7] Shang J S. Three decades of accomplishments in computational fluid dynamics[J]. Progress in Aerospace Sciences, 2004, 40(3): 173-197.

[8] MacCormack R. The effect of viscosity in hypervelocity impact cratering[J]. Journal of Spacecraft and Rockets, 2003, 40(5): 757-763.

[9] Shang J S, Hankey Jr W L. Numerical solution of the compressible Navier-Stokes equations for a three-dimensional corner[C]∥AIAA Aerospace Sciences Meeting. California: AIAA, 1977.

[10] Godunov S K. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[J]. Matematicheskii Sbornik, 1959, 89(3): 271-306.

[11] Steger J L, Warming R F. Flux vector splitting of the inviscid gas dynamic equations with application to finite-difference methods[J]. Journal of Computational Physics, 1981, 40(2): 263-293.

[12] van Leer B. Flux-vector splitting for the Euler equations[C]∥Eighth International Conference on Numerical Methods in Fluid Dynamics. Berlin: Springer Heidelberg, 1982: 507-512.

[13] Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[14] Roe P L. Characteristic-based schemes for the Euler equations[J]. Annual Review of Fluid Mechanics, 1986, 18(1): 337-365.

[15] Einfeldt B. On Godunov-type methods for gas dynamics[J]. SIAM Journal on Numerical Analysis, 1988, 25(2): 294-318.

[16] Toro E F, Spruce M, Speares W. Restoration of the contact surface in the HLL-Riemann solver[J]. Shock Waves, 1994, 4(1): 25-34.

[17] Harten A. High resolution schemes for hyperbolic conservation laws[J]. Journal of Computational Physics, 1983, 49(3): 357-393.

[18] Zhang H X. Non-oscillatory and non-free-parameter dissipation difference scheme[J]. Acta Aerodynamica Sinica, 1988, 6(2): 143-165 (in Chinese). 張涵信. 無波動(dòng)、無自由參數(shù)的耗散差分格式[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1988, 6(2): 143-165.

[19] Harten A, Engquist B, Osher S, et al. Uniformly high order accurate essentially non-oscillatory schemes, III[J]. Journal of Computational Physics, 1987, 131(2): 231-303.

[20] Liu X D, Osher S, Chan T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115(1): 200-212.

[21] Lele S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of Computational Physics, 1992, 103(1): 16-42.

[22] Zhong X L. Direct numerical simulation of hypersonic boundary-layer transition over blunt leading edges. I-A new numerical method and validation[C]∥AIAA, Aerospace Sciences Meeting & Exhibit. Reno, NV,USA: AIAA, 1997.

[23] Zhong X L. Direct numerical simulation of hypersonic boundary-layer transition over blunt leading edges. II-Receptivity to sound[C]∥AIAA, Aerospace Sciences Meeting & Exhibit. Reno, NV: AIAA, 1997.

[24] Zhong X. High-order finite-difference schemes for numerical simulation of hypersonic boundary-layer transition[J]. Journal of Computational Physics, 1998, 144(2): 662-709.

[25] Adams N A, Shariff K. A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems[J]. Journal of Computational Physics, 1996, 127(1): 27-51.

[26] Wang J, Wang L P, Xiao Z, et al. A hybrid numerical simulation of isotropic compressible turbulence[J]. Journal of Computational Physics, 2010, 229(13): 5257-5279.

[27] Davis R T. Numerical solution of the hypersonic viscous shock-layer equations[J]. AIAA Journal, 1970, 8(5): 843-851.

[28] Rubin S G, Lin A. Marching with the parabolized Navier-Stokes equations[J]. Israel Journal of Technology, 1980, 18: 21-31.

[29] Rubin S G, Tannehill J C. Parabolized/reduced Navier-Stokes computational techniques[J]. Annual Review of Fluid Mechanics, 1992, 24(1): 117-144.

[30] Prendergast K H, Xu K. Numerical hydrodynamics from gas-kinetic theory[J]. Journal of Computational Physics, 1993, 109(1): 53-66.

[31] Xu K, Prendergast K H. Numerical Navier-Stokes solutions from gas kinetic theory[J]. Journal of Computational Physics, 1994, 114(1): 9-17.

[32] Xu K. A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method[J]. Journal of Computational Physics, 2001, 171(1): 289-335.

[33] Harris S. An introduction to the theory of the Boltzmann equation[M]. New York: Dover Publications, Inc., 1994: 322.

[34] Ohwada T. On the construction of kinetic schemes[J]. Journal of Computational Physics, 2002, 177(1): 156-175.

[35] Ohwada T, Kobayashi S. Management of discontinuous reconstruction in kinetic schemes[J]. Journal of Computational Physics, 2004, 197(1): 116-138.

[36] Luo J, Xu K. A high-order multidimensional gas-kinetic scheme for hydrodynamic equations[J]. Science China Technological Sciences, 2013, 56(10): 2370-2384.

[37] Ohwada T, Adachi R, Xu K, et al. On the remedy against shock anomalies in kinetic schemes[J]. Journal of Computational Physics, 2013, 255: 106-129.

[38] Xu K, Huang J C. A unified gas-kinetic scheme for continuum and rarefied flows[J]. Journal of Computational Physics, 2010, 229(20): 7747-7764.

[39] Gnoffo P A. Updates to multi-dimensional flux recon-struction for hypersonic simulations on tetrahedral grids[C]∥48th AIAA Aerospace Sciences Meeting. Orlando, FL: AIAA, 2010.

[40] Xu K, Mao M, Tang L. A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow[J]. Journal of Computational Physics, 2005, 203(2): 405-421.

[41] Holden M S, Wadhams T P. A review of experimental studies for DSMC and Navier-Stokes code validation in laminar regions of shock/shock and shock boundary layer interaction including real gas effects in hyper-velocity flows[C]∥36th AIAA Thermophysics Conference. Orlando, FL: AIAA, 2003.

[42] Candler G V, Nompelis I, Druguet M C. Navier-Stokes predictions of hypersonic double-cone and cylinder-flare flow fields[C]∥AIAA Conference. Xi’an ,China: AIAA, 2001.

[43] Gaitonde D V, Canupp P W, Holden M S. Heat transfer predictions in a laminar hypersonic viscous/inviscid interaction[J]. Journal of Thermophysics and Heat Transfer, 2002, 16(4): 481-489.

[44] Knight D, Longo J, Drikakis D, et al. Assessment of CFD capability for prediction of hypersonic shock interactions[J]. Progress in Aerospace Sciences, 2012, 48: 8-26.

[45] Fukui S, Kaneko R. Analysis of ultra-thin gas film lubrication based on linearized Boltzmann equation: first report—derivation of a generalized lubrication equation including thermal creep flow[J]. Journal of Tribology, 1988, 110(2): 253-261.

[46] Li Q, Fu S, Xu K. Application of gas-kinetic scheme with kinetic boundary conditions in hypersonic flow[J]. AIAA Journal, 2005, 43(10): 2170-2176.

[47] Markelov G N, Kudryavtsev A N, Ivanov M S. Continuum and kinetic simulation of laminar separated flow at hypersonic speeds[J]. Journal of Spacecraft and Rockets, 2000, 37(4): 499-506.

[48] Chanetz B, Benay R, Bousquet J M, et al. Experimental and numerical study of the laminar separation in hypersonic flow[J]. Aerospace Science and Technology, 1998, 2(3): 205-218.

[49] Graur I A, Ivanov M S, Markelov G N, et al. Comparison of kinetic and continuum approaches for simulation of shock wave/boundary layer interaction[J]. Shock Waves, 2003, 12(4): 343-350.

[50] Li Q, Xu K, Fu S. A high-order gas-kinetic Navier-Stokes flow solver[J]. Journal of Computational Physics, 2010, 229(19): 6715-6731.

[51] Xuan L J, Xu K. A new gas-kinetic scheme based on analytical solutions of the BGK equation[J]. Journal of Computational Physics, 2013, 234: 524-539.

[52] Liu N. High-order accurate gas-kinetic schemes[D]. Beijing: Peking University, 2013 (in Chinese). 劉娜. 高精度氣體動(dòng)理學(xué)格式[D]. 北京: 北京大學(xué), 2013.

[53] Righi M. A modified gas-kinetic scheme for turbulent flow[J]. Communications in Computational Physics 2014, 16(1): 239-263.

[54] Kumar G, Girimaji S S, Kerimo J. WENO-enhanced gas-kinetic scheme for direct simulations of compressible transition and turbulence[J]. Journal of Computational Physics, 2013, 234: 499-523.

[55] Sun W, Jiang S, Xu K. Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[J]. Preprint, 2015.

Tel： 00852-23587433

E-mail: makxu@ust.hk

陳松澤 男， 博士研究生。主要研究方向：計(jì)算流體力學(xué)、稀薄氣體、氣體動(dòng)理學(xué)格式等。

Tel： 00852-59430791

E-mail: jacksongze@pku.edu.cn

*Corresponding author. Tel.： 00852-23587433 E-mail: makxu@ust.hk

Gas kinetic scheme in hypersonic flow simulation

XU Kun*, CHEN Songze

DepartmentofMathematics,TheHongKongUniversityofScienceandTechnology,HongKong,China

For hypersonic flow simulation, a review of computation fluid dynamics (CFD) and a summary of gas kinetic scheme are presented in this paper. The mechanism underlying the construction of gas kinetic scheme is clarified by comparing it with the traditional CFD method. The importance of direct modeling and the implementation of the physical laws in a discretized space are emphasized. Through some classical hypersonic applications in recent years, such as the shock/shock interaction, shock wave/boundary layer interaction, and hypersonic boundary layer separation problems, the advantages of the methodology are also demonstrated. As a trend of CFD, the gas kinetic scheme includes more fundamental physical laws in its algorithm construction, and the multiple scale nature makes the kinetic scheme feasible for the hypersonic applications. The principle of direct modeling and the methodology of constructing numerical schemes from mesoscopic or microscopic flow dynamics would benefit the development of reliable flow solvers, especially for the high speed flow.

hypersonic; gas kinetic scheme; slip boundary condition; shock/shock interaction; shock wave/boundary layer interaction

2014-07-25； Revised： 2014-09-16； Accepted： 2014-09-22； Published online： 2014-09-23 14:32

s： Hong Kong Research Grant Council (621011, 620813); State Key Laboratory of High-temperature Gas Dynamics Open Fund (2013KF03); Supported by State Key Laboratory for Turbulence and Complex Systems of Peking University

2014-07-25; 退修日期： 2014-09-16; 錄用日期： 2014-09-22; 網(wǎng)絡(luò)出版時(shí)間： 2014-09-23 14：32

www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0232.html

香港研資局基金(621011, 620813)； 高溫氣體動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金(2013KF03)； 北京大學(xué)湍流與復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室資助

Xu K, Chen S Z. Gas kinetic scheme in hypersonic flow simulation[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(1):135-146.徐昆，陳松澤．高超聲速計(jì)算中的氣體動(dòng)理學(xué)格式[J]．航空學(xué)報(bào), 2015, 36(1): 135-146.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2014.0232

V411.3; O354.4

A

1000-6893(2015)01-0135-12

徐昆 男， 博士， 教授。主要研究方向：計(jì)算流體力學(xué)、數(shù)值方法、高階格式、氣體動(dòng)理學(xué)格式和統(tǒng)一氣體動(dòng)理學(xué)格式等。

*通訊作者.Tel.： 00852-23587433 E-mail: makxu@ust.hk

URL： www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0232.html