曾嘉楠，李 琪，吳 雷

(南方科技大學(xué) 力學(xué)與航空航天工程系，深圳 518055)

0 引 言

經(jīng)典流體力學(xué)與空氣動力學(xué)建立在連續(xù)介質(zhì)假設(shè)之上，即宏觀微元中需包含足夠多的微觀分子。與之相反，稀薄氣體動力學(xué)描述的則是在此假設(shè)不成立時氣體的流動與熱力學(xué)等性質(zhì)?？伺瓟?shù)Kn是稀薄氣體研究中一個重要的無量綱參數(shù)，用來表征氣體的稀薄程度，通常定義為氣體分子平均自由程與流動特征長度的比值。因而，當(dāng)分子平均自由程較大（即氣體密度與壓強(qiáng)較?。r，或在流動特征長度極?。次⒓{米尺度下），氣體流動均可呈現(xiàn)出顯著的稀薄效應(yīng)。

由此，不難發(fā)現(xiàn)，稀薄氣體流動廣泛存在于航空航天、工業(yè)技術(shù)和能源環(huán)境等重要的國防、民生領(lǐng)域，如返回艙再入大氣層以及臨近空間飛行器高空域飛行時的繞流[1-5]、微機(jī)電系統(tǒng)中的氣體流動[6-9]以及頁巖氣開采中甲烷在復(fù)雜多孔介質(zhì)中的流動[10]等。相關(guān)科學(xué)技術(shù)的發(fā)展對我國國防、國民經(jīng)濟(jì)和社會發(fā)展意義重大，如下所述：

1）臨近空間飛行器既可以避免絕大部分的地面攻擊，同時也能夠有效實(shí)施對地攻擊和對航天器的打擊，因此臨近空間是進(jìn)行空中軍事活動的理想?yún)^(qū)域，發(fā)展?jié)摿O大，我國和不少其他國家已經(jīng)將其納入信息化武器裝備體系[11]。在70 km以上的飛行空域，復(fù)雜構(gòu)型的高超聲速飛行器表面存在極為復(fù)雜的流動，因而需要建立新的流動模型和計(jì)算方法，來預(yù)測復(fù)雜的稀薄氣體效應(yīng)、高溫真實(shí)氣體效應(yīng)以及解決激波/膨脹波與黏性作用的耦合干擾等問題。深入研究上述實(shí)際應(yīng)用中的稀薄氣體流動規(guī)律是優(yōu)化飛行器設(shè)計(jì)的關(guān)鍵所在，具有重要意義。

2）機(jī)電系統(tǒng)微小集成化是提高能源利用效率的重要手段，因?yàn)槲⑿』纱蠓葴p少傳熱傳質(zhì)阻力從而顯著提高工質(zhì)的相關(guān)輸運(yùn)能力。近年來，微流動系統(tǒng)在電子、化學(xué)、生物等領(lǐng)域已經(jīng)取得了巨大的發(fā)展并得到了廣泛的應(yīng)用，但是微尺度下的空氣動力學(xué)理論和模擬方法仍不完善，還存在許多亟待解決的問題[12]：如在電容式膜片壓力計(jì)中，當(dāng)被測氣體的溫度與壓力計(jì)的溫度不同時，稀薄氣體自動從低溫區(qū)向高溫區(qū)蠕動，從而導(dǎo)致壓力測量偏差，需要考慮稀薄氣體效應(yīng)進(jìn)行修正[13]；工作在低壓風(fēng)洞的測量設(shè)備也需要考慮相應(yīng)的稀薄氣體效應(yīng)修正以提高測量精度。

3）隨著自然資源逐漸枯竭和環(huán)境破壞愈發(fā)嚴(yán)重，能源轉(zhuǎn)型和能源利用效率提升成為各國科技發(fā)展的重要目標(biāo)。相較于傳統(tǒng)化石能源，以甲烷為主要成分的頁巖氣，因其具有彌補(bǔ)常規(guī)能源儲量及產(chǎn)量、減少溫室氣體排放的巨大潛力，在世界能源供給中變得愈發(fā)重要。值得一提的是，我國頁巖氣儲備位列全球第一，但復(fù)雜的地質(zhì)條件要求對其輸運(yùn)機(jī)制有更清晰深入的理解，才能實(shí)現(xiàn)經(jīng)濟(jì)有效的開采；雖然頁巖氣埋藏在地表深處，氣體壓力大，但是頁巖孔隙直徑在納米尺度，因而同樣需要考慮稀薄氣體效應(yīng)[14]。

圖1列舉了典型稀薄氣體流動中納維-斯托克斯方程的數(shù)值預(yù)測結(jié)果與實(shí)驗(yàn)測量存在明顯差距的例子。由此可以看出，基于連續(xù)性假設(shè)得到的流體力學(xué)方程不再適用于此類問題，其根本原因在于，在稀薄效應(yīng)較強(qiáng)的區(qū)域內(nèi)，由宏觀量的一階導(dǎo)數(shù)構(gòu)成的線性本構(gòu)關(guān)系不再準(zhǔn)確。另一方面，直接模擬所有氣體分子運(yùn)動與相互作用的分子動力學(xué)方法，提供了一個更為準(zhǔn)確地模擬稀薄氣體流動的途徑。然而，因其龐大的計(jì)算量，該方法僅能用于極小時空尺度下的氣體輸運(yùn)問題。因此，對于更為廣泛存在的介觀尺度下稀薄氣體流動，如何建立準(zhǔn)確的動理論模型方程并有效求解，則是目前最為關(guān)切與亟待解決的兩個問題。

圖1 典型稀薄氣體流動及納維-斯托克斯方程的預(yù)測誤差：第一列圖顯示，當(dāng)馬赫數(shù)大于2時，納維-斯托克斯方程低估正激波的厚度[15]；第二列圖顯示，當(dāng)聲波振動頻率接近氣體分子碰撞頻率時，納維-斯托克斯方程低估聲波接收器的聲強(qiáng)[16-17]；第三列圖顯示，低壓下納維-斯托克斯方程可能遠(yuǎn)遠(yuǎn)低估微納米通道的質(zhì)量流量[18]Fig. 1 Typical rarefied gas flows and the incapability of Navier-Stokes equations. Figures in the first column show that, when the Mach number is larger than 2, the Navier-Stokes equations underpredict the thickness of normal shock wave[15]. Figures in the second column show that the Navier-Stokes equations underpredict the sound pressure at the receiver[16-17]. Figures in the last column show that the Navier-Stokes equations might underestimate the mass flow rate significantly[18]

（1）研究現(xiàn)狀。

一般來說，傳統(tǒng)的納維-斯托克斯方程只能描述連續(xù)流區(qū)（Kn<0.001）的氣體流動。隨著克努森數(shù)的增加，氣體流動將分別經(jīng)歷滑移流區(qū)（0 .00110）。對于處于滑移區(qū)的流動，氣體的稀薄效應(yīng)集中于邊界處，因此采用適當(dāng)?shù)乃俣然坪蜏囟忍S邊界條件來求解傳統(tǒng)宏觀方程即可得到較為準(zhǔn)確的結(jié)果。過渡流區(qū)內(nèi)的稀薄氣體效應(yīng)由邊界擴(kuò)展到流動內(nèi)部；自由分子流區(qū)中的氣體相互碰撞極少，直接模擬蒙特卡羅方法（direct simulation Monte Carlo，DSMC）是一種常用的求解方法。然而，在許多實(shí)際問題中，氣體的稀薄程度可能存在跨尺度的現(xiàn)象。例如，中國空氣動力研究與發(fā)展中心的江定武等發(fā)現(xiàn)[19]，返回艙載入過程中會經(jīng)歷從自由分子流到近地面連續(xù)流的多流動尺度的變化，同時，當(dāng)?shù)乜伺瓟?shù)在飛行器迎風(fēng)面和背風(fēng)面有若干數(shù)量級的差別，因而同時存在混合多尺度的流動[20]。因此，一個能夠在全流域下統(tǒng)一描述稀薄氣體流動的模型方程與求解方法十分必要。歷史上，曾嘗試通過對速度分布函數(shù)及宏觀量的漸近展開以及不同階次取矩等方法，得到了十余個宏觀方程來描述稀薄氣體行為[17]。相較于納維-斯托克斯方程，這些方程的數(shù)量和復(fù)雜性大幅增加，而穩(wěn)定性大幅降低。例如，Burnett方程和矩方程都曾被提出用來描述稀薄氣體流動[21-23]，但由于方程本身的不穩(wěn)定性和邊界條件的難確定性，未能有效應(yīng)用于全流域稀薄流模擬。相反，基于氣體動理論方法的建模和數(shù)值方法則在全流域稀薄氣體模擬中取得了一定的成功。

作為最為重要的氣體動理論方程，玻爾茲曼方程描述了單原子氣體的速度分布函數(shù)在時間與空間中的演化規(guī)律。玻爾茲曼方程的求解，在航空航天、微機(jī)電系統(tǒng)和頁巖氣開發(fā)等領(lǐng)域中都發(fā)揮著重要作用。其中，基于統(tǒng)計(jì)隨機(jī)方法的DSMC計(jì)算，是模擬稀薄氣體流動的重要工具[24]。對于單原子氣體的模擬，DSMC方法已在數(shù)學(xué)上被證明在粒子數(shù)趨于無窮的條件下，收斂于玻爾茲曼方程的解[25-26]。玻爾茲曼方程中氣體分子的速度分布函數(shù)是定義在平均自由程尺度上，由于DSMC方法將氣體分子自由運(yùn)動與相互碰撞這兩個過程進(jìn)行了解耦，因此該方法的模擬中要求空間網(wǎng)格尺寸小于氣體平均自由程，且時間步長小于平均碰撞時間。而這直接導(dǎo)致了DSMC方法對近連續(xù)流模擬的低效性。此外，由于流場的宏觀量信息需要通過統(tǒng)計(jì)抽樣得到，因此DSMC模擬中需進(jìn)行大量系綜平均去除噪聲，特別地，當(dāng)宏觀物理量自身數(shù)值相對于其統(tǒng)計(jì)微觀量的漲落較小時，例如低速流動或較小溫差的系統(tǒng)，DSMC方法的時間成本將變得難以負(fù)擔(dān)。

相比于DSMC等統(tǒng)計(jì)模擬方法，確定性數(shù)值方法在近幾十年得到快速發(fā)展并展現(xiàn)出強(qiáng)大優(yōu)勢，我國學(xué)者在此方面貢獻(xiàn)良多。但由于玻爾茲曼方程碰撞項(xiàng)的復(fù)雜性，除了少數(shù)方法可以求解玻爾茲曼方程外[27-31]，大部分確定性方法求解的是簡化模型方程，如BGK模型[32]、ES-BGK模型[33-34]和Shakhov模型[35]等。中國空氣動力研究與發(fā)展中心的李志輝和張涵信[36]最早應(yīng)用氣體動理論統(tǒng)一數(shù)值算法求解模型方程，并建立了三維復(fù)雜飛行器高馬赫數(shù)繞流問題的統(tǒng)一算法大規(guī)模并行計(jì)算研究平臺，已應(yīng)用于返回艙回收和天宮一號離軌、再入、解體等航天工程[37-39]。香港科技大學(xué)徐昆首創(chuàng)的統(tǒng)一氣體動理論格式（unified gas kinetic scheme，UGKS）耦合處理了模型方程的遷移項(xiàng)和碰撞項(xiàng)，相對于傳統(tǒng)離散速度方法而言，在網(wǎng)格離散尺度遠(yuǎn)大于分子平均自由程時依然可以在連續(xù)流時恢復(fù)納維-斯托克斯方程的解，從而在計(jì)算多尺度問題時可以采用更粗?；臅r空離散[40]。華中科技大學(xué)的郭照立、西北工業(yè)大學(xué)的鐘誠文、新加坡國立大學(xué)舒昌等課題組在UGKS的基礎(chǔ)上做了進(jìn)一步發(fā)展[41-43]。最近，蘇微和吳雷等提出的合成迭代算法（general synthetic iterative scheme，GSIS）耦合了宏觀方程以實(shí)現(xiàn)全流場信息交換，具有良好的“漸近保持”和“快速收斂”性質(zhì)，不但可以使用大空間網(wǎng)格和時間步長，而且可以在短短幾十步迭代內(nèi)找到穩(wěn)態(tài)解，極大地提高了數(shù)值精度和效率[44-47]。更重要的是，該方法對碰撞項(xiàng)的具體形式?jīng)]有要求，可以直接針對玻爾茲曼方程進(jìn)行快速求解。

（2）存在的一些問題。

