楊鯉銘，李志輝，舒 昌,5

（1.南京航空航天大學(xué)航空學(xué)院，南京 210016；2.南京航空航天大學(xué)非定?？諝鈩恿W(xué)與流動控制工業(yè)和信息化部重點實驗室，南京 210016；3.中國空氣動力研究與發(fā)展中心超高速空氣動力研究所，綿陽 621000；4.國家計算流體力學(xué)實驗室，北京 100191；5.新加坡國立大學(xué)機械工程系，新加坡 119260）

飛船返回艙、HTV-2、東風(fēng)-17 等往返大氣層與空間軌道跨流域飛行器（圖1）的研究，一直以來都備受各軍事和航天強國的關(guān)注，在載人登月、深空探測和國防領(lǐng)域有著重要的應(yīng)用價值。然而，由于這類飛行器的飛行高度范圍非常大，經(jīng)常會跨越稀薄流和連續(xù)流整個流域范圍，對試驗研究和數(shù)值模擬都提出了嚴(yán)苛的考驗［1-5］。一方面，由于這類飛行器所需要考慮的相似參數(shù)眾多，流動相似準(zhǔn)則要求苛刻，地面風(fēng)洞實驗設(shè)備對于同時復(fù)現(xiàn)其所處的低雷諾數(shù)與高焓非平衡流動特征通常無能為力，而且自由飛試驗的代價極大，數(shù)據(jù)采集也非常困難。另一方面，由于飛行軌跡跨越不同的流域，對飛行器的流動模擬需要綜合考慮全流域范圍的計算精度和計算效率。此外，即便在某一環(huán)境克努森數(shù)下，飛行器周圍流場的局部克努森數(shù)也會跨越多個數(shù)量級，同時包含連續(xù)流、滑移流、過渡區(qū)域的高度稀薄流，致使該問題很難被準(zhǔn)確高效地求解［1-2，6-8］。

圖1 往返大氣層與空間軌道跨流域飛行器示意圖（圖片源于網(wǎng)絡(luò)）Fig.1 Schematic diagram of flight vehicles covering different flow regimes between the atmosphere and space orbit(pics.from the network)

相較于實驗研究，數(shù)值模擬無疑是解決該類問題的一種經(jīng)濟有效的手段。但由于稀薄效應(yīng)的存在，傳統(tǒng)的基于連續(xù)性假設(shè)的Navier-Stokes 方程不再適用于該類問題［9］。為了有效地研究該類問題，通常需要借助于從分子動理學(xué)理論（分子運動論）發(fā)展而來的Boltzmann 方程。該方程是位置空間、速度空間和時間的七維多相空間幾率密度分布函數(shù)方程［10］，同時其碰撞積分項是一個結(jié)構(gòu)復(fù)雜未有明確表達式的五重積分，因而直接求解極為困難?；谠摲匠蹋ㄟ^合理簡化，目前發(fā)展了兩類典型的數(shù)值求解算法。第一類是粒子方法，它通過模擬假想粒子的遷移和碰撞過程來實現(xiàn)對流場的仿真，建模過程等價于對Boltzmann 方程的輸運項和碰撞項進行解耦處理。典型代表為直接模擬蒙特卡羅（Direct simulation Monte Carlo，DSMC）方法［4，8，11-13］。第二類是離散速度或離散坐標(biāo)法（Discrete velocity or discrete ordinate method，DVM 或DOM［14-18］），其在位置空間和速度空間同時離散求解Boltzmann 方程，并且為了避免對碰撞積分不定式的計算，通常采用簡化模型來近似碰撞積分項。

