基于修正N-S方程本構關系的氣體動理論耦合方法

皮興才，李志輝，2，*，彭傲平，吳俊林，蔣新宇

(1. 中國空氣動力研究與發展中心 超高速空氣動力研究所，綿陽 621000；2. 國家計算流體力學實驗室，北京 100191)

0 引 言

近年來，飛行高度在40～70 km、馬赫數在5～20的近空間飛行器的研發受到重視[1]。在此環境下飛行普遍面臨著局部稀薄效應導致的連續介質假設失效的問題。因此，考慮氣體分子輸運統計分布函數演化特性的介觀方法獲得了廣泛的關注。Boltzmann方程是介觀方法的理論基礎[2]。在研究應用中，為了定量化表征處理Boltzmann方程復雜的碰撞項，很多簡化的模型方程被提出，包括BGK模型方程[3]、Shakhov模型方程[4]、ES-BGK模型方程[5]、Rykov模型方程[6]等，所能模擬的氣體流動介質也從最初的單原子理想氣體發展到多原子、多組分，考慮轉動、振動激發及熱化學非平衡效應的跨流域空氣動力學模擬體系日趨完善[7-17]。

模型方程是描述氣體分子速度分布函數關于時間、空間及速度相空間七個維度的偏微分方程。由于方程的數值求解涉及到對七個維度的離散，提高計算效率一直是亟待解決的難題。針對BGK模型、ESBGK模型方程，Mieusses[18]通過構造具有保守恒特性的離散平衡分布函數解析表達形式，建立了嚴格守恒的數值方法，使得該方法能以較少的離散速度點實現方程求解。Chen、Huang等[19-20]同樣立足于減少離散速度空間的大小，采用叉樹數據結構對UGKS的離散速度空間進行優化以提高計算效率。進一步地，2016年，江定武[21]針對顯式UGKS計算效率低下的問題，利用CFD領域的經典數值方法構造了隱式格式，并采用MPI+OpenMP并行計算策略以改善UGKS的計算效率及收斂速度。在耦合方法方面，2012年Alessandro等[22]采用DSMC方法對單原子ES-BGK方程進行數值求解，并通過區域剖分實現了與Euler方程的耦合求解，并將其應用到一維間斷問題分析模擬。2014年Rovenskaya等[23]采用有限差分法與離散速度坐標法結合，開展了Shakhov模型方程與NS模型方程的耦合求解，并將其應用于低速稀薄氣體流動分析，實現了任意壓比下的薄板稀薄繞流模擬。最近，基于多重網格技術及采用宏觀方程預測平衡分布函數的全隱UGKS成功應用于跨流域稀薄流模擬，大幅提升了傳統UGKS的計算效率[24]。

在氣體動理學理論框架下，N-S方程可以通過Boltzmann方程或模型方程的Chapman-Enskog展開一階近似。這體現了N-S方程與模型方程的理論可銜接性，也體現出N-S方程的質量、動量及能量守恒的基本特性即使在稀薄非平衡效應加重的近連續流區依舊成立，而對于稀薄效應更為嚴重、氣體分子速度分布函數呈現顯著非平衡特征的流動，N-S方程所遵從的本構關系難以正確表征模擬時，通常采用滑移邊界可解決邊界條件失效問題，但針對本構關系失效卻是難以修正。這也是二階滑移邊界條件，如Cercignan模型[25-26]、Deissler模型[27]等的模擬精度相較于一階滑移邊界條件的提升有限的根本原因[28-29]。為了描述更復雜的本構關系，肖洪[30]在廣義流體力學方程（GHE）的基礎上提出了非線性耦合本構關系（NCCR），并在激波結構、二維高超聲速稀薄繞流等問題模擬中獲得成功。Burnett方程、Grad矩方程同樣立足于通過理論建模的方式獲得非守恒量—應力張量及熱流的輸運關系，也取得了成功的應用[31]。