稀薄氣體研究主要涉及動理學(xué)建模和多尺度計(jì)算兩個方面，本文僅關(guān)注分子氣體動理論模型中存在的一些問題。眾所周知，目前大部分氣體動理論方程是針對單原子氣體提出的[48-50]，而實(shí)際氣體大多是分子氣體，即一個氣體分子包含兩個或兩個以上的原子，例如占據(jù)地球大氣絕大部分比例的N2和O2，占據(jù)火星大氣絕大部分比例的CO2等。相對于僅具有平動自由度的單原子氣體，分子氣體中還存在轉(zhuǎn)動與振動的內(nèi)能形式，并且各類型自由度的能量可以相互轉(zhuǎn)化。對于常見氣體，其轉(zhuǎn)動能在室溫下已充分激發(fā)，而振動能也在上千度高溫下逐漸激發(fā)。隨著溫度的繼續(xù)升高，分子氣體的離解與電離也將逐漸發(fā)生，與此同時，化學(xué)反應(yīng)也通常存在于此。此時，高溫真實(shí)氣體效應(yīng)顯著。而單原子玻爾茲曼方程無法包含這些過程，因此需要更復(fù)雜的模型方程來給出相應(yīng)的模擬。目前，大部分以此為目標(biāo)的建模工作是針對連續(xù)流提出的，例如對高超聲速飛行器周圍氣體的模擬。由于分子氣體的熱流由平動能量和內(nèi)部能量的輸運(yùn)與相互轉(zhuǎn)換共同作用導(dǎo)致，因此需要通過額外的能量方程來描述高溫氣體能量輸運(yùn)，并且保證能量交換速率和總熱導(dǎo)率與實(shí)驗(yàn)數(shù)值一樣，此類連續(xù)流下的模擬方法就可以較為準(zhǔn)確地模擬高溫真實(shí)氣體效應(yīng)，國內(nèi)外學(xué)者在這方面已經(jīng)做了大量工作方面已經(jīng)開展了大量工作[51]。

然而，對于臨近空間飛行器等類似問題而言，稀薄氣體效應(yīng)也非常顯著[19]，因此有必要將稀薄氣體效應(yīng)和高溫真實(shí)氣體效應(yīng)結(jié)合起來考慮。對此，王承書和烏倫貝克建立的WCU方程[52]給出了相對完備的描述，該方程以經(jīng)典力學(xué)方式描述平動能量，而以量子力學(xué)方式描述內(nèi)部（轉(zhuǎn)動與振動）能量，并對每個能級下的速度分布函數(shù)建立方程。顯然，該方法使得方程組自由度過于龐大，致使其求解變得相當(dāng)困難從而難以應(yīng)用于實(shí)際問題。因此建立相對簡化的模型方程尤為必要。與此同時，由于各類自由度的特征碰撞時間不同，導(dǎo)致稀薄氣體效應(yīng)的程度也不同，因此動理論建模難度也大為增加。

對于分子氣體的動理論建模， 如何在模型方程中準(zhǔn)確反應(yīng)氣體的各個輸運(yùn)系數(shù)是首要問題。相對于單原子氣體，由于更多內(nèi)部自由度的存在，分子氣體中出現(xiàn)了更多的輸運(yùn)過程及相應(yīng)參數(shù)。其一為體積黏性，其表征了氣體在膨脹壓縮的做功過程中能量與內(nèi)部自由度轉(zhuǎn)換的阻力，這一阻力來源于有限時間的能量弛豫過程；其二為轉(zhuǎn)動、振動熱導(dǎo)率，即內(nèi)部自由度同時參與了熱流的傳遞，與熱流的弛豫速率緊密相關(guān)。實(shí)際上，在分子氣體動理論的建模中，正確的體積黏性是容易實(shí)現(xiàn)的，只需保證正確的內(nèi)部能量（轉(zhuǎn)動、振動等）和平動能量的弛豫速率即可。但是，實(shí)現(xiàn)正確的平動、轉(zhuǎn)動和振動熱導(dǎo)率卻非常困難。本文以僅有轉(zhuǎn)動自由度的分子氣體為例說明問題所在。歷史上，Eucken首先將氣體總熱導(dǎo)率 κ與剪切黏性 μ的關(guān)系表示為：

其中：m是 分子質(zhì)量；kB是玻爾茲曼常數(shù)；dr是分子的轉(zhuǎn)動自由度的數(shù)量，如N2在常溫下具有兩個轉(zhuǎn)動自由度，則dr= 2； κtr和 κrot分別是氣體的平動和轉(zhuǎn)動熱導(dǎo)率；ftr、frot和feu分別是平動、轉(zhuǎn)動和總Eucken因子。對于單原子氣體，ftr非常接近2.5，但對于分子氣體，ftr和frot的數(shù)值很難確定。最初Eucken從簡單分子動理論出發(fā)，認(rèn)為：

后來，人們意識到氣體內(nèi)部熱導(dǎo)率與分子擴(kuò)散有關(guān)，將內(nèi)部Eucken因子改寫為：

其中 ρ為氣體密度，D為氣體自擴(kuò)散系數(shù)，Sc為斯密特?cái)?shù)。最后，Mason和Monchick根據(jù)王承書和烏倫貝克建立的WCU方程[52]，推導(dǎo)出Eucken因子的表達(dá)式[53]，但與實(shí)驗(yàn)數(shù)據(jù)比較仍然存在差距[54]。由此可知，目前并不存在簡明準(zhǔn)確的僅依賴于氣體參數(shù)與其他輸運(yùn)系數(shù)的各個Eucken因子的表達(dá)關(guān)系，因此各熱導(dǎo)率分量均需直接與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行匹配，才能保證相應(yīng)輸運(yùn)系數(shù)的準(zhǔn)確性。然而，雖然可以比較容易地測得總熱導(dǎo)率，但從實(shí)驗(yàn)獲取熱導(dǎo)率各分量本身便是相當(dāng)有難度的任務(wù)。

WCU方程的復(fù)雜性限制了其在實(shí)際問題中的應(yīng)用，而DSMC方法則相對容易地從單原子氣體模擬推廣至分子氣體的輸運(yùn)。即便如此，對于極為復(fù)雜的分子氣體碰撞過程，DSMC方法依然采用了適當(dāng)簡化的模型來處理，從而提高計(jì)算效率。目前通用的分子碰撞模型最早由Borgnakke和Larsen提出[55]。此后，人們提出各種各樣的改進(jìn)，主要目的均在于正確實(shí)現(xiàn)分子各能級間的能量交換速率，從而獲得正確的分子氣體的體積黏性[56-58]。然而，當(dāng)前的DSMC方法并沒有包含可以單獨(dú)調(diào)節(jié)熱導(dǎo)率的相應(yīng)機(jī)制與參數(shù)。與此同時，在實(shí)際中使用DSMC模擬分子氣體流動時，該碰撞模型對應(yīng)得到的平動和轉(zhuǎn)動熱導(dǎo)率是否正確尚未可知。因此，有必要對DSMC模擬分子氣體的相應(yīng)模型進(jìn)行檢查，這里的檢查有雙重意義：1）檢查DSMC模擬得到的總熱導(dǎo)率是否與實(shí)驗(yàn)測量數(shù)據(jù)一致，若不一致，則DSMC即便在近連續(xù)流下也會給出錯誤解。2）即使DSMC能給出正確的總熱導(dǎo)率，也還需要檢查ftr和frot的數(shù)值是否與實(shí)際氣體一致。

在近連續(xù)流條件下，能量在平動自由度與內(nèi)部自由度的分配趨于一致，因而平動溫度與轉(zhuǎn)動和振動溫度的差別極小，此時保證總體熱導(dǎo)率的正確性，即可得到相應(yīng)正確的結(jié)果。因而，在稀薄效應(yīng)較強(qiáng)的區(qū)域，保證各個熱導(dǎo)率分量的正確性就顯得尤為關(guān)鍵。例如，愛因斯坦等很早就發(fā)現(xiàn)，對于熱蠕動等稀薄氣體流動，溫度梯度引起的氣體流動只與ftr有關(guān)[59]。最近李琪和吳雷等發(fā)現(xiàn)，目前流行的、公認(rèn)的DSMC并不一定能做到以上兩點(diǎn)[54,60]，因此DSMC在分子氣體流動中的模擬結(jié)果并不可靠。同樣地，對于確定性方法，人們在原有單原子氣體模型方程的基礎(chǔ)上，提出了各種描述分子氣體稀薄效應(yīng)的動理論簡化模型方程，如ES-BGK模型[33-34]、Rykov模型[61]等。與DSMC方法相似的是，這類模型對于分子氣體的拓展同樣專注于多種能量相互轉(zhuǎn)換的速率（即氣體體積黏性的正確性），而忽略了熱流弛豫過程的速率（即熱導(dǎo)率分量的正確性）。具體而言，除了Rykov模型可以獨(dú)立調(diào)整平動和內(nèi)部Eucken因子以對應(yīng)真實(shí)參數(shù)外，其他模型只能調(diào)整總熱導(dǎo)率，因此不適合用于研究熱蠕動等稀薄氣體現(xiàn)象。

（3）本文目的和結(jié)構(gòu)。

分子氣體稀薄效應(yīng)的準(zhǔn)確刻畫是稀薄氣體動理學(xué)和高溫氣體動力學(xué)研究的前沿和熱點(diǎn)，在這個背景下，本文首先介紹分子氣體動理論建模的思路和最新研究進(jìn)展，并討論現(xiàn)有的分子氣體模型方程在稀薄氣體流動中的適用性與準(zhǔn)確性，為今后多組分分子氣體和化學(xué)反應(yīng)的非平衡流動的建模和模擬奠定堅(jiān)實(shí)的理論基礎(chǔ)。以下第1節(jié)介紹單原子氣體玻爾茲曼方程的基本性質(zhì)、輸運(yùn)系數(shù)的嚴(yán)格推導(dǎo)以及動理學(xué)建模的基本思路。特別地，我們指出，幾乎所有的簡化模型都可以從Gross-Jackson模型非線性化而來。第2節(jié)介紹王承書和烏倫貝克提出的WCU方程，推導(dǎo)分子氣體的體積黏性和平動、轉(zhuǎn)動熱導(dǎo)率系數(shù)，并指出DSMC碰撞模型的缺陷。第4節(jié)介紹WCU方程的線性化Hanson-Morse模型以及主流的非線性模型，并在經(jīng)典的稀薄氣體流動中測試模型的精度。第4節(jié)和第5節(jié)分別討論熱弛豫速率的不確定性對稀薄氣體流動的影響，以及從分子動力學(xué)模擬和實(shí)驗(yàn)中減小甚至消除熱弛豫速率不確定性的方法。最后做總結(jié)和展望。

1 單原子氣體玻爾茲曼方程及其簡化模型

在氣體動理論中，單原子氣體在相空間中的概率密度由速度分布函數(shù)f(t,x,v)表 示，是時間t、空間坐標(biāo)x=(x1,x2,x3) 和 分子速度v=(v1,v2,v3)的函數(shù)。定義f(t,x,v)dxdv是 體積為 dxdv的相空間上的分子數(shù)，則氣體的數(shù)密度n、 宏觀速度u、 溫度T、壓力張量pij和熱流q可以分別通過對速度分布函數(shù)求矩得到。

其中，c=v?u是氣體分子熱運(yùn)動速度，即氣體分子速度與當(dāng)?shù)睾暧^速度的矢量差，而c是熱運(yùn)動速率。注意本文其他矢量也如此表示：黑體表示矢量，斜體表示矢量的模。定義應(yīng)力偏量為：

其中，p=nkBT， δij為克羅內(nèi)克函數(shù)。

對于稀疏氣體（分子間距遠(yuǎn)遠(yuǎn)大于分子直徑），在外部施加的加速度a=(a1,a2,a3)的作用下，速度分布函數(shù)的演化由玻爾茲曼方程描述：

方程（6）左邊的三項(xiàng)分別表示速度分布函數(shù)在時間上的變化、在速度作用下在物理空間的變化以及在外力作用下在速度空間的變化，右邊表示使得分布函數(shù)趨于平衡態(tài)的氣體分子的碰撞過程。在玻爾茲曼方程中，兩體碰撞的形式為：

其中，v和v?分 別是碰撞前兩個分子的速度，而v′和v′?是碰撞后的速度。由于碰撞前后兩分子的距離足夠遠(yuǎn)，因此其相互作用可以忽略不計(jì)。根據(jù)動量和能量守恒定律，碰撞前后速度關(guān)系如下：碰撞示意圖見圖2。其中，碰撞前兩分子的相對速度為vr=v?v?， 碰撞后的相對速度為v′?v′?。 Ω為定義在單位球空間的矢量，其與碰撞后的相對速度同方向。于是相對速度的偏轉(zhuǎn)角 θ與碰撞前相對速度滿足如下關(guān)系：