DSMC 方法是現(xiàn)階段解決高超聲速飛行器在稀薄流域最有效實用的粒子類方法，由Bird［19］首次提出。該方法的關(guān)鍵是在小于當(dāng)?shù)胤肿悠骄鲎矔r間步長Δt內(nèi)將分子運動和碰撞解耦［3-4，8-9］，即讓所有分子運動一定距離并考慮邊界反射，然后計算此時間段內(nèi)具有代表性的分子間碰撞。在碰撞取樣數(shù)趨于無窮大時，DSMC 方法模擬得到的速度分布函數(shù)收斂到Boltzmann 方程的修正形式［20］。所以，為了確保模擬不失真，通常要求DSMC 方法的網(wǎng)格尺寸小于分子平均自由程，且時間步長小于分子平均碰撞時間。在稀薄流區(qū)域，由于分子數(shù)較少，分子平均自由程和平均碰撞時間較大，DSMC方法具有較高的計算效率和計算精度。但在連續(xù)和近連續(xù)流區(qū)域，由于分子平均自由程和平均碰撞時間較小，DSMC 方法的計算效率受到了極大的限制。此外DSMC 方法還存在統(tǒng)計噪聲問題，以及不便應(yīng)用于非定常流動模擬和不易于采用隱式算法加速等困難。為了應(yīng)用于跨流域問題的求解，一些學(xué) 者 提 出 了Navier-Stokes/DSMC 耦 合 算 法［21-23］，其思想是將計算域劃分為不同的子區(qū)域分別應(yīng)用DSMC 方法和Navier-Stokes 方程求解器進行求解。雖然該類耦合算法能在各自的區(qū)域獲得準(zhǔn)確高效的計算結(jié)果，但區(qū)域劃分和不同區(qū)域間的數(shù)據(jù)交互較為困難，而且一般還需要設(shè)置緩沖區(qū)。Torre等［24］發(fā)現(xiàn)，該類耦合算法的計算結(jié)果對區(qū)域劃分和緩沖區(qū)的設(shè)置非常敏感。近年來，一些改進的粒子類方法也被相繼提出用于克服原始DSMC 方法在低速流動統(tǒng)計噪聲過大的問題，以及在連續(xù)和近連續(xù)流區(qū)域網(wǎng)格尺寸和時間步長受限于分子平均自由程和分子平均碰撞時間的不足。例如，Sun 和Boyd［25］提出了DSMC 信息保存（Direct simulation Monte Carlo-information preservation，DSMC-IP）方法；Fei 等［26］提出了多尺度的統(tǒng)一粒子BGK（Unified stochastic particle-BGK，USP-BGK）方法等。

不同于粒子類方法，DVM 直接更新離散分布函數(shù)，有效避開了統(tǒng)計噪聲問題，并且可以方便地引入傳統(tǒng)計算流體力學(xué)的隱式加速算法來提高計算效率。在DVM 框架下，通過引入氣體分子碰撞松弛參數(shù)和當(dāng)?shù)仄胶鈶B(tài)速度分布函數(shù)來建立Boltzmann 方程碰撞積分統(tǒng)一模型，同時利用高效算子分裂方法和大規(guī)模并行技術(shù)來求解該模型方程，Li 和Zhang［27］提 出 了 氣 體 動 理 論 統(tǒng) 一 算 法（Gas kinetic unified algorithm，GKUA）；通過利用Boltzmann-BGK 方程的局部積分解來同時計算介觀 方 程 和 宏 觀 伴 隨 方 程 的 通 量，Xu 和Huang［28］提出了統(tǒng)一氣體動理學(xué)格式（Unified gas kinetic scheme，UGKS）；通過利用Boltzmann-BGK 方程的局部特征解來計算介觀方程的通量，Guo 等［29］提出了離散統(tǒng)一氣體動理學(xué)格式（Discrete unified gas kinetic scheme，DUGKS）。隨后，Chen 等［30］通過對介觀方程和宏觀伴隨方程的通量進行簡化，提出了簡化版本的UGKS；通過引入LU-SGS 技術(shù)和多重網(wǎng)格加速收斂技術(shù)，Zhu 等［31］改進了UGKS 的計算效率。此外，在采用傳統(tǒng)DVM 求解介觀方程的基礎(chǔ)上，通過引入宏觀伴隨方程并利用含碰撞項Boltzmann-BGK 方程的局部解來計算該方程通量，Yang 等［32］提 出 了 改 進 離 散 速 度 方 法（Improved discrete velocity method，IDVM）；Su 等［33］通過在計算宏觀伴隨方程通量時引入高階修正項，提出了合成迭代加速算法（General synthetic iterative scheme，GSIS）。由于在通量計算過程中同時考慮了分子遷移和碰撞的影響，基于DVM 框架發(fā)展的算法的網(wǎng)格尺度和時間步長不再受限于分子平均自由程和平均碰撞時間，因而有效地克服了DSMC 方法在連續(xù)和近連續(xù)流區(qū)域的計算困難，實現(xiàn)了全流域的統(tǒng)一求解。

但由于需要在位置空間和速度空間同時離散，DVM 所需的存儲量和計算量非常大，因而早期發(fā)展較為緩慢。近年來，隨著計算機內(nèi)存和運算速度的提升，基于DVM 框架發(fā)展的算法已取得了喜人的成績并逐漸被應(yīng)用于航天科技、微電子機械系統(tǒng)和真空技術(shù)等多個領(lǐng)域。本文首先對該類方法的研究進展進行回顧，尤其是GKUA 和UGKS 兩種算法的構(gòu)造途徑，并介紹他們在跨流域問題中的應(yīng)用。隨后，本文將介紹作者團隊近年來發(fā)展的IDVM 及其研究進展，并通過引入雙時間步格式構(gòu)造非定常隱式IDVM，用于非定?？缌饔騿栴}求解。最后，通過剖析基于DVM 框架發(fā)展的數(shù)值方法面臨的挑戰(zhàn)，展望未來跨流域問題的一些研究方向。