為了從數值層面構建本構關系，拓寬N-S方程在稀薄氣體動力學方面的模擬能力，規避復雜的宏觀非守恒量輸運方程帶來的求解難題，本研究以多原子氣體ES-BGK模型方程為基礎，采用氣體動理論統一算法（GKUA）[7-17]描述非線性本構關系及壁面速度滑移、溫度跳躍，修正N-S方程的線性本構關系及滑移邊界條件，拓展其在過渡流、滑移流域的應用能力。在后續內容中，本文首先通過Chapman-Enskog基于Maxwell平衡態分布漸近展開給出ES-BGK方程當地平衡態分布函數、本構關系的一階近似，其次，構造基于N-S方程本構關系修正及滑移邊界耦合的數值求解方法，并通過可壓縮平板邊界層流動、圓柱繞流經典案例來進行驗證。

1 理論模型

經典熱力學理論認為分子內能可以通過分子的自由度及溫度來表征。基于該理論，Brull等[5,16]在單原子ES-BGK模型的基礎上通過引入內能自由度eint作為分布函數的自變量，即分布函數f=f(t,x,ξ,eint)構造了一種適用于多原子分子氣體的ES-BGK類模型方程：

式中，

其中，δ為轉動自由度；Z為轉動松弛參數；I為單位矩陣；b∈[-0.5,1) ；氣體分子熱運動速度C是氣體分子隨機速度ξ與宏觀速度u之差，其平衡態分布函數為：

式中，

對ES-BGK模型方程進行基于Maxwell平衡態分布的Chapman-Enskog漸近分析。基于量綱原理將式(1)改寫為如下形式[26,32]：

式中，無量綱參數 ε是與Knudsen數同階的量，用于衡量松弛碰撞、對流及擴散等物理過程的時間尺度。

對分布函數及當地平衡態分布函數進行漸近展開，具有如下形式：

由Chapman-Enskog漸近分析需要遵循的約束條件：

可知f(0)=fM。將Ttr及Tint同樣進行多尺度展開：

式中，ΔTtr及ΔTint表示非平衡效應引起的平動及轉動溫度增量，它們滿足如下關系：

對于非守恒量—應力張量及熱流進行多尺度展開：

式中，

為了獲得當地平衡分布函數G的一階近似表達式，對式(2)中 det(Γ) 及 Γ-1進行多尺度展開分析。根據前述關系，整理Γ表達式可得：

式中：

將式(24)帶入式(29)，并保留d et(Γ) 及Γ-1的Taylor級數到O(ε)可得：

針對量熱完全氣體不考慮ΔTtr及 ΔTint對結果的影響，對當地平衡態G進行Taylor開展，其O(1)及O(ε)階項如下所示：

進一步地，將f及G的展開式(18)及式(19)代入式(17),可知：

等式左端等于[5]：

式中，

可得：

因此，利用式(9)可知分布函數的Chapman-Enskog漸近展開一階近似表達式為：

式中，

可見，ES-BGK模型的應力張量σ及 熱流q的一階漸近展開結果與N-S方程的本構關系相同。由此可知，ES-BGK方程與N-S方程在連續流域的物理描述上具有理論一致性。進一步說明，N-S方程在描述稀薄氣體流動問題時，其質量、動量及能量守恒的基本特性依舊有效，其線性本構關系是方程失效的核心。這也為研究者拓展N-S方程的應用范圍提供了理論支持。基于上述關系，本研究提出一種基于修正NS方程本構關系的耦合氣體動理論的方法，以保證耦合方法在質量、動量及能量守恒性，拓展N-S方程在局部稀薄效應氣體流動問題模擬的能力。

2 數值方法

本耦合方法基于課題組氣體動理論統一算法(GKUA)及開源代碼OpenCFD-EC2D的NS有限體積求解器構建。作為原理驗證研究，本文中耦合區域通過計算調整給定。其中，GKUA求解需要進行N-S方程本構關系修正的區域—壁面附近與平均分子自由程相當的區域，N-S求解器對全流區進行守恒方程求解，如圖1所示。

圖1 求解域分區示意圖Fig. 1 Schematic view of the computational domain decomposition

氣體動理論耦合策略大致可分如下幾步：