圖2 (左)兩體碰撞前后的速度分布：由于動量和能量守恒，碰撞前后的相對速度分布在球體上并且通過球心；(右上)中心力場作用下的經(jīng)典兩體碰撞示意圖，其中b 為 瞄準(zhǔn)距離，k為沿兩分子間最短距離方向的單位矢量；(右下) 直徑為σ的硬球分子的兩體碰撞Fig. 2 (Left) Velocity redistribution after the binary collision.Due to the conservation of momentum and energy, the pre- and post-collision relative velocities fall in a sphere and pass through the sphere center. (Top right) Classical binary collision between molecules with central force, where b is the aiming distance, andk is the unit vector along the minimum distance between two colliding molecules. (Bottom right) Binary collision between HSmolecules of the diameter σ

最后，碰撞核B(θ,vr)是碰撞偏轉(zhuǎn)角度和相對速度的函數(shù)，具體形式取決于分子間的作用力。

1.1 偏轉(zhuǎn)角和微分散射截面

假設(shè)分子間通過中心力場作用，作用勢 φ(r)已知，其中r為分子間距，則偏轉(zhuǎn)角可以通過經(jīng)典力學(xué)和量子力學(xué)兩種方式求解。若氣體溫度不低（如氦氣的溫度高于100 K），兩種方式得到的輸運(yùn)系數(shù)（如黏性和熱導(dǎo)率）相同[62]。這里僅介紹經(jīng)典力學(xué)的計(jì)算結(jié)果，如圖2所示，偏轉(zhuǎn)角可以表示為：

其中，W=b/r為瞄準(zhǔn)距離b與 分子間距r的比值，積分上限W1對應(yīng)于最短分子間距，即式（10）中括號內(nèi)表達(dá)式等于零的方程的正根。

在氣體動理論中，經(jīng)?？紤]如下形式的逆冪律分子作用勢：

其中，K表征分子間相互作用的強(qiáng)度。從式（10）可知，偏轉(zhuǎn)角只與變量s有 關(guān)，即 θ=Θ(s) ：

微分散射截面定義為：

而碰撞核為：

當(dāng)式（11）中 η=5時，即為麥克斯韋分子，此時碰撞核與相對碰撞速度無關(guān)，記為對于硬球分子模型，可看作式（11）中 η=∞。從圖2可以看出，偏轉(zhuǎn)角可通過如下公式確定：b=σcos(θ/2) ，其中 σ為硬球直徑。因此微分散射截面為 σD=σ2/4， 而碰撞核為B(θ,vr)=σ2vr/4。對于一般的氣體，其行為介于麥克斯韋分子和硬球分子之間。

1.2 熵增原理

對于任意關(guān)于分子速度的函數(shù)Ψ (v)，對玻爾茲曼碰撞項(xiàng)（式（7））在速度空間積分，具有如下對稱性：

其 中 Δ[Ψ]=Ψ(v?)+Ψ(v)?Ψ(v′?)?Ψ(v′)。 若 Ψ滿 足則稱其為碰撞不變量。根據(jù)質(zhì)量、動量和能量守恒，可知Ψ =1、v、v2為碰撞不變量，而碰撞不變量的線性組合也是碰撞不變量。

定義H函數(shù)為：

H是與氣體系統(tǒng)的熵相關(guān)的標(biāo)量。在無外力的情況下，H的時間導(dǎo)數(shù)可寫為：

式中，n為 系統(tǒng)表面微元 dS的外法線方向。對于孤立系統(tǒng)，上式右端第一項(xiàng)為零；利用方程（15），右端項(xiàng)第二項(xiàng)可寫為：

因?yàn)榕鲎埠薆非負(fù)，且對于任意的兩個正整數(shù)a和b， 不等式 (a?b)(lna?lnb)≥0恒成立，所以有：

式（17）表明孤立系統(tǒng)的H函數(shù)不會減小，這就是著名的玻爾茲曼的熵增原理。H隨時間單調(diào)增加，但存在有限的上界，上界對應(yīng)式（17）中等號成立時的平衡態(tài)，即 lnf(v?)+lnf(v)=lnf(v′?)+lnf(v′)。 因此， lnf也是碰撞不變量，可以表示為五個基本碰撞不變量的線性組合： lnf=α1+α2·v+α3v2；給定系統(tǒng)的密度、速度和溫度，參數(shù)α1、α3和α2可以唯一確定。此時，麥克斯韋平衡態(tài)速度分布函數(shù)為：

1.3 線性化玻爾茲曼方程

線性化玻爾茲曼方程在氣體動理論中具有重要地位。第一，線性化玻爾茲曼碰撞項(xiàng)的本征值與本征函數(shù)，不僅在漸近展開推導(dǎo)納維-斯托克斯方程的過程中至關(guān)重要，而且是發(fā)展簡化動理學(xué)模型方程的源頭理論。第二，在許多微機(jī)電系統(tǒng)中，氣體的壓力梯度與溫度梯度非常小，因此使用線性化方程可以高效準(zhǔn)確地模擬微流動。第三，雖然玻爾茲曼方程是單粒子速度分布函數(shù)確定性演化的平均方程，但是在某些問題中（例如：第5.2節(jié)中的瑞利-布里淵散射）可以通過線化玻爾茲曼方程研究粒子漲落帶來的影響。

通常將速度分布函數(shù)在全局平衡態(tài)

下展開為：

其中擾動量 φ滿足約束條件 |φ|?1。不考慮外力項(xiàng)，只保留 φ的一次項(xiàng)，玻爾茲曼方程（6）可線性化為如下形式：

這里的速度和空間坐標(biāo)分別通過最概然速率vm和特征尺寸 ?進(jìn)行無量綱化，時間用 ?/vm歸一化；無量綱分子速度為：

其中：

為氣體分子的最概然速率，即與麥克斯韋速率分布的極大值對應(yīng)的速率。

對于麥克斯韋分子，線性化玻爾茲曼碰撞項(xiàng)的本征值與本征函數(shù)為[52]：

式中，Pl(x) 是勒讓德多項(xiàng)式，是索南多項(xiàng)式，為 ξ方向的球諧函數(shù)。

前三項(xiàng)本征值為 λ00=λ01=λ10=0，分別代表三大守恒律。另外兩個非常重要的本征值為：

這兩個本征值決定了剪切黏性系數(shù)和熱導(dǎo)率的大小。若歸一化的速度 ξ 在z軸 的投影為 ξz=ξcosθ（其中 θ為極距角），與上述本征值對應(yīng)的本征函數(shù)為：

其與氣體的分子數(shù)密度、速度、溫度、壓力和熱流的擾動密切相關(guān)：

式（26）中，下標(biāo)中的尖括號表示該張量是無跡張量，即ξ〈iξj〉=ξiξj?ξ2δij/3。為避免使用過多符號，在線性化問題中，本文使用與原物理量相同的符號表示歸一化的擾動物理量，如這里的擾動數(shù)密度n應(yīng)理解為偏離參考數(shù)密度n0的那部分，再除以參考數(shù)密度；擾動溫度T為偏離參考溫度T0的部分，再除以參考溫度；在全局平衡態(tài)下，參考速度、應(yīng)力偏量和熱流全部為零，因此上式中它們分別以T0下 的最概然速率vm、n0kBT0和n0kBT0vm歸一化。

1.4 輸運(yùn)系數(shù)和弛豫率：表象與本質(zhì)

玻爾茲曼方程和納維-斯托克斯方程分別在介觀和宏觀尺度上描述氣體系統(tǒng)的演化，它們之間的關(guān)系一直是數(shù)學(xué)家和力學(xué)家研究的重要課題[63]。對玻爾茲曼方程求守恒量m、mc、mc2/2的矩，可以得到如下宏觀方程組：

其中E=3kBT/2m+u2/2。但是該方程組并不封閉，因?yàn)閴毫埩亢蜔崃鞯谋磉_(dá)式不能用密度、速度和溫度表示。有三種主要方法嘗試封閉宏觀方程組，它們分別是希爾伯特展開法[64]、Chapman-Enskog展開法[21,65-66]和矩方法[22,67]。這里介紹最常用的第二種方法。

在Chapman和Enskog[65-66]方法中，首先將玻爾茲曼方程改寫為：

然后將分布函數(shù)和玻爾茲曼碰撞項(xiàng)展開為 ε的冪級數(shù)：

值得注意的是，ε是一個小的形式參數(shù)，主要用于監(jiān)測展開的階數(shù)。展開完成后，取 ε=1。

同時，壓力和熱流也相應(yīng)展開如下：

但是，碰撞不變量對應(yīng)的宏觀量CM={ρ,u,T}僅由分布函數(shù)的零階展開決定，即：

而各個高階展開對CM的貢獻(xiàn)均為零，這稱為兼容性條件。

將式（30）代入式（27），可知時間導(dǎo)數(shù)亦可表示為ε級數(shù)[67]，

且對于零階近似，可得到：

從式（28）和式（29）易知Q(f(0))=0。因此，速度分布函數(shù)的零階展開為：

而速度分布函數(shù)的一階近似f(1)=Feqφ滿足以下方程[67]：

其中，方程的推導(dǎo)過程中用到

以及式（33）。注意式（21）中的feq應(yīng)改寫為Feq。

積分方程（35）的解可分為齊次部分與非齊次部分。齊次部分滿足J(φ)=0， 這說明 φ必然是碰撞不變量的線性組合，而兼容性條件要求所有線性組合的系數(shù)都為零。非齊次部分滿足如下形式：

其中，A(ξ) 和B(ξ)的解滿足：

一旦得到A(ξ)和B(ξ)，通過一階展開就可以恢復(fù)牛頓黏性定理與傅里葉熱傳導(dǎo)定理：

且剪切黏性與熱導(dǎo)率分別為：

一般將A(ξ) 和B(ξ)展開為索南多項(xiàng)式級數(shù)：

然后通過式（37）和索南多項(xiàng)式的正交性求解輸運(yùn)系數(shù)。對于麥克斯韋分子，取第一項(xiàng)展開可得到準(zhǔn)確的剪 切 黏 性 和 熱 導(dǎo) 率，即 取A(ξ)=a1(ξ2?5/2)和B(ξ)=b0；而對于其他相互作用勢，僅考慮第一項(xiàng)展開得到的輸運(yùn)系數(shù)的相對誤差也不超過2%。定義等效黏度截面為：

則剪切黏性為：

而熱導(dǎo)率為：

從而單原子氣體的普朗特?cái)?shù)為：

對于逆冪律分子間相互作用勢（式(11)），可以得出氣體黏性和溫度的關(guān)系：

其中 ω是黏性溫度冪指數(shù)。對于麥克斯韋分子和硬球分子，ω分別取1和0.5；對于其他氣體，ω一般介于0.5和1之間。

將式（38）代回式（27），即可獲得納維-斯托克斯方程。若繼續(xù)計(jì)算Chapman-Enskog展開的高階近似，可以得到Burnett、super-Burnett等宏觀方程組[21]。但是可能由于應(yīng)力和熱流與宏觀量CM展開的不匹配，導(dǎo)致高階方程組的線性不穩(wěn)定性。同時，由于邊界條件數(shù)量隨著方程組階數(shù)升高而增加，但對附加的邊界條件提法缺少充分的研究，因此高階方程組很少被實(shí)際應(yīng)用。另一方面，即使這些高階方程組在某些情況下能得到精確解，解的精度也不一定隨著展開階數(shù)的增加而增加[17]。

值得一提的是，對于麥克斯韋分子，在空間均勻系統(tǒng)中，本征值 λ02和 λ11與應(yīng)力偏量和熱流的弛豫速率相關(guān)，而這些弛豫速率決定著剪切黏性系數(shù)和熱傳導(dǎo)的大?。?/p>

可以認(rèn)為這些弛豫過程是本質(zhì)，而黏性系數(shù)和熱導(dǎo)率則為表象：在稀薄氣體流動中，本構(gòu)關(guān)系式（38）和等效輸運(yùn)系數(shù)會隨著克努森數(shù)的改變而改變，但一旦分子作用勢確定，弛豫系數(shù)就會確定不變。

利用弛豫時間和輸運(yùn)系數(shù)的關(guān)系，可以大大簡化模型方程中輸運(yùn)系數(shù)的推導(dǎo)過程。即，如果應(yīng)力偏量的弛豫時間為τ，則剪切黏性為pτ。如果熱流的弛豫時間為 τ/A， 則普朗特?cái)?shù)為A。

1.5 模型方程

玻爾茲曼方程碰撞項(xiàng)的復(fù)雜性給數(shù)值求解帶來極大挑戰(zhàn)[31]，因此，學(xué)者建議使用簡化的碰撞模型來代替。在簡化玻爾茲曼方程的碰撞項(xiàng)時，一般認(rèn)為需要滿足以下幾點(diǎn)：

1）動理學(xué)模型必須保證質(zhì)量、動量、能量守恒；