1 Boltzmann模型方程和DVM框架

1.1 Boltzmann 模型方程和速度空間離散

通過引入氣體分子速度分布函數(shù)f，則氣體動力學(xué)中感興趣的各種宏觀量便可以通過對f求矩而得到。由此可見，除了采用常見的宏觀守恒律方程，流體系統(tǒng)亦可由分布函數(shù)f的演化方程來描述。1872 年，Boltzmann 給出了分布函數(shù)對位置空間、速度空間和時間的變化率的關(guān)系，即Boltzmann 方程

式中：c=ξ-u為分子的最可幾熱運動速度；u為宏觀速度；ρ為密度；T為溫度；Rg為氣體常數(shù)。除特殊說明外，本文中約定用黑體來表示矢量，用對應(yīng)的白斜體來表示矢量的長度。

依據(jù)速度分布函數(shù)f的定義以及分子碰撞過程中的質(zhì)量、動量和能量守恒關(guān)系（即相容性條件），氣體動力學(xué)中常見的宏觀量可以表示為

由于方程（1）的碰撞項過于復(fù)雜且對于流體力學(xué)計算難以形成明確的數(shù)學(xué)表達式，在速度空間直接離散求解該方程對于實際問題的表征難有意義［34］。因此，各類簡化的碰撞模型被提出用于近似該碰撞項，例如BGK 模型［35］，ES-BGK 模型［36］和Shakhov-BGK 模型［37］等。BGK 模型由Bhatnagar，Gross 和Krook 提出，它可以滿足質(zhì)量、動量和能量守恒，滿足熵增條件，并且導(dǎo)出的平衡態(tài)即為Maxwell 分布。但是，該模型對應(yīng)的普朗特數(shù)為1，與通過原始方程（1）推導(dǎo)得到的正確值2/3 不一致，因此不能同時保證正確的黏性系數(shù)與熱導(dǎo)率。ES-BGK模型由Holway 提出，通過在平衡態(tài)分布函數(shù)中引入應(yīng)力修正項來實現(xiàn)對BGK 模型普朗特數(shù)的修正；Shakhov-BGK 模型則通過在平衡態(tài)分布函數(shù)中引入熱流修正項來實現(xiàn)對BGK 模型普朗特數(shù)的修正。由于可以恢復(fù)正確的普朗特數(shù)，ES-BGK 模型和Shakhov-BGK 模型都能同時保證正確的黏性系數(shù)與熱導(dǎo)率。但大多數(shù)情況下Shakhov-BGK 模型的精度優(yōu)于ES-BGK 模型，因此Shakhov-BGK 模型獲得了更廣泛的應(yīng)用［38］。采用Shakhov-BGK 模型作為碰撞項的Boltzmann 方程可以改寫為

式中Pr為普朗特數(shù)。相較于方程（1），采用方程（7）可使得求解過程大為簡化。因此，目前基于DVM 框架發(fā)展的算法大都是直接求解方程（7）。

雖然Boltzmann 模型方程的碰撞項已大為簡化，但分布函數(shù)f仍然是時間、位置空間和速度空間的七維變量，需要離散化方能數(shù)值求解。DVM首先在速度空間對Boltzmann 模型方程進行離散，獲得離散速度形式的Boltzmann 模型方程

由于速度空間的網(wǎng)格量直接影響DVM 的計算量和存儲量，為了盡可能優(yōu)化速度空間的網(wǎng)格分布，Gutnic 等［43］、Kolobov 和Arslanbekov［44］、Chen等［45］發(fā)展了速度空間自適應(yīng)DVM。該方法的內(nèi)存需求相較于采用均勻網(wǎng)格的Newton-Cotes 求積可以降低1～2 個數(shù)量級，但由于采用自適應(yīng)算法之后不同位置空間網(wǎng)格的分布函數(shù)對應(yīng)的離散速度不一致，需要頻繁的插值運算，程序處理較為復(fù)雜。最近，Zhao 等［46］采用降階模型對計算結(jié)果進行模態(tài)分析，通過提取主導(dǎo)模態(tài)對應(yīng)的離散速度點，并僅更新這些離散點的分布函數(shù)，有效減少了DVM 的計算量。但由于不同來流參數(shù)對應(yīng)的主導(dǎo)模態(tài)不一致，該方法僅適用于來流參數(shù)變化較小的流動問題求解。為了應(yīng)對工程實際問題中通常需要計算一系列不同來流參數(shù)工況的情形，Yang等［47］進一步提出了基于參數(shù)化的降階離散速度方法。首先，從全部工況中挑選少數(shù)工況作為預(yù)計算工況，并利用IDVM 計算其結(jié)果。其次，基于這些預(yù)計算工況得到的離散分布函數(shù)，采用奇異值分解求出主導(dǎo)模態(tài)對應(yīng)的降階子空間，并利用對數(shù)和指數(shù)映射，在Grassmann 流型及其切空間上插值計算其余工況所對應(yīng)的降階子空間。最后，通過離散經(jīng)驗插值算法搜尋降階子空間對應(yīng)的離散速度點，從而構(gòu)成縮減的離散速度空間?；诖?，其余工況便可以直接在對應(yīng)的縮減離散速度空間上進行求解，從而有效地減少了計算量。