1）通過有限體積法全流域求解N-S方程；

2）通過GKUA方法求解壁面區域，耦合邊界處的分布函數根據特征線n指向給定分布函數的邊界值。其中， (n·ξ＞0)的離散分布函數采用外推方式給定；(n·ξ≤0)的離散分布函數由N-S方程計算給出的宏觀物理量，通過式（38）計算離散分布函數給定。壁面分布函數按完全漫反射條件給出；

3）求解GKUA計算域的應力張量 σGKUA、熱流矢量qGKUA及 壁面速度滑移uslip及溫度跳躍Tjump條件；

4）利用獲得的 σGKUA、qGKUA替換GKUA域的NS計 算 結果 σNS及qNS， 利 用uslip及Tjump給 定N-S求解域邊界條件；

5）以修正后的本構關系及壁面邊界進行下一時間步N-S方程及GKUA求解。

為了降低對計算資源的需求，采用離散速度坐標法直接數值求解Boltzmann模型方程通常對分布函數的自變量進行約化[7-8]，引入變量：

對于二維流動問題，引入下式對ES-BGK模型方程進行進一步的約化：

在空間格式構造方面，本文采用二階NND格式對N-S方程及ES-BGK模型方程進行界面重構[7-9]，采用Steger-Warming格式進行界面通量計算，并且采用LU-SGS隱式格式對時間項進行離散以提高效率。

3 算例驗證

為了驗證構建的基于修正本構關系的氣體動理論耦合方法，本文分別進行了可壓縮平板邊界層流動模擬及稀薄圓柱繞流模擬。其中，DSMC計算結果由DS2V軟件[33]計算得到。

3.1 可壓縮平板邊界層

本研究首先模擬了連續流域不同馬赫數可壓縮平板邊界層流動，以驗證介觀方法及傳統CFD在連續流域模擬的一致性。計算域x×y= [0,1 m]×[0,0.5 m]，網格量：100×100，壁面第一層網格高度：0.000 1 m。來流狀態為：高度H= 55 km，靜壓p= 42.5 Pa，靜溫T=206.7 K，密度ρ= 5.68×10-4kg/m3，馬赫數Ma= 2、4、8，各 狀 態 的 雷 諾 數Re分 別 為：2.22×104、4.45×104、8.89×104。壁面溫度等于來流靜溫。分別采用N-S及GKUA進行了Ma= 2、4、8流場的數值模擬，并將其速度場和溫度場結果與文獻[34]給出的可壓縮層流邊界層方程的理論結果進行了對比，見圖2。根據邊界層相似理論，其橫坐標 η定義如下：

對比可知，采用單一方法計算連續流域可壓縮平板邊界層流動給出的速度場及溫度場結果與文獻吻合良好，驗證了本文的GKUA及N-S的可靠性。

隨著高度的增加，稀薄效應將逐漸顯現，壁面附近將出現速度滑移及溫度跳躍。本文對來流條件：高度H= 80 km，靜壓p= 1.05 Pa，靜溫T= 198.6 K，密度ρ= 1.845 8×10-5kg/m3，馬赫數Ma= 2，雷諾數Re= 7.9×102的可壓縮平板邊界層流動進行了模擬，壁面溫度等于來流靜溫。數值方法包括GKUA方法、無滑移NS方法、DSMC方法及采用前文所述理論構建的耦合方法—“Kinetic Boundary”。其中，耦合方法的本構關系修正區域給定為壁面附近2倍平均分子自由程高度范圍內的區域。由圖3所示壁面物理量對比可以看出，構建的耦合方法的壁面壓力計算結果與GKUA、N-S計算結果吻合良好，略大于DSMC結果；耦合方法計算的溫度跳躍結果和GKUA吻合良好，并和DSMC基本吻合；滑移速度及壁面剪切應力與GKUA、DSMC計算結果吻合良好。并且，由于平板前緣的速度滑移現象較平板后部嚴重，采用耦合方法、GKUA及DSMC方法計算的壁面剪切應力在平板前緣較N-S小；在平板后緣，各種結果相互吻合良好。