2）理想氣體的局部平衡態(tài)速度分布函數(shù)為麥克斯韋分布；

3）從模型方程導(dǎo)出的黏性系數(shù)與熱導(dǎo)率應(yīng)與玻爾茲曼方程導(dǎo)出的結(jié)果保持一致；

4）滿足熵增定理，當(dāng)且僅當(dāng)孤立系統(tǒng)處于平衡態(tài)時熵增為零，并達(dá)到其最大值。

其中，前兩項(xiàng)是基本要求，第三項(xiàng)要求在連續(xù)流區(qū)域通過Chapman-Enskog展開恢復(fù)出納維-斯托克斯方程。而第四項(xiàng)僅為從物理角度出發(fā)的考量，并不能保證模型方程的精度：如果一個動理學(xué)模型方程與玻爾茲曼方程完全一致，那么熵增速率也應(yīng)該相同，然而這通常是不可能做到的。實(shí)際上，大多數(shù)時候，僅在線性化情況下滿足熵增定理的Shakhov模型的精度優(yōu)于完全滿足熵增定理的ES-BGK模型。

由于玻爾茲曼碰撞項(xiàng)在形式上可記為Q=Q+?f/τ，因此，動理學(xué)模型方程的碰撞項(xiàng)一般采取如下形式：

其中：fr為參考速度分布函數(shù)；τ為特征碰撞弛豫時間，表征速度分布函數(shù)趨向參考分布函數(shù)的快慢，為簡單計(jì)，通常假設(shè)與分子速度無關(guān)。

1.5.1 BGK模型

BGK模型是最著名的氣體動理學(xué)簡化模型，由Bhatnagar、Gross和Krook提出[32]，其參考速度分布函數(shù)是局部麥克斯韋分布：

易見此模型滿足三大守恒律，且給出了在平衡態(tài)下的正確解f=Feq。在空間均勻系統(tǒng)中，由于 lnFeq是碰撞不變量的線性組合，有，所以可以證明BGK模型滿足熵增原理：

空間均勻系統(tǒng)中應(yīng)力偏量和熱流的弛豫時間均為τ，因此根據(jù)式（46）隨后段落的描述，可知剪切黏性為 μ=pτ，而普朗特?cái)?shù)為1。然而，通過玻爾茲曼方程推導(dǎo)得到的普朗特?cái)?shù)正確值為 2/3。這說明BGK模型不可能同時得到正確的剪切黏性和熱導(dǎo)率。對于黏性主導(dǎo)的流動，一般選取 τ以恢復(fù)剪切黏性；對于熱傳導(dǎo)主導(dǎo)的流動，則以恢復(fù)熱導(dǎo)率為目標(biāo)來選取τ。

1.5.2 ES-BGK模型

Holway[33]提出了ES-BGK模型以修正BGK模型中錯誤的普朗特?cái)?shù)，其中參考速度分布函數(shù)是在給定質(zhì)量、動量、能量和應(yīng)力張量條件下，通過最大化熵函數(shù)式（16）而來：

其中：

若b=0， 二階張量 λij中僅有主對角元素非零，模型恢復(fù)為BGK模型。

容易證明，應(yīng)力偏量的弛豫時間為 τ/(1?b)，而熱流的弛豫時間為τ。因此，剪切黏性和普朗特?cái)?shù)分別為：

為確保弛豫時間為正，則有b<1。 另外，有界的充要條件是b≥?1/2。因此，ES-BGK模型的普朗特?cái)?shù)范圍為 [2/3,+∞)。為了恢復(fù)單原子氣體的普朗特?cái)?shù)2/3，取

可以看出，當(dāng)碰撞項(xiàng)為零時，ES-BGK方程給出f=，當(dāng)且僅當(dāng) σij=0時，速度分布函數(shù)退化為麥克斯韋平衡態(tài)分布。而這并不違背動理學(xué)建模的第二個要求，因?yàn)楫?dāng)局部平衡態(tài)為f=時，應(yīng)力偏量的演化滿足如下形式 ?σij/?t=?pσij/μ，即應(yīng)力偏量將隨時間推進(jìn)逐漸衰減至零。因此，ES-BGK模型中趨于平衡態(tài)的過程可以解釋為：當(dāng)速度分布函數(shù)朝著橢球分布（式（50））松弛，橢球分布也朝著麥克斯韋分布演化，最終實(shí)現(xiàn)f=Feq。

Andries等[34]證明了ES-BGK模型滿足熵增定理，使之成為唯一滿足熵增定理且能正確恢復(fù)各輸運(yùn)系數(shù)的動理學(xué)模型，因此，ES-BGK模型受到廣泛關(guān)注。

1.5.3 Shakhov模型

不同于ES-BGK模型在參考速度分布函數(shù)中引入應(yīng)力修正項(xiàng)，Shakhov模型通過在參考速度分布函數(shù)中引入熱流修正項(xiàng)來恢復(fù)正確的輸運(yùn)系數(shù)：

此模型同樣滿足三大守恒律。從偏應(yīng)力和熱流的弛豫時間可以看出，該模型的黏性為μ =pτ，且能實(shí)現(xiàn)正確的普朗特?cái)?shù)。此外，由于熱流弛豫過程滿足?q/?t=?(2p/3μ)q， 在速度分布函數(shù)向frS靠近的過程中，參考速度分布函數(shù)也將逐漸演化為麥克斯韋分布，從而使Shakhov模型在平衡態(tài)時滿足f=Feq。

線性化的Shakhov模型滿足熵增定理。而對于非線性流動，熵增定理未被證明但也無從證偽；另外，速度分布函數(shù)可能出現(xiàn)非物理的負(fù)值。盡管如此，Shakhov模型仍能在許多問題中給出準(zhǔn)確結(jié)果，因此得到了廣泛應(yīng)用。

1.5.4 Gross-Jackson模型及其非線性化

針對麥克斯韋分子，根據(jù)線化玻爾茲曼碰撞項(xiàng)的本征值與本征函數(shù)，Gross和Jackson提出的以下形式的動理學(xué)模型去逼近線性化玻爾茲曼碰撞項(xiàng)J(φ)：

式中， λst是 一個任意負(fù)數(shù)， λnl是線性化玻爾茲曼項(xiàng)的本征值（式（24）），θ′為v和v′的夾角。該模型的物理意義在于保證模型方程和線性化玻爾茲曼方程各階矩的弛豫率相同。

Gross和Jackson將L(h)的N階近似表示為：

即保留了所有階數(shù)小于或等于N的多項(xiàng)式，并且取λst=λ0N，因此該截?cái)嗄Ｐ偷淖畲蟪谠ヂ蕿?λ0N。

在線性化情況下，即把速度分布函數(shù)在全局平衡態(tài)下根據(jù)式（20）展開，BGK、ES-BGK和Shakhov模型有如下形式：

其中擾動宏觀物理量的定義可參考式（26）。它們都可以看作是Gross-Jackson模型的特例。當(dāng)N=2時，從式（56）可以得到=LBGK；令N=3，根據(jù)式（24）有：

從而式（56）給出：

雖然，這既不符合線性化ES-BGK模型，也不符合線性化Shakhov模型，但是可將其視為這兩個模型的線性組合。另一方面，若給定約束 2n+l≤3，當(dāng)λst=?2p/3μ時，式（55）恰好展開為線性化ES-BGK模型；而當(dāng) λst=?p/μ時，式（55）展開為線性化Shakhov模型。

因此，通過Gross-Jackson模型和線性化BGK、ES-BGK、Shakhov模型的關(guān)系，可以反推出各種形式的非線性動理學(xué)模型。例如，陳松澤等將ES-BGK模型與Shakhov模型線性組合為一個復(fù)合模型，通過改變兩模型所占的比例和普朗特?cái)?shù)，復(fù)合模型能夠分別恢復(fù)到BGK、ES-BGK和Shakhov模型式[68]。類似地，吳雷在博士論文也提出一個帶有自由參數(shù)b的復(fù)合模型[69]：

該模型實(shí)際上是在Gross-Jackson模型（55）中取2n+l≤N=3和 λst=1/τ，并通過將如下的線性化碰撞項(xiàng)變?yōu)榉蔷€性碰撞項(xiàng)得來：

需要注意到，含有應(yīng)力偏量的項(xiàng)既可以被吸收到指數(shù)函數(shù)里使參考分布函數(shù)變成類似式（50）的形式，也可以放在外面（此時對b的 限制僅為b<1），兩者對數(shù)值計(jì)算的結(jié)果影響不大[69]。通過將Gross-Jackson模型中的其他階本征函數(shù)非線性化為相應(yīng)的形式，可以構(gòu)造更加準(zhǔn)確的動理學(xué)模型。

2 分子氣體的特性

玻爾茲曼方程是僅針對單原子氣體提出的，為模擬具有兩個或多個原子的分子結(jié)構(gòu)的氣體的動力學(xué)行為，則需要在氣體動理論層面重新建模。對于分子氣體，兩體碰撞保持動量守恒，但是能量卻可以在各形式自由度和能級之間進(jìn)行交換。因此，相比于僅具有剪切黏性的單原子氣體，分子氣體還存在體積黏性。另外，除了由分子平動引起的平動熱導(dǎo)率外，還有由分子轉(zhuǎn)動和振動引起的內(nèi)部熱導(dǎo)率。復(fù)雜的物理過程和更多的輸運(yùn)參數(shù)都給分子氣體動理論建模帶來挑戰(zhàn)。

分子的內(nèi)能可以表示為轉(zhuǎn)動能量、振動能量以及電離能之和，分子氣體內(nèi)部自由度導(dǎo)致了更加復(fù)雜的碰撞項(xiàng)。王承書和烏倫貝克從量子力學(xué)角度出發(fā)建立了WCU方程，定義fi(t,x,v) 為 第i個 能量為ei的能級的速度分布函數(shù)，并對每個能級的分布函數(shù)都寫出一個玻爾茲曼方程，其中遷移項(xiàng)保持不變，而碰撞項(xiàng)Qi是玻爾茲曼碰撞項(xiàng)的唯象擴(kuò)展，

碰撞后的分子速度為：

碰撞核和躍遷概率的具體表達(dá)式是非常復(fù)雜的。這里以N2為例，在只存在轉(zhuǎn)動自由度的情況下，對WCU方程的計(jì)算復(fù)雜性進(jìn)行說明。在Lennard-Jones作用勢下，轉(zhuǎn)動能可表達(dá)為：

而轉(zhuǎn)動能級的簡并度為：

對照分子動力學(xué)模擬數(shù)據(jù)[70-71]，將躍遷概率擬合為：

其中，P0為一常數(shù)，顯然，如果總能級數(shù)為N，式（61）的計(jì)算量為單原子玻爾茲曼碰撞項(xiàng)（式（7））的N4倍，計(jì)算時間成本巨大。

根據(jù)碰撞過程中平動動能是否守恒，WCU方程中的碰撞項(xiàng)可以分為彈性碰撞與非彈性碰撞。若處在i能 級和j能級的兩個分子發(fā)生碰撞，碰撞后仍然能分別處在i能 級和j能級，則為彈性碰撞，否則為非彈性碰撞：

對于彈性碰撞和非彈性碰撞，可以分別引入兩個等效的弛豫時間[72]：

而將WCU碰撞項(xiàng)形式地改寫為：

這兩個等效的碰撞時間，表征著能量在各種自由度或能級間的弛豫速度，因而與下面將討論的分子氣體輸運(yùn)系數(shù)、體積黏性、平動/轉(zhuǎn)動熱導(dǎo)率等，有著直接的關(guān)聯(lián)。

2.1 體積黏性

體積黏性表征了氣體無剪切膨脹或壓縮時所受的阻力，對于稀疏氣體，這種阻力來自于氣體分子平動動能與內(nèi)能的交換過程：在膨脹或者壓縮過程中，壓力做功會立刻改變平動動能，但在一定的延遲后才會通過非彈性碰撞來改變內(nèi)能。為方便起見，在下述推導(dǎo)中，假設(shè)只存在一個轉(zhuǎn)動能級。

引入非彈性碰撞所對應(yīng)的平動能與轉(zhuǎn)動能交換的弛豫時間 τr， 在空間均勻系統(tǒng)中，轉(zhuǎn)動溫度Tr的弛豫過程由Jeans-Landau方程描述：

式中，Tt和Tr分別為平動和轉(zhuǎn)動溫度，而平衡態(tài)溫度T可表達(dá)為：

忽略剪切黏性和熱傳導(dǎo)的影響，通過膨脹壓縮過程中能量守恒可得到：

式中， D/Dt為物質(zhì)導(dǎo)數(shù)，推導(dǎo)得到：

將Tt?Tr基 于小量 τr展 開為Tt?Tr=Ttτr+O(τ2r)，可知式（67）左側(cè)的時間導(dǎo)數(shù)項(xiàng)較右側(cè)為更高階小量，于是得到：