此外，由于DVM 中采用數(shù)值求積來計算宏觀量，其結(jié)果必然會與直接積分存在誤差，因而導(dǎo)致相容性條件不能精確滿足，影響結(jié)果精度和數(shù)值穩(wěn)定 性。為 了 解 決 這 一 問 題，Titarev［48］，江 定 武等［49］，Yang 等［50］發(fā)展了守恒型DVM，通過引入內(nèi)迭代強制滿足相容性條件。采用該類方法可以在較稀的速度空間網(wǎng)格上得到網(wǎng)格無關(guān)解，減少的計算量最大可達2/3。

1.2 GKUA 研究進展

針對航天飛行器再入各流域多尺度非平衡繞流特點，為了反映再入過程不同流域氣體分子相互作用、稀薄程度、分子模型與內(nèi)部能量變化關(guān)系，通過引入氣體分子碰撞松弛參數(shù)ν和當(dāng)?shù)仄胶鈶B(tài)速度分布函數(shù)fS來模型化表征Boltzmann 方程碰撞積分項對航天飛行器再入跨流域氣動力/熱繞流特性，可確立描述各流域全局馬赫數(shù)復(fù)雜流動輸運現(xiàn)象統(tǒng)一的Boltzmann 模型速度分布函數(shù)方程，其無量綱形式為

由此，可將最新的計算流體力學(xué)中的有限差分格式等引入到基于離散速度坐標(biāo)形式的Boltzmann 方程中，在位置空間和時間方向?qū)υ摲匠踢M行求解?；诖?，Li 等［51-52］提出了GKUA用于模擬高超聲速跨流域流動問題。其中，采用算子分裂思想結(jié)合二階顯式Runge-Kutta 方法和無波動、無自由參數(shù)、具有耗散性（Nonoscillatory nonfree dissipative，NND）差分格式［53］，方程（10）可以在位置空間和時間方向進一步離散為

式 中：LS、Lξ、Lη、Lζ分 別 為 如 下 方 程 的 二 階 差 分算子

完成分布函數(shù)更新后，便可利用方程（3）和方程（6）計算新時刻的守恒量和熱通量，并利用方程（2）和（8）計算新時刻的當(dāng)?shù)仄胶鈶B(tài)分布，從而開始下一時刻推進求解。該方法中，時間步長Δt僅由差分格式穩(wěn)定性條件決定，與當(dāng)?shù)貧怏w分子平均碰撞時間無關(guān)［3，20，34，52］。另外，由于其宏觀流動量是通過離散速度分布函數(shù)對所有離散速度坐標(biāo)點進行離散積分歸約求和來更新，與常規(guī)計算流體力學(xué)位置空間具體網(wǎng)格尺度無關(guān)，因而位置空間流場網(wǎng)格劃分不受任何限制，可在大尺度網(wǎng)格快速計算收斂，確保了GKUA 計算復(fù)雜飛行器高超聲速氣動力/熱繞流特性，解決實際應(yīng)用問題的準(zhǔn)確可行性［54-59］。