圖2 不同馬赫數下平板邊界層模擬結果Fig. 2 Simulation results of a flat-plate boundary layer at different Mach numbers

圖3 Ma = 2平板邊界層壁面物理量Fig. 3 Distributions of physical quantities along the streamwise direction on the surface of a flat-plate boundary layer at Ma = 2

3.2 稀薄圓柱繞流

為了進一步驗證采用氣體動理論耦合方法的有效性，本文以高度H= 90 km，靜壓p= 0.183 Pa，靜溫T= 186.8 K，密度ρ= 3.416 31×10-6kg/m3，馬赫數Ma=2、4，其雷諾數Re分別為150、300作為來流條件，模擬了圓柱超聲速稀薄繞流，并對流場及壁面物理量進行了對比分析。壁面溫度等于來流靜溫。其中，耦合方法的N-S方程本構關系修正區域給定為壁面附近2倍平均分子自由程高度范圍內的區域。圓柱半徑1 m，遠場半徑10 m，網格量：周向60×徑向80，壁面第一層網格高度0.000 5 m。

圖4 Ma = 2圓柱稀薄繞流流場Fig. 4 Rarefied flow fields around a circular cylinder at Ma = 2

圖5 圓柱稀薄繞流壁面物理量Fig. 5 Distributions of physical quantities along the circumferential direction on the surface of a circular cylinder in rarefied flow

圖4給出了氣體動理論耦合方法計算的流場馬赫數、靜壓、密度云圖的與DSMC計算結果的對比情況。圖5給出了圓柱壁面物理量計算結果對比，其中，橫坐標θ= 180°對應于圓柱迎風面駐點。從圖4可以看出，本文構建的耦合方法獲得的流場與DSMC獲得的流場在圓柱迎風側吻合良好，在圓柱背風回流區存在一定的差異。結合圖5的壁面剪切應力及滑移速度對比圖可看出，本文構建的耦合方法獲得的回流區小于DSMC獲得的分離區，是導致圓柱背風區流場云圖存在差異的主要原因。

從圖5可知，各方法獲得的壁面壓力吻合良好。壁面速度滑移結果顯示：DSMC、GKUA及氣體動理論耦合方法在圓柱迎風一側(θ∈[90°,180°])吻合良好；在背風區，由于計算的分離區的大小差異導致結果存在一定偏差，但總體規律基本吻合。進一步地，雖然圖中給出的Ma= 4流動無量綱壁面速度滑移小于Ma= 2流動，但其有量綱的絕對量是大于Ma2流動的壁面速度滑移。壁面剪切應力對比結果顯示，DSMC、GKUA及氣體動理論耦合方法計算結果基本吻合，而無滑移NS計算結果和其他結果存在一定差異。并且，該差異整體上隨壁面滑移速度的增大而增加，流動分離點(θ≈ 40°)附近的流動較復雜，其差異最大。壁面熱流計算結果顯示無滑移N-S結果略大于其他方法計算結果，但各方法計算結果總體吻合。

4 結 論

本研究通過對多原子氣體ES-BGK方程進行Chapman-Enskog漸近分析，構建了分布函數、應力張量、熱流矢量的漸近一階表達式，證明了ES-BGK方程與NS方程在連續流域的物理描述上具有理論一致性。在此基礎上，提出通過對稀薄效應明顯的流動區域進行應力張量及熱流矢量修正，構造了氣體動理論耦合方法，避免了離散速度坐標法引入的質量、動量及能量守恒型誤差對耦合方法計算結果的影響。

進一步地，本研究通過引入約化變量，構造了分布函數一階表達式的二維約化形式，并通過有限體積法對二維耦合方程進行數值求解。通過二維可壓縮平板邊界層及圓柱繞流算例的數值模擬及結果對比，初步驗證了本研究提出氣體動理論耦合方法的有效性，預期進一步研究發展，可得到一套適于近連續滑移過渡流區高效可靠的氣體動理論耦合方法。

致謝：感謝中國科學院力學研究所李新亮研究員團隊提供的幫助。