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

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

飛船返回艙、HTV-2、東風(fēng)-17 等往返大氣層與空間軌道跨流域飛行器（圖1）的研究，一直以來都備受各軍事和航天強(qiáng)國的關(guān)注，在載人登月、深空探測(cè)和國防領(lǐng)域有著重要的應(yīng)用價(jià)值。然而，由于這類飛行器的飛行高度范圍非常大，經(jīng)常會(huì)跨越稀薄流和連續(xù)流整個(gè)流域范圍，對(duì)試驗(yàn)研究和數(shù)值模擬都提出了嚴(yán)苛的考驗(yàn)［1-5］。一方面，由于這類飛行器所需要考慮的相似參數(shù)眾多，流動(dòng)相似準(zhǔn)則要求苛刻，地面風(fēng)洞實(shí)驗(yàn)設(shè)備對(duì)于同時(shí)復(fù)現(xiàn)其所處的低雷諾數(shù)與高焓非平衡流動(dòng)特征通常無能為力，而且自由飛試驗(yàn)的代價(jià)極大，數(shù)據(jù)采集也非常困難。另一方面，由于飛行軌跡跨越不同的流域，對(duì)飛行器的流動(dòng)模擬需要綜合考慮全流域范圍的計(jì)算精度和計(jì)算效率。此外，即便在某一環(huán)境克努森數(shù)下，飛行器周圍流場(chǎng)的局部克努森數(shù)也會(huì)跨越多個(gè)數(shù)量級(jí)，同時(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í)驗(yàn)研究，數(shù)值模擬無疑是解決該類問題的一種經(jīng)濟(jì)有效的手段。但由于稀薄效應(yīng)的存在，傳統(tǒng)的基于連續(xù)性假設(shè)的Navier-Stokes 方程不再適用于該類問題［9］。為了有效地研究該類問題，通常需要借助于從分子動(dòng)理學(xué)理論（分子運(yùn)動(dòng)論）發(fā)展而來的Boltzmann 方程。該方程是位置空間、速度空間和時(shí)間的七維多相空間幾率密度分布函數(shù)方程［10］，同時(shí)其碰撞積分項(xiàng)是一個(gè)結(jié)構(gòu)復(fù)雜未有明確表達(dá)式的五重積分，因而直接求解極為困難?；谠摲匠蹋ㄟ^合理簡化，目前發(fā)展了兩類典型的數(shù)值求解算法。第一類是粒子方法，它通過模擬假想粒子的遷移和碰撞過程來實(shí)現(xiàn)對(duì)流場(chǎng)的仿真，建模過程等價(jià)于對(duì)Boltzmann 方程的輸運(yùn)項(xiàng)和碰撞項(xiàng)進(jìn)行解耦處理。典型代表為直接模擬蒙特卡羅（Direct simulation Monte Carlo，DSMC）方法［4，8，11-13］。第二類是離散速度或離散坐標(biāo)法（Discrete velocity or discrete ordinate method，DVM 或DOM［14-18］），其在位置空間和速度空間同時(shí)離散求解Boltzmann 方程，并且為了避免對(duì)碰撞積分不定式的計(jì)算，通常采用簡化模型來近似碰撞積分項(xiàng)。