除了有限差分格式，吳俊林等［55］還發(fā)展了基于有限體積格式的GKUA。另外，為了提高計算效 率 和 穩(wěn) 定 性，Peng 等［56］發(fā) 展 了 隱 式 版 本 的GKUA。由于僅求解介觀Boltzmann 型速度分布函數(shù)方程并且在同一時間步內(nèi)每個離散速度坐標(biāo)點處的分布函數(shù)的求解是相互獨立的，GKUA 的計算過程較為簡單，且非常易于在離散速度空間分塊并行求解。基于HPF、MPI+OpenMP、MPI+OpenACC 程序構(gòu)架，各種并行版本的GKUA 被提出并應(yīng)用于跨流域工程實際問題求解。此外，考慮到高溫多原子氣體必然會存在轉(zhuǎn)動和振動自由度，為了準(zhǔn)確模擬高溫高馬赫數(shù)下多原子氣體內(nèi)能激發(fā)對跨流域非平衡流動的影響，李志輝等［57］發(fā)展了同時考慮平動-轉(zhuǎn)動能松弛特性的GKUA，彭傲平等［58］發(fā)展了同時考慮平動-轉(zhuǎn)動-振動能松弛特性的GKUA。 應(yīng)用發(fā)展的算法，李志輝等［3，6，27，34，52，57］構(gòu)建了適用于返回艙從外層空間自由分子流再入到近地面連續(xù)流的跨流域空氣動力學(xué)一體化模擬平臺。圖2 展示了采用該平臺計算得到的返回艙再入120～70 km 流場壓力分布。最近，結(jié)合空間站建設(shè)與運維，為了研究服役期滿大型航天器如空間實驗室、貨運飛船等再入氣動力/熱環(huán)境致結(jié)構(gòu)變形/失效解體的非線性力學(xué)響應(yīng)行為，李志輝等［54，59］將瞬態(tài)熱傳導(dǎo)方程與材料熱彈性動力學(xué)方程引入GKUA 中，構(gòu)建了跨流域動態(tài)力/熱耦合計算平臺，實現(xiàn)了再入強氣動力/熱環(huán)境致材料變形/失效解體的統(tǒng)一求解，如圖3 所示。研究表明，GKUA 作為一種基于氣體分子速度分布函數(shù)演化更新實時捕捉宏觀氣體流動特征的介觀模擬方法，可以較好地求解全流域氣體流動問題，在航空航天領(lǐng)域得到了成功運用。目前算法考慮了多種非平衡效應(yīng)，包括平動-轉(zhuǎn)動非平衡、平動-轉(zhuǎn)動-振動非平衡，建立了可靠模擬大型復(fù)雜結(jié)構(gòu)航天器（如我國天宮一號目標(biāo)飛行器）服役期滿離軌隕落再入解體跨流域高超聲速氣動力/熱非平衡繞流問題的氣體動理論統(tǒng)一算法大規(guī)模并行計算應(yīng)用研究平臺。

圖2 返回艙以第一宇宙速度7.5 km/s 再入跨流域（120～70 km）氣動模擬軟件系統(tǒng)實時計算與姿態(tài)配平繞流壓力分布[3]Fig.2 Real-time computing pressure distribution during the spacecraft capsule re-entry (120—70 km) with the first cosmic velocity of 7.5 km/s[3]

圖3 類天宮飛行器在H = 120～100 km、v∞= 7.6 km/s 條件下強氣動力/熱環(huán)境致結(jié)構(gòu)溫度變化和穩(wěn)態(tài)變形[54]Fig.3 Temperature distributions and steady-state deformation of structure at H = 120—100 km and v∞=7.6 km/s around Tiangong type spacecraft by strong aerodynamic heating[54]

1.3 UGKS 研究進展

與GKUA 不同，UGKS 同步求解了介觀的Boltzmann 模型方程和對應(yīng)的宏觀伴隨方程。宏觀伴隨方程實際上對應(yīng)于Boltzmann 模型方程在速度空間求守恒矩

式中：xij為單元界面位置；f(xij，ξα，t)為單元界面上的離散分布函數(shù)；f(xij-ξαt，ξα，0)為單元界面周圍的初始離散分布函數(shù)；fS(xij-ξα(tt′)，ξα，t′)為t′時 刻 單 元 界 面 周 圍 的 平 衡 態(tài) 分 布 函數(shù)。實際計算中，f(xij-ξαt，ξα，0)可由n時刻單元中心的離散分布函數(shù)插值計算得到，而平衡態(tài)分布fS(xij-ξα(t-t′)，ξα，t′)則 可 由 宏 觀 物 理 量 及其導(dǎo)數(shù)來計算。最終，f(xij，ξα，t)可寫為

由于f(xij，ξα，t)是時間的函數(shù)，通量計算時需要取其在(0，Δt)積分的時間平均

由方程（18）可知，當(dāng)分子遷移時間遠小于分子平均碰撞時間時（t?τ），單元界面分布函數(shù)中初始分布f(xij-ξαt，ξα，0)占主導(dǎo)，表現(xiàn)為自由分子流情形的完全自由遷移現(xiàn)象；而當(dāng)分子遷移時間遠大于分子平均碰撞時間時（t?τ），單元界面周圍的平衡態(tài)分布fS(xij-ξα(t-t′)，ξα，t′)占主導(dǎo)，表現(xiàn)為連續(xù)流情形的分子頻繁碰撞現(xiàn)象。由于實際計算中分子的遷移時間即為網(wǎng)格尺度對應(yīng)的時間步長Δt，方程（18）等效于在網(wǎng)格尺度上的流動建模，將主導(dǎo)稀薄流的分子自由遷移機制和主導(dǎo)連續(xù)流的分子碰撞機制的權(quán)重與網(wǎng)格尺度有機聯(lián)系起來。由于通量計算時分子的遷移和碰撞過程相互耦合，UGKS 的推進時間步長不受限于當(dāng)?shù)胤肿拥钠骄鲎矔r間，且其網(wǎng)格尺度也不受限于當(dāng)?shù)胤肿悠骄杂沙蹋瑥亩沟迷摲椒ㄔ谌饔蚓色@得準(zhǔn)確高效的計算結(jié)果。

