離散速度方法及其在跨流域問題中的應用研究進展

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

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

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

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

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

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

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

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

1 Boltzmann模型方程和DVM框架

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

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

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

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

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

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

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

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

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

1.2 GKUA 研究進展

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

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

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

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

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

圖2 返回艙以第一宇宙速度7.5 km/s 再入跨流域（120～70 km）氣動模擬軟件系統實時計算與姿態配平繞流壓力分布[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 條件下強氣動力/熱環境致結構溫度變化和穩態變形[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 模型方程和對應的宏觀伴隨方程。宏觀伴隨方程實際上對應于Boltzmann 模型方程在速度空間求守恒矩

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

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

由方程（18）可知，當分子遷移時間遠小于分子平均碰撞時間時（t?τ），單元界面分布函數中初始分布f(xij-ξαt，ξα，0)占主導，表現為自由分子流情形的完全自由遷移現象；而當分子遷移時間遠大于分子平均碰撞時間時（t?τ），單元界面周圍的平衡態分布fS(xij-ξα(t-t′)，ξα，t′)占主導，表現為連續流情形的分子頻繁碰撞現象。由于實際計算中分子的遷移時間即為網格尺度對應的時間步長Δt，方程（18）等效于在網格尺度上的流動建模，將主導稀薄流的分子自由遷移機制和主導連續流的分子碰撞機制的權重與網格尺度有機聯系起來。由于通量計算時分子的遷移和碰撞過程相互耦合，UGKS 的推進時間步長不受限于當地分子的平均碰撞時間，且其網格尺度也不受限于當地分子平均自由程，從而使得該方法在全流域均可獲得準確高效的計算結果。

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

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

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

1.4 IDVM 研究進展

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

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

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

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

圖6 三維頂蓋驅動流問題的計算結果[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 及其應用

2.1 非定常隱式IDVM

為了實現非定常跨流域流動問題的準確高效求解，本文將先前發展的隱式定常IDVM 進一步拓展到了非定常情形?；陔p時間步LU-SGS 迭代的思想，在Boltzmann 模型方程和宏觀伴隨方程中同時增加了偽時間導數項

2.2 算例驗證

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

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

為了與文獻條件一致，本文將該算例的克努森數取為Kn= 0.05。分子速度空間采用包含28×28 個離散點的Gauss-Hermite 積分，位置空間采用12 800 個均勻四邊形非結構網格進行離散。無量綱的外迭代時間步長取為Δt=0.01。圖9 展示了t= 1.5 時刻水平和垂直中心線溫度、x方向速度分量和y方向速度分量分布，當前計算結果與UGKS結果較好吻合。另外，圖10 展示了t= 1.5 時刻上下平板的壓力、熱通量和切應力分布，當前結果也與GUKS 結果吻合良好。該算例驗證了提出的方法在非定常跨流域流動的計算能力。

圖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 總結與展望

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

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

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