DSMC 方法是現(xiàn)階段解決高超聲速飛行器在稀薄流域最有效實(shí)用的粒子類方法，由Bird［19］首次提出。該方法的關(guān)鍵是在小于當(dāng)?shù)胤肿悠骄鲎矔r(shí)間步長Δt內(nèi)將分子運(yùn)動(dòng)和碰撞解耦［3-4，8-9］，即讓所有分子運(yùn)動(dòng)一定距離并考慮邊界反射，然后計(jì)算此時(shí)間段內(nèi)具有代表性的分子間碰撞。在碰撞取樣數(shù)趨于無窮大時(shí)，DSMC 方法模擬得到的速度分布函數(shù)收斂到Boltzmann 方程的修正形式［20］。所以，為了確保模擬不失真，通常要求DSMC 方法的網(wǎng)格尺寸小于分子平均自由程，且時(shí)間步長小于分子平均碰撞時(shí)間。在稀薄流區(qū)域，由于分子數(shù)較少，分子平均自由程和平均碰撞時(shí)間較大，DSMC方法具有較高的計(jì)算效率和計(jì)算精度。但在連續(xù)和近連續(xù)流區(qū)域，由于分子平均自由程和平均碰撞時(shí)間較小，DSMC 方法的計(jì)算效率受到了極大的限制。此外DSMC 方法還存在統(tǒng)計(jì)噪聲問題，以及不便應(yīng)用于非定常流動(dòng)模擬和不易于采用隱式算法加速等困難。為了應(yīng)用于跨流域問題的求解，一些學(xué) 者 提 出 了Navier-Stokes/DSMC 耦 合 算 法［21-23］，其思想是將計(jì)算域劃分為不同的子區(qū)域分別應(yīng)用DSMC 方法和Navier-Stokes 方程求解器進(jìn)行求解。雖然該類耦合算法能在各自的區(qū)域獲得準(zhǔn)確高效的計(jì)算結(jié)果，但區(qū)域劃分和不同區(qū)域間的數(shù)據(jù)交互較為困難，而且一般還需要設(shè)置緩沖區(qū)。Torre等［24］發(fā)現(xiàn)，該類耦合算法的計(jì)算結(jié)果對(duì)區(qū)域劃分和緩沖區(qū)的設(shè)置非常敏感。近年來，一些改進(jìn)的粒子類方法也被相繼提出用于克服原始DSMC 方法在低速流動(dòng)統(tǒng)計(jì)噪聲過大的問題，以及在連續(xù)和近連續(xù)流區(qū)域網(wǎng)格尺寸和時(shí)間步長受限于分子平均自由程和分子平均碰撞時(shí)間的不足。例如，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)計(jì)噪聲問題，并且可以方便地引入傳統(tǒng)計(jì)算流體力學(xué)的隱式加速算法來提高計(jì)算效率。在DVM 框架下，通過引入氣體分子碰撞松弛參數(shù)和當(dāng)?shù)仄胶鈶B(tài)速度分布函數(shù)來建立Boltzmann 方程碰撞積分統(tǒng)一模型，同時(shí)利用高效算子分裂方法和大規(guī)模并行技術(shù)來求解該模型方程，Li 和Zhang［27］提 出 了 氣 體 動(dòng) 理 論 統(tǒng) 一 算 法（Gas kinetic unified algorithm，GKUA）；通過利用Boltzmann-BGK 方程的局部積分解來同時(shí)計(jì)算介觀 方 程 和 宏 觀 伴 隨 方 程 的 通 量，Xu 和Huang［28］提出了統(tǒng)一氣體動(dòng)理學(xué)格式（Unified gas kinetic scheme，UGKS）；通過利用Boltzmann-BGK 方程的局部特征解來計(jì)算介觀方程的通量，Guo 等［29］提出了離散統(tǒng)一氣體動(dòng)理學(xué)格式（Discrete unified gas kinetic scheme，DUGKS）。隨后，Chen 等［30］通過對(duì)介觀方程和宏觀伴隨方程的通量進(jìn)行簡化，提出了簡化版本的UGKS；通過引入LU-SGS 技術(shù)和多重網(wǎng)格加速收斂技術(shù)，Zhu 等［31］改進(jìn)了UGKS 的計(jì)算效率。此外，在采用傳統(tǒng)DVM 求解介觀方程的基礎(chǔ)上，通過引入宏觀伴隨方程并利用含碰撞項(xiàng)Boltzmann-BGK 方程的局部解來計(jì)算該方程通量，Yang 等［32］提 出 了 改 進(jìn) 離 散 速 度 方 法（Improved discrete velocity method，IDVM）；Su 等［33］通過在計(jì)算宏觀伴隨方程通量時(shí)引入高階修正項(xiàng)，提出了合成迭代加速算法（General synthetic iterative scheme，GSIS）。由于在通量計(jì)算過程中同時(shí)考慮了分子遷移和碰撞的影響，基于DVM 框架發(fā)展的算法的網(wǎng)格尺度和時(shí)間步長不再受限于分子平均自由程和平均碰撞時(shí)間，因而有效地克服了DSMC 方法在連續(xù)和近連續(xù)流區(qū)域的計(jì)算困難，實(shí)現(xiàn)了全流域的統(tǒng)一求解。

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

1 Boltzmann模型方程和DVM框架

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

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

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

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

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

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

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

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

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

1.2 GKUA 研究進(jìn)展

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

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

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

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

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

圖2 返回艙以第一宇宙速度7.5 km/s 再入跨流域（120～70 km）氣動(dòng)模擬軟件系統(tǒng)實(shí)時(shí)計(jì)算與姿態(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 條件下強(qiáng)氣動(dòng)力/熱環(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 研究進(jìn)展

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

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

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

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

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

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

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

1.4 IDVM 研究進(jìn)展

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

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

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

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

圖6 三維頂蓋驅(qū)動(dòng)流問題的計(jì)算結(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

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

2.2 算例驗(yàn)證

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

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

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

圖9 墻壁束縛瑞利流在Kn = 0.05 和t= 1.5 時(shí)水平中心線和垂直中心線流動(dòng)變量分布圖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 時(shí)上下表面流動(dòng)變量分布圖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 模型方程求解算法，展示了它們?cè)诘退俚礁叱曀?、連續(xù)到稀薄全流域范圍的應(yīng)用。同時(shí)，本文還將筆者前期發(fā)展的定常隱式IDVM 擴(kuò)展到了非定常情形，并通過非定?？缌饔騿栴}進(jìn)行了驗(yàn)證。這類方法由于直接演化分布函數(shù)，可以有效避免粒子類方法的統(tǒng)計(jì)噪聲問題，同時(shí)可以采用傳統(tǒng)計(jì)算流體力學(xué)中的加速技術(shù)來提高計(jì)算效率。在微電子機(jī)械和真空技術(shù)等低速跨流域問題以及不含離解和電離化學(xué)非平衡效應(yīng)的高速跨流域問題中，該類算法已經(jīng)展現(xiàn)出一定的優(yōu)勢(shì)。但是對(duì)于航天領(lǐng)域涉及的高溫高速和真實(shí)氣體效應(yīng)的非平衡流動(dòng)，該類算法也還存在許多丞待解決的問題，需要持續(xù)不斷地改進(jìn)完善。

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

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