基于上述優(yōu)勢，UGKS 已成功應(yīng)用于從連續(xù)到稀薄流的許多流動計算中。例如，Liu 等［61］發(fā)展了同時考慮平動-轉(zhuǎn)動能松弛特性的UGKS；Wang等［62］發(fā)展了同時考慮平動-轉(zhuǎn)動-振動能松弛特性的UGKS；Sun 等［63］將UGKS 應(yīng) 用 于 輻 射 傳 熱 問題；Xiao 等［64］將UGKS 應(yīng)用于多組分流問題；Liu和Xu［65］將UGKS 應(yīng)用于等離子體問題；Tan 等［66］將UGKS 應(yīng) 用于中子輸 運問題；Liu［67］將UGKS 應(yīng)用于顆粒流問題等。同時，為了進一步提高UGKS的 性 能，Chen 等［68］和Yang 等［69］發(fā) 展 了 內(nèi) 存 縮 減UGKS 用于定?？缌饔蛄鲃訂栴}的求解；Zhu等［70-71］發(fā)展了隱式版本的UGKS 用于定常和非定常流動問題求解，使得計算效率提高了1～2 個數(shù)量 級。應(yīng)用UGKS，Li 等［72］模擬了類X38 高超聲速飛行器在不同克努森數(shù)時的氣動力/熱問題并與DSMC 結(jié) 果 進 校 了 比 較，如 圖4 所 示。Chen 等［45］結(jié)合動網(wǎng)格技術(shù)，應(yīng)用UGKS 模擬了稀薄環(huán)境下噴流推進問題，如圖5 所示。該問題中，噴管內(nèi)部為連續(xù)流動，而噴管外部為稀薄流動，因此需要保證算法在連續(xù)到稀薄整個流域范圍都能獲得可靠的計算結(jié)果。

圖4 采用UGKS 和DSMC 計算的類X38 飛行器溫度場比較(Argon, α=20°)[72]Fig.4 Comparison of temperature contours of X38-like model[72]

圖5 稀薄環(huán)境下噴管噴流計算結(jié)果[45]Fig.5 Numerical results of gas injection to the rarefied atmosphere[45]

1.4 IDVM 研究進展

GKUA 僅求解Boltzmann 型速度分布函數(shù)方程，而且同一時間步內(nèi)每個離散分布函數(shù)的演化是相互獨立的。因此，該算法實施較為簡單，且非常易于在離散速度空間分塊并行求解。由于其與分子層級的遷移和碰撞解耦無關(guān)，GKUA 的時間步長僅由所使用的差分格式的穩(wěn)定性條件決定。不僅時間步長與分子平均碰撞時間無關(guān)，而且位置空間網(wǎng)格劃分尺度也與氣體分子平均自由程無關(guān)，確保了GKUA 在大網(wǎng)格尺度計算復(fù)雜高超聲速飛行器氣動力/熱非平衡繞流問題的高可靠性與工程置信度。然而，由于GKUA 的宏觀量由當(dāng)?shù)馗麟x散速度坐標(biāo)點的分布函數(shù)實時歸約求和進行更新，Boltzmann 型速度分布函數(shù)方程的碰撞項和對流項同步顯式或隱式離散差分計算，對連續(xù)流問題的求解按所用差分格式穩(wěn)定性條件確定的時間步長會使得收斂變慢。相比較而言，UGKS 同步求解了Boltzmann 模型方程和相應(yīng)的宏觀伴隨方程，并將分子的遷移和碰撞過程耦合處理。同時，由于同步求解了宏觀伴隨方程，其結(jié)果可用于預(yù)估新時刻的平衡態(tài)，從而實現(xiàn)了Boltzmann 模型方程中碰撞項的全隱式離散，保證了連續(xù)流區(qū)域的計算效率。但是，由于通量計算時所采用的Boltzmann-BGK 方程的局部積分解較為復(fù)雜，稀薄流區(qū)域的計算效率受到了一定的影響，而且為了獲得單元界面的積分解，需要聯(lián)立界面周圍所有的離散分布函數(shù)來計算界面的平衡態(tài)分布。