這說明，內(nèi)能弛豫的影響將壓力從原始壓力p0=nkBT改變?yōu)槠絼訅毫t=nkBTt。結(jié)合式（65）可得：

式中速度散度前的系數(shù)即為體積黏性 μb。在此，可以引入無量綱轉(zhuǎn)動碰撞數(shù)Z：

其中由分子平動引起的平均碰撞時間定義為τ=μ/n0kBTt。于是體積黏性與剪切黏性的比值有如下關(guān)系：

從式（70）和式（71）可知，體積黏性正比于能量交換時間。對于單原子氣體，若不考慮電離，不存在內(nèi)部能量的交換，因此體積黏性為零。而對于一些分子氣體，例如CO2，體積黏性可以是剪切黏性的上千倍[73]。但是，需要強(qiáng)調(diào)的是，式（68）成立的條件是τr很小，即能量弛豫的時間遠(yuǎn)遠(yuǎn)小于系統(tǒng)的特征時間。如果系統(tǒng)特征頻率 ?過高，則等效的體積黏性一般與系統(tǒng)頻率有如下近似關(guān)系[74-77]：

也就是說，在 ?τr→∞極限下，體積黏性趨向零而不是無窮大，這是因?yàn)槟芰拷粨Q時間過長以至于系統(tǒng)無法感知，此時可以認(rèn)為內(nèi)部自由度是凍結(jié)的。因此，與式（46）中提到的剪切黏性和熱導(dǎo)率的弛豫速率是本質(zhì)而剪切黏性和熱導(dǎo)率是表象一樣，這里Jeans-Landau方程里的能量交換也是本質(zhì)，而體積黏性只是連續(xù)流下的一個表象。

2.2 熱導(dǎo)率

分子氣體的熱流由兩部分組成，一部分是由分子的平動引起的動能傳遞，一部分是由分子的擴(kuò)散引起的內(nèi)能的傳遞。由于氣體分子遷移過程中存在動能和內(nèi)能交換，這一非彈性碰撞改變了能量在不同自由度下的分布，從而間接改變了動能傳遞與內(nèi)能隨擴(kuò)散傳遞的結(jié)果。因而，平動熱流和轉(zhuǎn)動熱流的弛豫過程是相互耦合的。于是，在空間均勻系統(tǒng)中，類似于式（46），平動熱流qt和 轉(zhuǎn)動熱流qr的弛豫過程一般可以寫成如下形式：

而通過Chapman-Enskog展開得到的連續(xù)流下的平動熱導(dǎo)率 κtr和 轉(zhuǎn)動熱導(dǎo)率 κrot可表示為：

將方程（74）代入方程（1），就可以通過熱弛豫系數(shù)求解如下方程，得到對應(yīng)的Eucken系數(shù)：

基于WCU方程（61），Manson和Manchick通過修正的Chapman-Enskog展開推導(dǎo)出四個無量綱熱弛豫系數(shù)的具體形式：

其中D′是有效擴(kuò)散系數(shù)，如果極性分子發(fā)生強(qiáng)共振碰撞，該系數(shù)可能與自擴(kuò)散系數(shù)D有明顯差異，例如HCl的有效擴(kuò)散系數(shù)滿足 ρD′/μ=0.89， 而 ρD/μ=1.33。根據(jù)式（75），忽略1 /Z的高階項(xiàng)，可將平動Eucken因子和轉(zhuǎn)動Eucken因子寫為[53]：注意到，ftr隨著轉(zhuǎn)動碰撞數(shù)Z的減小而減小，即內(nèi)能交換越頻繁，平動熱導(dǎo)率降低。這是因?yàn)樵诜肿託怏w的碰撞過程中，部分平動動能轉(zhuǎn)化為內(nèi)能，當(dāng)D′增加，Z減小時，能量轉(zhuǎn)化效率增強(qiáng)。這在物理上是合理的。

Mason與Monchick的解析解依然存在的一個問題，即無法保證通過式（76）推導(dǎo)而來的無量綱熱弛豫系數(shù)A與實(shí)際氣體參數(shù)保持一致。其原因在于，在黏性和擴(kuò)散系數(shù)確定的條件下，轉(zhuǎn)動碰撞數(shù)Z作為公式中唯一用來調(diào)節(jié)Eucken因子的參數(shù)，只能使得總熱導(dǎo)率可能達(dá)到與實(shí)驗(yàn)數(shù)據(jù)相近的數(shù)值，但無法分別獨(dú)立調(diào)整平動與內(nèi)能Eucken因子。與此同時，由于Z也是決定分子氣體體積黏性的唯一調(diào)節(jié)參數(shù)，參考式（71），則總熱導(dǎo)率與體積黏性的準(zhǔn)確性通常無法同時被滿足。

2.3 DSMC中的熱弛豫速率

對于單原子氣體，在粒子數(shù)和空間網(wǎng)格足夠多的情況下，DSMC已被證明收斂于玻爾茲曼方程的解[25]。實(shí)際上，由于DSMC中兩體碰撞后粒子速度的更新方式與玻爾茲曼方程一致，式（43）在可接受的誤差范圍內(nèi)恒成立，即DSMC可確保平動Eucken因子ftr非 常接近 2.5。但對于分子氣體，DSMC不一定能夠準(zhǔn)確恢復(fù)熱導(dǎo)率。在DSMC的碰撞模型中，可變軟球模型的提出是為了獲得準(zhǔn)確的剪切黏性與自擴(kuò)散系數(shù)；唯象的Larsen-Borgnakke模型[55]可以通過恢復(fù)能量交換速率[56-58]從而準(zhǔn)確獲取體積黏性。然而迄今為止，尚無研究證明DSMC能夠在分子氣體模擬中正確恢復(fù)熱導(dǎo)率。

我們通過提取DSMC中的熱弛豫速率A來 分析其恢復(fù)熱導(dǎo)率的能力。具體而言，在Larsen-Borgnakke模型中選取不同的轉(zhuǎn)動碰撞數(shù)，對非極性氣體N2與極性氣體HCl進(jìn)行研究。參照Bird書中的公式（4.50）和公式（5.67），DSMC中非彈性碰撞概率 Λ與本文中定義的轉(zhuǎn)動碰撞數(shù)Z通過Jeans-Landau方程（64）聯(lián)系在一起[24]：

式中，ω為黏性溫度冪指數(shù)，見式（45）；α是可變軟球模型中與斯密特?cái)?shù)相關(guān)的參數(shù)：

使用開源軟件SPARTA在空間均勻系統(tǒng)中提取DSMC方法中的熱弛豫系數(shù)矩陣A。將一百萬個DSMC模擬粒子均勻地分布在邊長為10 nm的立方體內(nèi)，所有的邊界設(shè)置為周期性邊界條件。氣體數(shù)密度為n0=2.69×1025m?3， 溫度為T0=300 K。如圖3（a）所示，在模擬初始時刻，沿x軸正方向運(yùn)動的粒子的速度分布函數(shù)為200 K溫度下的麥克斯韋分布，沿x軸負(fù)方向運(yùn)動的粒子速度分布函數(shù)為 400 K溫度下的麥克斯韋分布。類似地，沿x軸正方向運(yùn)動的粒子的轉(zhuǎn)動能滿足200 K溫度下的麥克斯韋分布，而沿反方向運(yùn)動的粒子轉(zhuǎn)動能滿足400 K溫度下的麥克斯韋分布，見圖3（b）。這樣的設(shè)置是為了建立初始平動和轉(zhuǎn)動熱流分布。隨后記錄平動和轉(zhuǎn)動熱流隨時間的演化，直至其在系統(tǒng)中按照熵增原理達(dá)到熱平衡后消失。圖3（c）展示的是通過100次系綜平均得到的平滑計(jì)算結(jié)果。

圖3 從DSMC模擬中提取熱流的弛豫速率：(a)初始狀態(tài)用于激發(fā)平動熱流的速度分布函數(shù)，橫坐標(biāo)由最概然速率v m歸一化； (b)初始狀態(tài)用于轉(zhuǎn)動熱流的轉(zhuǎn)動能分布函數(shù)，橫坐標(biāo)由k BT0歸一化；(c, d)氮?dú)馄絼訜崃骱娃D(zhuǎn)動熱流的弛豫過程，模擬參數(shù)為：式(78)中的非彈性碰撞概率Λ =0.25 ， 施密特?cái)?shù)S c=1/1.33， 轉(zhuǎn)動自由度d r=2 ， 溫度黏性指數(shù)ω =0.74. 圖中數(shù)據(jù)來源于文獻(xiàn)[60]Fig. 3 Extraction of heat flux relaxation rates from DSMC simulations. The initial distribution of (a) molecular velocity and (b) rotational energy of nitrogen molecules in DSMC, where the abscissas are normalized by and k BT0, respectively. (c, d) The evolution of heat fluxes and their time derivatives in Nitrogen. Λ =0.25 in Eq. (78), the Schmidt number is S c=1/1.33, the rotational degree of freedom is d r=2 , and the viscosity index is ω =0.74. The figures are modified from figure 1 in Ref. [60]

圖3（d）給出N2的熱流變化率的演化過程，可以看出平動熱流的時間導(dǎo)數(shù)隨時間先增加后減小，這可以解釋為平動熱流與轉(zhuǎn)動熱流之間存在強(qiáng)耦合效應(yīng)，即在平動熱流衰減的同時得到了轉(zhuǎn)動熱流的補(bǔ)償，且根據(jù)式（73）可知，Atr必定小于零。

圖4顯示通過最小二乘法求解線性回歸問題（式73）得到的熱流弛豫速率A?？梢钥闯?，在相同的轉(zhuǎn)動碰撞數(shù)Z和斯密特?cái)?shù)下，N2與HCl氣體的熱弛豫速率幾乎相同。另外，當(dāng)Z→∞ 時 ，有Att=2/3，Arr=Sc，Atr=Art=0，可以得到與式（3）一致的結(jié)果。圖中虛線為Mason和Manchick的結(jié)果（式（76）），由于忽略了1 /Z的高階項(xiàng)，其結(jié)果在1 /Z較大時與DSMC的結(jié)果出現(xiàn)明顯偏差。

圖4 根據(jù)理論公式(76)從DSMC模擬中提取的熱弛豫速率. 施密特?cái)?shù)為S c=1/1.33， 轉(zhuǎn)動自由度為d r=2; 對N2和HCl，溫度黏性指數(shù)ω 分別取0.74和1，圖中數(shù)據(jù)來源于文獻(xiàn)[60]Fig. 4 The extracted rates A in the relaxation of heat fluxes from the DSMC simulation, as per Eq. (76). The Schmidt number is Sc=1/1.33, and the rotational degree of freedom is d r=2; the viscosity index for Nitrogen and hydrogen chloride is ω = 0.74 and 1,respectively. The figures are modified from figure 1 in Ref. [60]

圖5給出根據(jù)式（74）和式（75）計(jì)算得到的Eucken式因子，并與常溫下的實(shí)驗(yàn)結(jié)果對比。對于N2，當(dāng)施密特?cái)?shù)取1 /1.33時 ，通過選擇特定的碰撞數(shù)1 /Λ≈3.5，DSMC可以恢復(fù)實(shí)驗(yàn)測量得到的總熱導(dǎo)率，但由于缺少相應(yīng)實(shí)驗(yàn)數(shù)據(jù)，我們無法分辨其是否能正確恢復(fù)平動與轉(zhuǎn)動熱導(dǎo)率分量。若此溫度下N2的體積黏性要求 1 /Λ?3.5，則DSMC給出的總熱導(dǎo)率必然存在誤差。對于極性氣體HCl，實(shí)驗(yàn)中常溫下測得的總Eucken因子約為1 .7，但是從圖5（a）中可見，如果取施密特?cái)?shù)為1 /1.33，從DSMC中導(dǎo)出的總Eucken因子最小只能達(dá)到1 .8。因此，無論如何選擇碰撞數(shù)，DSMC都不能恢復(fù)總熱導(dǎo)率。

圖5 室溫下Eucken因子的解析解(77)與DSMC數(shù)值解的對比. DSMC中使用可變軟球模型，自擴(kuò)散系數(shù) D取值為 D'，菱形、正方形、三角形分別表示平動、內(nèi)能和總Eucken因子，實(shí)心點(diǎn)為解析解，空心點(diǎn)為DSMC結(jié)果，虛線表示300 K溫度下實(shí)驗(yàn)測量的總Eucken因子，圖中數(shù)據(jù)來源于文獻(xiàn)[54]Fig. 5 The thermal conductivity obtained from the analytical solution (77) and the DSMC at room temperature, as a function of the inelastic collision probability Λ. The variable-soft-sphere model is used and the self-diffusion coefficient D in DSMC takes the value of D' . Filled diamonds, squares and triangles represent the translational, internal and total Eucken factors from Eq. (77),respectively, while the open symbols are the corresponding results from DSMC. Dashed lines show the total Eucken factor obtained from experiments at a temperature of 300 K.The figures are from Ref. [54]