在上述兩種算法的啟發(fā)下，為了保證從稀薄到連續(xù)整個流域范圍的計算精度和計算效率，同時使算法的實施較為簡單，Yang 等［73］發(fā)展了IDVM。與UGKS 類似，該算法也同步求解了Boltzmann 模型方程和相應(yīng)的宏觀伴隨方程，以便實現(xiàn)Boltzmann 模型方程中碰撞項的全隱式離散。但不同于UGKS，IDVM 在求解Boltzmann 模型方程時仍然將分子的遷移和碰撞過程解耦處理，以便保持同一時間步內(nèi)每個離散分布函數(shù)的演化相互獨立的優(yōu)勢。此外，為了提高連續(xù)流區(qū)域的計算精度和計算效率，IDVM 在求解宏觀伴隨方程時同時考慮了分子的遷移和碰撞過程對通量計算的影響。具體地，IDVM 將Boltzmann 模型方程和相應(yīng)的宏觀伴隨方程離散為如下形式

式中：fDVM為不含碰撞項Boltzmann 方程在單元界面的局部解，即方程（25）；fNS為對應(yīng)于連續(xù)流極限的分布函數(shù)。將方程（26）代入方程（4），即可得到宏觀伴隨方程的通量

由方程（25）和方程（26）可知，IDVM 在計算Boltzmann 模型方程通量和宏觀伴隨方程通量時所采用的單元界面分布函數(shù)是不一樣的。原因在于，采用不考慮分子碰撞影響的分布函數(shù)來計算Boltzmann 模型方程的通量雖然會影響該方程在連續(xù)流區(qū)域的計算精度和計算效率，但卻可以保證同一時間步內(nèi)每個離散分布函數(shù)的演化是相互獨立的，以保持算法的簡單性。實際上，在連續(xù)流區(qū)域由宏觀伴隨方程主導(dǎo)流場的解，而且此時宏觀伴隨方程的通量即為Navier-Stokes 方程的通量，因此只要引入合理的宏觀伴隨方程即可獲得該區(qū)域的準(zhǔn)確計算結(jié)果。綜上所述，該策略既保持了原始DVM 的簡單性優(yōu)勢，也提高了連續(xù)流區(qū)域的計算精度和計算效率。但是，由于Boltzmann 模型方程的通量和宏觀伴隨方程的通量計算不一致，該方法在理論上還存在不自洽的問題，其具體影響和改進措施還有待進一步研究。最近，通過引入LU-SGS迭代，Yang 等［74］發(fā)展了三維隱式IDVM 用于定?？缌饔蛄鲃訂栴}的求解。如圖6 所示，在連續(xù)流區(qū)域，IDVM 相較于傳統(tǒng)DVM 計算精度更高，且計算效率提高了1～2 個數(shù)量級。另外，通過引入內(nèi)迭代技術(shù)，Yang 等［75］實現(xiàn)了定常IDVM 計算效率的進一步提升，并模擬了高超聲速稀薄流狀態(tài)下的獵戶座飛船返回艙的氣動力/熱問題，如圖7 所示。

圖6 三維頂蓋驅(qū)動流問題的計算結(jié)果[74]Fig.6 Numerical results for 3D lid-driven cavity flow[74]

圖7 獵戶座飛船返回艙高超聲速稀薄流的密度和壓力分布[75]Fig.7 Density and pressure contours for hypersonic rarefied flow around an Orion crew module[75]

2 非定常隱式IDVM 及其應(yīng)用

2.1 非定常隱式IDVM

為了實現(xiàn)非定?？缌饔蛄鲃訂栴}的準(zhǔn)確高效求解，本文將先前發(fā)展的隱式定常IDVM 進一步拓展到了非定常情形?；陔p時間步LU-SGS 迭代的思想，在Boltzmann 模型方程和宏觀伴隨方程中同時增加了偽時間導(dǎo)數(shù)項

2.2 算例驗證

為了驗證上述非定常隱式IDVM 的性能，本文模擬了墻壁束縛瑞利流，并將計算結(jié)果與非定常隱式UGKS 進行比較。如圖8 所示，初始時刻兩個間隔1 m 長4 m 的平行平板之間充滿靜止的氬氣。平板的溫度和氬氣的溫度均為273 K，氬氣的分子質(zhì)量為6.63×10-26kg，對應(yīng)的氣體常數(shù)為Rg=208.13 J/（kg·K）。啟動瞬間，左右兩側(cè)的平板以10 m/s 的速度向上運動，平板溫度為373 K。由于幾何對稱性，實際計算中僅考慮左半側(cè)區(qū)域流場。

圖8 墻壁束縛瑞利流示意圖[71]Fig.8 Schematic for wall-bounded Rayleigh flow[71]

為了與文獻條件一致，本文將該算例的克努森數(shù)取為Kn= 0.05。分子速度空間采用包含28×28 個離散點的Gauss-Hermite 積分，位置空間采用12 800 個均勻四邊形非結(jié)構(gòu)網(wǎng)格進行離散。無量綱的外迭代時間步長取為Δt=0.01。圖9 展示了t= 1.5 時刻水平和垂直中心線溫度、x方向速度分量和y方向速度分量分布，當(dāng)前計算結(jié)果與UGKS結(jié)果較好吻合。另外，圖10 展示了t= 1.5 時刻上下平板的壓力、熱通量和切應(yīng)力分布，當(dāng)前結(jié)果也與GUKS 結(jié)果吻合良好。該算例驗證了提出的方法在非定?？缌饔蛄鲃拥挠嬎隳芰?。

圖9 墻壁束縛瑞利流在Kn = 0.05 和t= 1.5 時水平中心線和垂直中心線流動變量分布圖Fig.9 Flow variables along the horizontal and vertical central lines for wall-bounded Rayleigh flow at Kn = 0.05 and t= 1.5

圖10 墻壁束縛瑞利流在Kn = 0.05 和t= 1.5 時上下表面流動變量分布圖Fig.10 Flow variables along the upper and bottom surfaces for wall-bounded Rayleigh flow at Kn = 0.05 and t= 1.5

3 總結(jié)與展望

本文首先回顧了基于速度空間離散的幾類Boltzmann 模型方程求解算法，展示了它們在低速到高超聲速、連續(xù)到稀薄全流域范圍的應(yīng)用。同時，本文還將筆者前期發(fā)展的定常隱式IDVM 擴展到了非定常情形，并通過非定?？缌饔騿栴}進行了驗證。這類方法由于直接演化分布函數(shù)，可以有效避免粒子類方法的統(tǒng)計噪聲問題，同時可以采用傳統(tǒng)計算流體力學(xué)中的加速技術(shù)來提高計算效率。在微電子機械和真空技術(shù)等低速跨流域問題以及不含離解和電離化學(xué)非平衡效應(yīng)的高速跨流域問題中，該類算法已經(jīng)展現(xiàn)出一定的優(yōu)勢。但是對于航天領(lǐng)域涉及的高溫高速和真實氣體效應(yīng)的非平衡流動，該類算法也還存在許多丞待解決的問題，需要持續(xù)不斷地改進完善。

（1）由于基于DVM 框架發(fā)展的算法需要在位置空間和速度空間同時進行離散，因而實際工程問題的計算量和內(nèi)存開銷非常大。尤其是對于高超聲速跨流域問題，由于流動速度和溫度非常高，速度空間截斷區(qū)域也需要相應(yīng)擴大，方能保證求矩計算宏觀量的精度。最近，Liu 等［76］和Xu［77］提出了統(tǒng)一氣體動理學(xué)波粒子（Unified gas kinetic wave particle，UGKWP）方法，采用粒子來表征稀薄效應(yīng)部分的貢獻，而連續(xù)效應(yīng)部分的貢獻直接采用解析的方式來計算。在連續(xù)流區(qū)域，UGKWP 方法自動退化為Navier-Stokes 方程求解器，而在稀薄流區(qū)域，該方法也具有粒子在速度空間自適應(yīng)分布的優(yōu)勢。但是，由于引入了粒子的計算，UGKWP 方法同樣也存在統(tǒng)計噪聲問題，以及不便應(yīng)用于非定常流動求解和不易于采用傳統(tǒng)計算流體力學(xué)的加速技術(shù)來提高計算效率的問題。

（2）在Boltzmann 模型方程中，每增加一種能量松弛方式和化學(xué)組分，都需要增加相應(yīng)的分布函數(shù)來對其進行描述。對于高超聲速稀薄流動問題，由于溫度通常高達10 000 K，不僅需要考慮分子的轉(zhuǎn)動松弛和振動松弛，還需要考慮離解和電離化學(xué)非平衡效應(yīng)，這給基于DVM 框架發(fā)展的算法帶來了極大的困難和挑戰(zhàn)。一方面是如何建立合理的模型方程，另一方面是如何設(shè)計針對模型方程的高效數(shù)值算法。文獻［78］從量子力學(xué)的角度出發(fā)建立了WCU 方程，對每個能級的分布函數(shù)寫出了一個Boltzmann 方程，其遷移項保持不變，而碰撞項是單組分Boltzmann 方程碰撞項的唯象擴展。該模型方程有望實現(xiàn)對真實氣體效應(yīng)的模擬，但由于求解困難，目前文獻中尚未發(fā)現(xiàn)可工程實用的相應(yīng)算法。