從式（77）可得，當(dāng)氣體種類確定時，轉(zhuǎn)動碰撞數(shù)決定了體積黏性，因此可以通過調(diào)整DSMC中的擴(kuò)散系數(shù)以恢復(fù)總熱導(dǎo)率，見圖5（b）。然而，即便在體積黏性與總熱導(dǎo)率均與實(shí)驗(yàn)值相同的情況下，方程（77）的近似解與使用Larsen-Borgnakke模型的DSMC方法得到的平動熱導(dǎo)率ftr并不一致，且兩者的準(zhǔn)確性都尚不確定。對于HCl氣體，極性分子的共振相互作用[51]導(dǎo)致 ρD′/μ從 1 .33降 低至 0.89，為了恢復(fù)實(shí)驗(yàn)測量的總熱導(dǎo)率，從圖5（b）可知，DSMC的非彈性碰撞概率約等于1，達(dá)到了平動能與內(nèi)能交換的極限，因此并不是一個滿足物理要求的值。若選擇合理的非彈性碰撞概率，則將要求DSMC中的有效擴(kuò)散系數(shù)遠(yuǎn)小于自擴(kuò)散系數(shù)。這將帶來另一個問題：對于多組分氣體，擴(kuò)散系數(shù)與熱導(dǎo)率同等重要，而Larsen-Borgnakke模型難以同時正確恢復(fù)這兩個輸運(yùn)系數(shù)，亟需改進(jìn)。

3 分子氣體的動理論模型及精度

WCU方程的每一個內(nèi)能能級都對應(yīng)一個速度分布函數(shù)，因而計(jì)算代價巨大。對于強(qiáng)激波問題，N2的能級數(shù)可達(dá)到70，給數(shù)值求解帶來了極大的困難[71]。因此對WCU方程進(jìn)行簡化是非常必要的，而簡化的基本原則與單原子氣體相同，即需要滿足必要的守恒律和平衡態(tài)分布，還需要盡量使各物理量的輸運(yùn)系數(shù)與WCU方程保持一致，如擴(kuò)散系數(shù)、剪切黏性、體積黏性、熱導(dǎo)率及其在各自由度下的分量等。更接近本質(zhì)的說法是，各輸運(yùn)系數(shù)對應(yīng)的物理量的弛豫速率與實(shí)際氣體保持一致。學(xué)者們提出了許多描述分子氣體稀薄效應(yīng)的動理學(xué)模型方程，如Rykov模型和ES-BGK模型等。與在單原子氣體中Gross-Jackson模型可以視為所有模型的源頭類似，在分子氣體中可以認(rèn)為Hanson-Morse模型[78]是所有動理論模型的源頭。

3.1 Hanson-Morse線性化模型

此模型在線性化情況下逼近WCU方程（61）的碰撞項(xiàng)。這種構(gòu)造方法類似于單原子氣體的Gross-Jackson模型[79]，即通過對玻爾茲曼碰撞項(xiàng)進(jìn)行正交多項(xiàng)式展開，并保證各多項(xiàng)式對應(yīng)的宏觀量的弛豫速率與WCU方程一致而獲得。對于分子氣體，關(guān)于分子速度的多項(xiàng)式仍由麥克斯韋分子的本征函數(shù)（24）給出；關(guān)于內(nèi)部能量的多項(xiàng)式的前兩項(xiàng)取為1和ε=ei?eˉ， 其中eˉ為平均內(nèi)部能量[52]。

分子氣體的平衡態(tài)速度分布函數(shù)可表示為fi=oiFeq，其中

表征了內(nèi)能為ei的分子比例。在線性化問題中，通常考慮宏觀速度為零的全局平衡態(tài)，并將速度分布函數(shù)寫為如下形式：

其中每個內(nèi)部能級對應(yīng)的擾動分布函數(shù) φi滿足|φi|?1。注意到為與上文的單原子線性化模型保持一致，這里的速度和空間坐標(biāo)分別通過最概然速率vm和 特征尺寸 ? 進(jìn)行無量綱化，時間通過 ?/vm歸一化，加速度通過v2m/?無 量綱化，內(nèi)能ei通 過kBT0無量綱化。歸一化后的擾動物理量可以通過對擾動分布函數(shù)求矩得到：

若加速度為小量，WCU方程的線性化形式為：

其中Ji的形式與式（21）類似，不再贅述。根據(jù)Hanson-Morse模型[78]，該線性化碰撞項(xiàng)可近似為[78]：

Hanson-Morse模型假設(shè)各內(nèi)部能級的弛豫時間相同，則與式（76）對應(yīng)的弛豫率為：

在自發(fā)瑞利-布里淵散射中，如果模型中包含應(yīng)力偏量 σαβ，則該模型稱為Tenti-S7模型[80-81]，其缺點(diǎn)在于不能得到實(shí)驗(yàn)散射光譜。Tenti-S6模型通過移除應(yīng)力偏量 σαβ，即，令模型中的 σαβ=0，獲得與實(shí)驗(yàn)一致的散射光譜。實(shí)際上，如果不考慮內(nèi)部能級，Tenti-S7和Tenti-S6模型分別退化為單原子ES-BGK模型和Shakhov模型，而在大多數(shù)情況下Shakhov模型的精度高于ES-BGK模型。

當(dāng)在Hanson-Morse模型中移除應(yīng)力偏量，為了復(fù)現(xiàn)應(yīng)力偏量、能量弛豫、熱流的弛豫過程，式（84）中的J030應(yīng) 該改成 ? 1/τ 。同時，將系數(shù)J011修改為以下形式[80,82]，從而保證正確恢復(fù)熱導(dǎo)率。

3.2 Hanson-Morse模型的約化及非線性化

為了減小計(jì)算量，引入如下兩個約化速度分布函數(shù)：

消去內(nèi)能自由度，則Hanson-Morse模型可以簡化為：

其中：

至此，注意到單原子氣體中Gross-Jackson模型與BGK、ES-BGK、Shakhov等模型的關(guān)系，可以通過類似于式（60）中展示的非線性化程序構(gòu)造出分子氣體的動理學(xué)模型。即，首先寫出擾動量g和r對應(yīng)的速度分布函數(shù)，

并定義宏觀量如下：

總熱流定義為q=qt+qr，平動壓力、轉(zhuǎn)動壓力和總壓力分別定義為pt=nkTt、pr=nkTr和p=nkT。 注意此處G可以線性化為G=feq(1+g)， 而R應(yīng)線性化為R=(drkBT/2)feq(1+g+r)。其次，把速度分布函數(shù)的演化過程形式地寫為：

最后，根據(jù)式（88）和式（89）中的碰撞項(xiàng)，類比式（60），寫出四個參考分布函數(shù)Gt、Gr、Rt和Rr的形式。

例如，考慮Hanson-Morse模型中不出現(xiàn)應(yīng)力偏量。式（88）中的第一行和第二行即為線性化Shakhov模型的碰撞項(xiàng)，可以非線性化為：

結(jié)合非線性模型（式92）中的非彈性碰撞項(xiàng)(Gr?Gt)/Zτ，可以選取

以在線性化情況下退化為式（88）中的后兩項(xiàng)，其中能量交換由Feq(T)和Feq(Tt)的差距而來，熱流松弛由q′決定。

3.3 非線性模型

然而，歷史并沒有選擇非線性化路線，而是由許多學(xué)者各自提出了非線性動理學(xué)模型[33-34,50,61,83-87]。下面先介紹主要的Rykov模型和ES-BGK模型，并揭示與Hanson-Morse模型的內(nèi)在聯(lián)系，最后介紹吳模型。

3.3.1 Rykov模型

原始的Rykov模型[61]僅適用于無振動模態(tài)的雙原子氣體，但是構(gòu)造思想可以直接拓展至多原子分子氣體建模[83]。在該模型中，四個參考速度分布函數(shù)定義為：

選取彈性碰撞弛豫時間 τ=μ/pt以及適當(dāng)?shù)霓D(zhuǎn)動碰撞數(shù)Z，可保證恢復(fù)正確的剪切黏性和體積黏性。四個熱流弛豫系數(shù)為：

因此通過式（75）可以得到Eucken因子的表達(dá)式如下：

這說明，在Rykov模型中，可通過調(diào)整自由參數(shù) ω0和ω1修改Eucken系數(shù)。

當(dāng)Z→∞，該模型退化為單原子氣體的Shakhov模型。如果將其線性化，可以得到類似于第3.2節(jié)中的約化Hanson-Morse模型。不同之處在于空間均勻系統(tǒng)中，Rykov模型中平動和轉(zhuǎn)動熱流的弛豫不會相互干擾，這是因?yàn)槭剑?4）中Atr=Art=0。

3.3.2 ES-BGK模型

Holway通過最大熵原理提出了ES-BGK模型，該模型保留了標(biāo)準(zhǔn)BGK模型的數(shù)學(xué)簡便性，并且能夠恢復(fù)正確的普朗特?cái)?shù)，同時滿足熵增定理[33-34]。速度分布函數(shù)的演化方程為：

其中，弛豫溫度Trel定義為轉(zhuǎn)動溫度與總溫度的加權(quán)平均值：

當(dāng)Z→∞，轉(zhuǎn)動和平動能量交換消失，因此模型方程退化為單原子氣體的ES-BGK方程。

取彈性碰撞弛豫時間 τ=μ/ptPr及適當(dāng)?shù)霓D(zhuǎn)動碰撞數(shù)Z，則分子氣體的剪切和體積黏性自動滿足要求。其中普朗特?cái)?shù)可以通過調(diào)整b的取值得到：

為使參考分布函數(shù)有界且在物理上有意義，要求Z≥1和?1/2≤b≤1。 因此，普朗特?cái)?shù)范圍為2 /3≤Pr≤+∞。對于大部分雙原子分子來說，普朗特?cái)?shù)大約為 5 /7，因此可以保證要求。

與Rykov模型一樣，ES-BGK模型的平動熱流和轉(zhuǎn)動熱流的弛豫過程互不干涉；從其弛豫率可知，Eucken因子可以表達(dá)為：

注意到，對于多原子ES-BGK模型，平動Eucken因子與轉(zhuǎn)動Eucken因子的比值始終為 5/3，并且Eucken因子完全由普朗特?cái)?shù)決定。

3.3.3 吳模型

上述氣體動理學(xué)模型的弛豫時間 τ都與分子速度無關(guān)，與玻爾茲曼方程的弛豫時間不相符。因此，在正激波問題中，Rykov模型計(jì)算得到的結(jié)果可能會出現(xiàn)波前溫升的誤差，且隨著馬赫數(shù)的增加，這一偏離顯著增長。由此，為了提高模型方程準(zhǔn)確性，吳雷等首先將Rykov模型中的彈性碰撞項(xiàng) (Gt?G)/τ替換為單原子玻爾茲曼方程的碰撞項(xiàng)Q(G)，以使彈性碰撞弛豫時間是分子速度的函數(shù)[83]；然后，又考慮引入更為準(zhǔn)確的熱流弛豫過程[60]，從而在相同參數(shù)條件下，得到更為符合DSMC的模擬結(jié)果。

吳模型的演化方程為：

式中，R′t=(dr/2)kBTr(τQ+G)， 弛豫時間為 τ=μ/pt，對應(yīng)的參考速度分布函數(shù)為：

為了恢復(fù)熱流弛豫過程（式（73）），選取q′和q′′如下：

當(dāng)Z→∞，吳模型完全退化為單原子玻爾茲曼方程。另一方面，如果將玻爾茲曼碰撞項(xiàng)換回Shakhov碰撞模型，在線性化的情況下，吳模型則完全與第3.2節(jié)中的約化Hanson-Morse模型一致。

3.4 模型方程精度比較

本節(jié)以正激波結(jié)構(gòu)、庫特流、熱蠕動、二維熱驅(qū)動問題為例，比較不同模型描述稀薄氣體流動的精度。模擬中，分子氣體為N2，轉(zhuǎn)動自由度dr=2，模型中的相關(guān)參數(shù)選為Z=2.6671、ω=0.74，施密特?cái)?shù)Sc=1/1.33。 吳模型中的熱弛豫系數(shù)矩陣A來自第2.3節(jié)中的DSMC模擬，從而保證吳模型的Eucken因子為feu=1.993、ftr=2.365和frot=1.436，與DSMC保持一致；此外，進(jìn)一步假設(shè)這些參數(shù)不隨溫度變化。Rykov模型可以通過式（97）調(diào)整 ω0和 ω1恢復(fù)正確的Eucken系數(shù)。從式（102）可知，ES-BGK模型只有一個自由度來恢復(fù)總Eucken系數(shù)，即ES-BGK模型沒有分辨平動Eucken系數(shù)和轉(zhuǎn)動Eucken系數(shù)的機(jī)制。幸運(yùn)的是，對于N2，ES-BGK模型給出的平動和轉(zhuǎn)動熱導(dǎo)率與DSMC非常接近，見表1中關(guān)于Eucken參數(shù)的匯總。再次強(qiáng)調(diào)，ES-BGK和Rykov模型中平動和轉(zhuǎn)動熱流的弛豫過程是相互獨(dú)立的。

表1 不同動理學(xué)模型和DSMC模擬中所使用的Eucken因子與熱弛豫速率. N2的Eucken因子為feu=1.993Table 1 Eucken factors and heat flux relaxation rates used in different gas kinetic models and DSMC. The total Eucken factor of the Nitrogen at room temperature is feu=1.993

引入?yún)⒖济芏萵0、 參考溫度T0、最概然速率參考壓力n0kBT0、 參考熱流n0kBT0vm。宏觀量通過參考值進(jìn)行無量綱化，且無量綱克努森數(shù)定義為：

用于表征流動的稀薄程度，其中 ?為流動特征長度。

在以下各研究問題中，氣體在固體表面的反射均用完全漫反射作為邊界條件，即氣體反射速度分布為滿足固體表面宏觀量的麥克斯韋分布。在數(shù)值解法上，使用離散速度法求解動理學(xué)模型，快速譜方法求解玻爾茲曼碰撞項(xiàng)[31]，開源軟件SPARTA進(jìn)行DSMC模擬。

3.4.1 正激波

圖6給出馬赫數(shù)為7的一維N2正激波的結(jié)構(gòu)?？臻g坐標(biāo)用激波上游的分子平均自由程對其進(jìn)行無量綱化，密度與溫度通過 (Q?Qu)/(Qd?Qu) 歸 一化，其中Q=n、Tt、Tr，下標(biāo)u和d分別代表上游和下游的物理量。所有動理學(xué)模型的密度分布均與DSMC計(jì)算結(jié)果大體上相符，只有ES-BGK模型給出的密度在波前稍低一些。Shakhov模型和ES-BGK模型的平動溫度在波前存在“過早溫升”現(xiàn)象，其中ES-BGK的平動與轉(zhuǎn)動溫度分布偏離最大，這是因?yàn)檫@兩個模型方程的碰撞時間與分子速度無關(guān)。而吳模型和DSMC中的碰撞時間都是分子速度的單調(diào)遞增函數(shù)，因而得到最為接近的結(jié)果。圖6（c、d）給出了壓力與熱流分布的計(jì)算結(jié)果，可見吳模型與DSMC符合得更好，而其他模型均存在顯著誤差，尤其是熱流分布的對比。

圖6 不同動理學(xué)模型預(yù)測的馬赫數(shù)為7的氮?dú)庹げǖ拿芏?、溫度、?yīng)力偏量σ22和熱流的對比Fig. 6 Comparisons in the density, temperature, stress, and heat flux between different kinetic models for a normal shock in Nitrogen gas with Mach number 7

3.4.2 庫特流動

平板庫特流動的上下兩板溫度相同，處在x2=L/2的 上板速度是vw=vm， 處在x2=?L/2的下板速度相反。圖7給出各模型和DSMC的對比結(jié)果。當(dāng)Kn=0.1時，所有動理學(xué)模型都給出了良好的密度分布。ES-BGK模型高估了平動溫度和平動熱流，卻低估了轉(zhuǎn)動熱流。Rykov模型給出良好的溫度和轉(zhuǎn)動熱流，卻高估了平動熱流。吳模型與DSMC的對比結(jié)果非常理想。當(dāng)Kn=1時，吳模型同樣能給出與DSMC較為一致的結(jié)果。Rykov模型方程的密度分布與DSMC的結(jié)果稍有差別，在兩板中心區(qū)域略低。ES-BGK模型的轉(zhuǎn)動溫度也略低于其他模型。

圖7 平板庫特流動中的密度、溫度和垂直于流動方向的熱流分布. 第一行和第二行的克努森數(shù)分別為K n=0.1、1； 基于對稱性，另一半?yún)^(qū)域? 0.5≤x1≤0沒有顯示Fig. 7 Profiles of the density, temperature, and heat flux (perpendicular to the plates) in the planar Couette flow of the Knudsen number Kn = 0.1 (first row) and 1 (second row). Due to symmetry, the other half region ? 0.5≤x1≤0 is not shown

3.4.3 麥克斯韋妖驅(qū)動的微流動

考察處在兩平行靜止的無限長平板間的分子氣體，每個氣體分子受到與其速度相關(guān)的水平外力的作用：

在此例中，取 2a0L/v2m=0.0718， 以保證系統(tǒng)僅略微偏離平衡態(tài)。可以看出，當(dāng)分子速率大于時，其加速度為正，反之則為負(fù)。因此該加速度可以看出是一種特殊形式的麥克斯韋妖，見圖8。圖8給出克努森數(shù)分別為0.1和1的模擬結(jié)果對比。吳模型與DSMC非常符合，速度剖面與熱流分布的最大誤差僅為 2%。值得注意的是，ES-BGK模型和Rykov模型的轉(zhuǎn)動熱流值遠(yuǎn)大于吳模型和DSMC的結(jié)果。其原因在于，這兩個模型方程中的平動熱流弛豫與轉(zhuǎn)動熱流弛豫相解耦，即Atr=Art=0，即便熱導(dǎo)率相近，但熱弛豫系數(shù)的差別仍可導(dǎo)致顯著的宏觀量的差異。這正說明了采用正確的熱弛豫速率的重要性。

圖8 麥克斯韋妖驅(qū)動下的稀薄氣體熱蠕動, 其中速度和熱流均沿 x 1方向. 第二行和第三行圖對應(yīng)的克努森數(shù)分別為Kn = 0.1和1，基于對稱性，0 .5≤x2≤1的區(qū)域沒有顯示Fig. 8 Profiles of the vertical velocity and heat flux in the thermal creep flow driven by the Maxwell demon, where the Knudsen number in the second and third rows is Kn = 0.1 and 1, respectively. Due to symmetry, the other half region 0 .5≤x2≤1 is not shown

3.4.4 封閉方腔中的熱蠕動

最后，考察二維封閉方腔內(nèi)的熱蠕動問題。方腔寬度與高度比值為5，固定左板和右板的溫度分別為200 K和400 K，而上下兩板的溫度在此間線性分布。由于此問題以x2=0.5為軸具有對稱性，實(shí)際計(jì)算區(qū)域?yàn)榉角坏南掳氩糠?。取參考溫度?00 K，參考長度為方腔高度，模擬克努森數(shù)為Kn=0.6的氣體流動。

圖9中給出各個動理學(xué)模型方程與DSMC結(jié)果的比較。吳模型的熱流云圖與速度剖面均與DSMC結(jié)果相符，而ES-BGK模型熱流與速度剖面的解與DSMC的解有非常大的偏離，并且ES-BGK模型計(jì)算得到的中心渦更偏向右側(cè)的高溫區(qū)，且結(jié)構(gòu)范圍更小。另外，Rykov模型在預(yù)測熱流方面存在一定誤差。

4 熱弛豫速率不確定度量化

為了準(zhǔn)確地預(yù)測稀薄分子氣體的非平衡行為，氣體動理學(xué)的建模過程應(yīng)實(shí)現(xiàn)正確的剪切黏性、體積黏性、平動熱導(dǎo)率、轉(zhuǎn)動熱導(dǎo)率。分子氣體的熱導(dǎo)率比剪切黏性和體積黏性更加復(fù)雜，實(shí)驗(yàn)通常只能確定總熱導(dǎo)率，但是難以分辨其平動與轉(zhuǎn)動分量。雖然在連續(xù)流極限下，氣體動力學(xué)行為取決于氣體黏性與總熱導(dǎo)率，但是在稀薄流區(qū)，分子氣體的平動與內(nèi)能熱導(dǎo)率可能對非平衡效應(yīng)有各自不同的影響。

上文的算例表明，熱弛豫速率在稀薄氣體流動中也起著重要作用。盡管在DSMC和其他分子氣體動理學(xué)模型中[33,72,88-89]，已有大量關(guān)于溫度弛豫效應(yīng)（與體積黏性相關(guān)）的研究工作，卻從未有過關(guān)于熱弛豫速率矩陣A對稀薄氣體非平衡效應(yīng)的影響的相關(guān)研究。通過實(shí)驗(yàn)可以直接測量總熱導(dǎo)率，并且在某些情況下，還可以提取平動熱導(dǎo)率[54,59,90-91]。然而，根據(jù)式（73），弛豫系數(shù)矩陣A中至少有兩個元素不能確定，即不同的系數(shù)矩陣A可以對應(yīng)一致的熱導(dǎo)率及其分量。因此有必要量化在輸運(yùn)系數(shù)都一致的前提下，由系數(shù)A的變化所帶來的宏觀量的不確定性。

在DSMC中，因碰撞模型的原因，固定體積黏性后，熱弛豫系數(shù)也隨之確定。在吳模型中，可以通過獨(dú)立修改Aij以量化熱弛豫系數(shù)的影響。上文已驗(yàn)證吳模型的準(zhǔn)確性，這里使用吳模型來量化這種不確定性。首先，在約束ftr和frot不變的前提下修改Aij，以研究熱弛豫速率的影響。其次，在保持總熱導(dǎo)率不變前提下，通過改變ftr和frot來量化熱導(dǎo)率分量的影響。

4.1 正激波

給定ftr和frot，即體系中所有的熱導(dǎo)率分量都被確定，若此時其余的輸運(yùn)系數(shù)也保持不變，則稱此系統(tǒng)為宏觀意義上的確定系統(tǒng)。然而，不同的Aij的選擇將可能導(dǎo)致結(jié)果的不確定性。

考察一維正激波問題，Atr、Art的取值范圍分別為[?5/(6Z),0]和 [?1/(3Z),0]， 而Att和Arr的值根據(jù)式（73）調(diào)整，以保持ftr和frot不變。在DSMC中非對角線元素 為Atr=?0.201,Art=?0.059，本 例 中，給 定Z=2.6671，Art和Atr的 最小值分別為? 0.3124 和 ? 0.1250，數(shù)值大小約為DSMC對應(yīng)參數(shù)的一到兩倍。研究發(fā)現(xiàn)，熱弛豫系數(shù)A的變化導(dǎo)致轉(zhuǎn)動溫度和轉(zhuǎn)動熱流稍有變化，而對密度、速度和應(yīng)力偏量幾乎沒有影響。

接著，保持總熱導(dǎo)率不變，且選擇Atr=?5/(6Z),Art=?1/(3Z)，研究熱導(dǎo)率分量的改變對流動的影響。圖10展示了馬赫數(shù)為4時不同平動Eucken因子下吳模型的數(shù)值結(jié)果。在實(shí)際中，較小數(shù)值的平動Eucken因子ftr是可能存在的，特別是在極性氣體中，真實(shí)的平動Eucken因子可能遠(yuǎn)小于 2.5。例如，水（H2O）的平動Eucken因子為ftr=1.78，甲醇分子（CH3OH）[53]的平動Eucken因子為ftr=0.41。對比可見，所有的宏觀量都隨ftr改變而有明顯的變化。首先，較大的ftr使得平動溫度更快上升到其最大值，并在波后更快趨于下游溫度，對于應(yīng)力偏量和總熱流也有相同的趨勢。其次，ftr變化所導(dǎo)致的轉(zhuǎn)動溫度變化則集中于激波結(jié)構(gòu)的中心處，較小的ftr導(dǎo)致更高的轉(zhuǎn)動溫度。最后，較大的ftr使得密度升高更快。

圖10 總Eucken因子不變時，不同平動Eucken因子對馬赫數(shù)為4的正激波結(jié)構(gòu)的影響. 圖中數(shù)據(jù)來源于文獻(xiàn)[60]Fig. 10 Influence of the translational Eucken factor on the structure of normal shock wave with Mach number 4. The figures are from Ref. [60]

由此，在正激波問題中，由熱弛豫系數(shù)A的變化引起的不確定性較小，而熱導(dǎo)率分量的變化則會帶來顯著的不確定性。例如，ftr=1.5的轉(zhuǎn)動熱流大約是ftr=2.5的兩倍，幾乎與平動熱流相當(dāng)。

4.2 麥克斯韋妖驅(qū)動的微流動

使用與正激波算例中相同的系數(shù)矩陣A，使ftr和frot均保持不變，來考察麥克斯韋妖驅(qū)動的微流動中，熱弛豫系數(shù)對流動速度和熱流的影響。圖11第一行圖片數(shù)據(jù)為Kn=0.2的結(jié)果?？梢钥闯?，不同于正激波問題，此問題中不同熱弛豫系數(shù)A對應(yīng)的流動速度與熱流有相當(dāng)明顯的差別。流動速度和平動熱流的最大相對不確定性分別為1 6.7%和 1 7.6%。同時可以看出，流動的中間區(qū)域出現(xiàn)了較大的不確定性，而壁面附近的速度滑移與熱流幾乎不變。

為進(jìn)一步研究平動Eucken因子的影響，考察Eucken因子ftr=1.5,2.0,2.5的算例，其他參數(shù)為Kn=0.2，fu=1.993，Atr=?5/6Z，Art=?1/3Z。如圖11第二行圖片數(shù)據(jù)所示，速度曲線與平動熱流隨著平動Eucken因子的變化而劇烈變化。ftr=2.5時的流動速度和熱流比ftr=1.5時 增加了大約 68%。在第一行圖中，在確定的ftr下，壁面速度滑移與熱流不受平動Eucken因子的影響。而在第二行圖中，速度與熱流依賴于平動Eucken因子，例如壁面附近的速度與熱流均與ftr呈正相關(guān)。這說明，在此問題中，平動Eucken因子ftr起主導(dǎo)作用。

圖11 (第一行)熱流弛豫速率對麥克斯韋妖驅(qū)動的熱蠕動速度和熱流分布的影響，以及吳模型中改變熱弛豫速率的結(jié)果對比. 紅色實(shí)線對應(yīng)的 A數(shù) 值來自于DSMC，藍(lán)色陰影部分來自吳模型計(jì)算結(jié)果，對應(yīng)的取值范圍為 A rt∈[?0.3124,0.0]、 A tr∈[?0.1250,0.0]，其他參數(shù)為Kn=0.2 , ftr=2.365 以 及 frot=1.435； (第二行)總Eucken因子 feu保 持不變，改變平動Eucken因子的結(jié)果對比，K n=0.2.圖中數(shù)據(jù)來源于文獻(xiàn)[60]Fig. 11 (First row) Influence of the thermal relaxation rates in the creep flow driven by the Maxwell demon, and comparison of the results of varying the thermal relaxation rate in the Wu model. Red solid lines are the results with A extracted from DSMC, blue shade region shows the results from the kinetic model (103), with A rt∈[?0.3124,0.0] and A tr∈[?0.1250,0.0]. Other parameters are K n=0.2,ftr=2.365 and frot=1.435. (Second row) Influence of the translational Eucken factor, while feu is fixed. The Knudsen number is K n=0.2.The figures are from Ref. [60]

平動Eucken因子在此問題中的重要性可得到如下解釋：從方程（73）可以看出，熱弛豫系數(shù)Atr和Art與平動能與轉(zhuǎn)動能之間的能量交換有關(guān)。當(dāng)Atr與Art為零時，平動熱流與轉(zhuǎn)動熱流的弛豫過程互相解耦。例如，更多的計(jì)算結(jié)果顯示，在Art=0的情況下，轉(zhuǎn)動熱流恒為零。這是因?yàn)辂溈怂鬼f妖驅(qū)動的熱蠕動流動中，外力項(xiàng)僅能直接影響平動能，而轉(zhuǎn)動能和熱流則需通過能量交換才能發(fā)生改變，具體取決于Atr、Art和Z。由于矩陣的非對角元素相對于其他元素為小量，這導(dǎo)致qrot≈0 ，而平動效應(yīng)（qtr或者ftr）占主導(dǎo)。

5 確定熱流弛豫速率的方法

上文的研究結(jié)果揭示了熱流弛豫速率不確定性對稀薄氣體流動的影響，因此有必要通過獲取并使用正確的相應(yīng)系數(shù)，來減少或消除模型的不確定性。我們認(rèn)為可以嘗試以下兩種方法。

5.1 分子動力學(xué)模擬

從物理上來說，只要知道了分子間的相互作用勢，根據(jù)牛頓第二運(yùn)動定律，就可以通過分子動力學(xué)方法模擬氣體的行為。由此，輸運(yùn)系數(shù)能夠通過Green-Kubo公式獲得[92-93]：

其中，V和T分別是模擬系綜的體積和溫度，τ為時間，尖括號表示系綜平均；Jp=mcici?pVI，Jq=Eici，I是單位矩陣。對于單原子氣體，Ei=mcici/2得到平動熱導(dǎo)率。對于分子氣體，取Ei=εi為 第i個分子的內(nèi)部能量，則式（110）給出內(nèi)部熱導(dǎo)率。

然而，對于如何獲得熱流弛豫矩陣系數(shù)尚不清楚。一個可行的方法是仿照第2.3節(jié)中介紹的DSMC方法，即在分子動力學(xué)中設(shè)置帶有初始熱流的狀態(tài)，模擬低密度氣體的演化過程，然后通過擬合公式（73）獲得弛豫系數(shù)。雖然目前尚未有相關(guān)工作的報(bào)道，但是原理上是可行的。

5.2 瑞利-布里淵散射

瑞利-布里淵散射指的是光在氣體中傳播時，被氣體分子的自發(fā)密度漲落散射。由于散射光譜包含氣體信息，因此可以用來無損測量聲速、流速、溫度和體積黏性[94-96]。例如，歐洲航天局于2018年發(fā)射的ADM-Aeolus衛(wèi)星，測量從地球表面到30 km高度的全球風(fēng)速，用于提高天氣預(yù)報(bào)精度。

如圖12所示，波矢為ki的入射光被散射，散射角為θ，則散射波數(shù)為k=|ki?ks|=2|ki|sin(θ/2)。散射光強(qiáng)度與光譜密度函數(shù)S(fs)相關(guān)，即密度擾動相關(guān)函數(shù)G(r,t) 的 時空傅里葉變換。假設(shè)散射波沿x2方向傳播，相關(guān)函數(shù)和譜密度函數(shù)的解[97-98]為：

圖12 自發(fā)瑞利-布里淵散射示意圖[54]，在氣體中傳播的光被氣體分子自發(fā)的密度漲落散射Fig. 12 Schematic of the spontaneous Rayleigh-Brillouin scattering[54], where the light is scattered by the spontaneous density fluctuations in the gas

因此，頻譜為：

式中，i為虛數(shù)單位，n(x2,t)為 密度擾動，fs為散射頻率。

對于短波長和高散射頻率，應(yīng)該使用氣體動理論來描述密度擾動的演化過程。以單原子氣體為例說明瑞利-布里淵散射的計(jì)算。密度擾動量與式（20）中分布函數(shù)擾動量 φ有關(guān)，初始密度脈動可定義為[83,97]：φ(t=0,x2,ξ)∝δ(x2)。 取特征長度為散射波長 ?=2π/k，通過時間上的拉普拉斯變換與空間上的傅里葉變換

求解線性化玻爾茲曼方程（21），可得：

其中方程右端的1來自于初始密度擾動對應(yīng)的拉普拉斯變換，(v) 分 別對應(yīng) φ的時空坐標(biāo)下的拉普拉斯和傅里葉變換[17]。散射頻譜可通過以下公式計(jì)算：

其中 ? 表示復(fù)數(shù)的實(shí)部。

此例中，我們根據(jù)荷蘭自由大學(xué)的實(shí)驗(yàn)數(shù)據(jù)[91]來提取N2的體積黏性和平動Eucken因子。其中，激光波長為 366.8 nm，散射角為90°，因此等效波長為?=259.4 nm。實(shí)驗(yàn)條件下N2的轉(zhuǎn)動自由度與振動自由度分別為dr=2和0。剪切黏性和熱導(dǎo)率根據(jù)Sutherland公式計(jì)算。通過調(diào)整吳模型中非彈性碰撞數(shù)和平動熱導(dǎo)率，使得實(shí)驗(yàn)數(shù)據(jù)和數(shù)值模擬得到的頻譜偏差最小。表2給出了不同溫度下氣體物性參數(shù)和提取的碰撞數(shù)與熱導(dǎo)率分量。

表2 自發(fā)瑞利-布里淵散射實(shí)驗(yàn)中，N2相關(guān)參數(shù)匯總以及根據(jù)吳模型從實(shí)驗(yàn)數(shù)據(jù)中提取的體積黏性與平動/轉(zhuǎn)動Eucken因子(剪切黏性與總熱導(dǎo)率均使用國際單位制)Table 2 Experimental conditions in the spontaneous Rayleigh-Brillouin scattering, and the extracted bulk viscosity,translational/rotational Eucken factors from the Wu model. The SI units are adopted for the shear viscosity and thermal conductivity

實(shí)驗(yàn)數(shù)據(jù)與吳模型最佳預(yù)測的比較如圖13所示。實(shí)驗(yàn)[99]測得N2在標(biāo)準(zhǔn)溫度和壓力下的轉(zhuǎn)動弛豫時間為 τr=7.4×10?10s。在 296.7 K下，測得的N2體積黏性為 μb=1.22×10?5kg·m?1·s?1。根據(jù)式（71），計(jì)算得到N2的體積黏性為：

圖13 基于實(shí)驗(yàn)頻譜(圓圈)和吳模型預(yù)測結(jié)果(線條)，提取N2的體積黏性與平動Eucken因子[54]：第一行的實(shí)線為計(jì)算光譜，通過計(jì)算給定 ftr 和 Z 范圍內(nèi)的最小誤差；第二行給出了實(shí)驗(yàn)光譜與理論光譜之間的誤差.實(shí)驗(yàn)與計(jì)算中所用到的參數(shù)匯總在表2中，頻率 fs由 v m/?歸一化Fig. 13 Extraction of the bulk viscosity and translational Eucken factor of Nitrogen from the experimental spectra (circles)[54]. Lines in the first row show the spectra obtained from the Wu model, while those in the second row show the corresponding residuals between the experimental and theoretical spectra. Experimental conditions and extracted parameters are summarized in Table 2. The frequency fs is normalized by vm/?

這一結(jié)果與實(shí)驗(yàn)值吻合較好。另外，提取的平動Eucken因子也在合理的范圍。

6 總結(jié)與展望

本文回顧了國內(nèi)外學(xué)者在稀薄效應(yīng)動理學(xué)建模中所取得的進(jìn)展，包括單原子氣體與分子氣體。由于高溫真實(shí)氣體效應(yīng)，當(dāng)分子內(nèi)自由度增加時，分子間復(fù)雜的相互作用通常伴隨著更加新穎、復(fù)雜的物理現(xiàn)象，這對稀薄分子氣體動理學(xué)建模提出了很高的要求。我們指出，應(yīng)力偏量、能量交換、熱流的弛豫過程是輸運(yùn)現(xiàn)象的本質(zhì)，而剪切黏性、體積黏性、平動/內(nèi)部熱導(dǎo)率等輸運(yùn)系數(shù)則是表象，即，隨著克努森數(shù)的變化，這些輸運(yùn)系數(shù)對應(yīng)的本構(gòu)關(guān)系會失效，而相應(yīng)的弛豫過程和速率不會改變。為了準(zhǔn)確預(yù)測稀薄分子氣體的非平衡行為，氣體動理學(xué)的建模過程應(yīng)盡可能準(zhǔn)確地捕捉應(yīng)力偏量、能量交換、熱流的弛豫過程和速率。

本文還介紹了DSMC中提取熱弛豫速率的方法，并討論了DSMC方法恢復(fù)熱導(dǎo)率的能力。雖然經(jīng)典的Larsen-Borgnakke模型可以通過調(diào)整非彈性碰撞概率改變體積黏性，但是很難同時調(diào)整熱導(dǎo)率及其分量。而這些分量的比例以及本質(zhì)的熱流弛豫時間對稀薄氣體流動的影響巨大。因此，非常有必要對DSMC碰撞模型進(jìn)行修改和優(yōu)化。

非常有意思的是，基于王承書等的工作，分子氣體的線性化理論和方法已經(jīng)非常成熟。而非線性模型雖然五花八門，但都與線性化模型存在著本質(zhì)聯(lián)系，即所有的非線性模型都可以通過對線性化Gross-Jackson模型和Hanson-Morse模型的非線性化而來。本文首次給出一個非線性化方案。

對于臨近空間高超聲速飛行器而言，其周圍的空氣動力學(xué)非常復(fù)雜。如何建立描述多組分氣體、化學(xué)反應(yīng)的稀薄氣體動理學(xué)模型是個極大的挑戰(zhàn)。本文提出的弛豫過程本質(zhì)論、輸運(yùn)系數(shù)表象論、從線性化模型中反推非線性模型的方法，以及從分子動力學(xué)和相關(guān)實(shí)驗(yàn)中提取弛豫速率的設(shè)想，將為解決此類問題提供極大助力。

致謝：感謝愛丁堡大學(xué)蘇微博士的有益討論。