基于宏觀方程數值本構關系的氣體動理論加速收斂方法*

皮興才 朱煉華 李志輝3)? 彭傲平 張勇豪

1) (中國空氣動力研究與發展中心超高速所, 綿陽 621000)

2) (英國思克萊德大學機械與航空工程系, 詹姆斯·維爾流體實驗室, 格拉斯哥 G1 1XJ)

3) (國家計算流體力學實驗室, 北京 100191)

1 引 言

在超聲速稀薄氣體動力學數值計算領域, 基于 Boltzmann方程[1,2]或 Shakhov方程[3]、ESBGK方程[4]等模型方程的數值方法獲得了廣泛的關注, 如基于離散速度坐標法[5,6]的氣體動理論統一算法 (gas kinetic unified algorithm, GKUA)[7?10]、統一氣體動理學格式 (unified gas kinetic scheme,UGKS)[11], 以及漸近保持格式[12,13]、矩方法[14]、離散 Boltzmann 方法 (discrete Boltzmann method,DBM)[15,16]等. 特別是 GKUA, 在大型航天器、飛船返回艙再入過程模擬等工程問題中獲得了成功應用[17?21]. 基于離散速度坐標法的各類算法需要求解包括時間、位置空間及速度相空間的高維偏微分方程組, 導致其計算量大. 為了使該類方法獲得更廣泛應用, 提升計算效率是核心研究主題.

在提高模型方程數值計算效率方面, 目前的研究工作主要分為如下幾類: 1)減少速度相空間的離散速度點; 2) 提高數值格式的穩定性, 增大時間步長; 3) 與其他類型方程進行物理空間分區耦合求解, 減少氣體動理論方法的計算區域; 4) 構造并同步求解宏觀守恒方程, 利用宏觀守恒方程加速信息傳遞過程. 其中, 第一類(如文獻[22])通過構造新的離散平衡態速度分布函數解析表達形式, 建立了嚴格保守恒的數值方法, 使得該方法能以較少的離散速度點實現方程求解. 文獻[23]采用叉樹數據結構對UGKS的離散速度空間進行優化以提高計算效率; 第二類(如文獻[24])對模型方程進行了隱式處理, 并引入演化時間平均界面通量, 通過對控制方程矩陣進行近似LU分解實現了隱式UGKS,還包括GKUA基于LU-SGS構造的有限體積隱式數值格式, 成功應用于航天器再入解體物的繞流分析[25?27]; 第三類 (如文獻 [28])開展了 Shakhov 模型方程與N-S方程的耦合求解, 并將其應用于低速稀薄氣體流動分析, 實現了任意壓力比下的狹縫稀薄流模擬; 第四類 (如文獻 [29])對求解線化Boltzmann模型方程的離散速度坐標法進行Fourier穩定性分析, 獲得了方法在連續流域收斂變慢的原因, 并進一步構造了不同階數的Hermite矩方程用于算法加速收斂. 另外, 文獻[30,31]采用UGKS獲得網格界面的守恒量數值通量, 并用數值通量構造宏觀演化方程, 通過同步求解構造的宏觀演化方程獲得宏觀量并用于預測當地平衡態速度分布函數, 大幅度提升了傳統UGKS的效率.

通過Chapman-Enskog漸近分析方法對Boltzmann模型方程進行一階展開可獲得N-S方程. 另一方面, 對Boltzmann模型求碰撞不變量的矩可獲得關于質量、動量及能量的宏觀守恒方程. 對比可知, 宏觀守恒方程相較于N-S方程的本質區別在于本構關系, 即非守恒量——應力、熱流與守恒量的關系. 目前, 基于宏觀輸運方程的稀薄氣體動力學理論建模有相當多的工作都聚焦在建立準確、數學性質好、可數值計算的應力、熱流輸運方程或本構關系方面, 如 Burnett方程的正則化[32]、Grad13方程的正則化[33]以及在廣義流體力學方程的基礎上提出了非線性耦合本構關系(nonlinear coupled constitutive relations, NCCR)[34,35]. 這些理論模型成功開展了一維激波結構、二維及三維高超繞流等問題的模擬分析[36]. 值得關注的是, 通過數值的方式封閉非守恒量輸運方程也是描述近空間飛行出現的非線性稀薄效應的一個新思路. 例如, 文獻[37]基于求解輻射輸運方程的加速算法思想[38]構造了求解線化Boltzmann方程的通用協同迭代格式 (general synthetic iteration scheme, GSIS).該方法通過協同迭代宏觀輸運方程及Boltzmann模型方程, 實現了庫葉特流動、方腔稀薄流動等的快速收斂. GSIS方法的核心在于通過對Boltzmann模型方程求高階矩, 建立宏觀守恒量及非守恒量的輸運方程. 通過Chapman-Enskog展開將非守恒量宏觀輸運方程的應力、熱流高階項剝離出來, 并利用數值求解Boltzmann模型方程給定高階項的值, 由此實現非守恒量宏觀輸運方程的封閉并加速收斂過程. 最近, 文獻[39]將GSIS拓展到非線性氣體動理論方程的加速收斂算法應用研究中, 并實現了低速流動問題100迭代步收斂,Kn= 0.01的超聲速圓柱繞流問題提升10余倍計算效率的效果. 為了更全面地驗證協同迭代方法的有效性,基于相同的思路, 本文在氣體動理論統一算法(GKUA)框架下, 構造了耦合宏觀方程本構關系的加速收斂方法. 在控制方程方面, 采用Boltzmann模型方程描述強非線性稀薄流動問題, 且簡化了宏觀輸運方程中應力、熱流高階項的構造方式; 在數值格式方面, 改進了具有捕捉強間斷能力的隱式氣體動理論數值格式; 并且, 利用多塊對接網格技術提高處理復雜問題的能力. 本文的研究工作將進一步驗證方法的有效性并展示了其向工程應用領域拓展的能力.

本文第2節將構建加速收斂方法的理論模型;第3節將對數值求解方法進行描述; 第4節通過方腔流、超聲速圓柱繞流及雙圓柱干擾繞流模擬, 對本文理論算法模型進行驗證分析; 第5節將對研究內容進行總結.

2 理論模型

本研究以Shakhov模型方程為例, 構造加速收斂方法的理論模型. 此構造方式同樣適用于其他模型方程. Shakhov模型方程的速度分布函數f(t,x,ξ)滿足關系[3]

其中t為時間,x為物理空間,υ為碰撞頻率,ρ為密度,T為溫度,R為氣體常數,q為熱流, 壓力p=ρRT. 對于單原子氣體, 普朗特數Pr=2/3 .氣體分子熱運動速度C是氣體分子隨機速度ξ與宏觀速度U之差. 密度、溫度、宏觀速度及熱流的定義如下:

以二維流動問題為例, 構造求解模型方程的有限體積隱式格式. 引入約化速度分布函數對方程(1)進行速度相空間降維,

采用離散速度坐標法對約化速度分布函數方程進行速度空間數值離散, 可得

在物理空間對方程(3)構造格心型有限體積格式, 可得

式中I為物理空間網格單元索引, 符號“?”表示物理量在網格單元的均值,N(I) 為網格I的鄰近網格 單 元 索 引,ξJ,n為ξJ在 界 面(I,K) 的 法 向 分 量,SI,K為界面 (I,K) 的面積,為分布函數在界面的重構量,?I為第I網格單元的體積 (二維: 面積). 分布函數的重構可以采用經典的CFD重構格式, 如NND格式, 也可采用考慮氣體分子碰撞遷移物理過程, 具有漸近保持特性的重構方式, 如DUGKS的重構方式[40,41]. 為了簡化問題, 本研究采用經典的CFD重構格式進行界面重構. 界面數值通量的計算采用具有迎風特性的Steger-Warming矢通量分裂方法:

顯式地求解方程(4)將面臨對流項數值穩定性條件——CFL條件, 以及當流動趨近于連續流時的剛性源項導致的時間步長限制. 采用LUSGS進行隱式格式構造可以顯著增加時間步長,提升計算效率[7,19,26]. 方程(4)在n+1 時刻具有如下形式:

在等式兩邊同時加上n時刻殘差:

的控制方程:

方程(9)左端項的差量的界面重構方式不影響計算結果. 因此, 一階迎風格式是最好的選擇, 即依據離散速度坐標點的方向, 選擇界面迎風側的網格單元均值作為界面值. 由此可獲得關于差量的線性方程組, 可方便地使用LU-SGS方法數值求解上述方程組, 獲得n+1 時刻的差量解進而通過數值積分獲得n+1 時刻的流場宏觀物理量.

這樣的隱式處理方式提高了算法的計算效率,然而隨著努森數Kn的降低, 其迭代步數依舊會增加, 收斂速度變慢. 文獻[29]對求解線化Boltzmann模型方程的離散速度坐標法開展誤差傳遞分析, 得出努森數Kn降低將導致誤差傳遞矩陣的譜半徑逐漸趨近于1. 這是導致收斂變慢的主要原因. 另外, 有研究認為, 隱式格式構造中引入的簡化假設對收斂具有負面影響[30].

為了提高算法的收斂速度, 借助輻射輸運方程數值模擬采用的協同迭代求解的思路[37?39], 本研究通過構造并協同迭代求解一套與Boltzmann模型方程相容的宏觀守恒方程來實現流場擾動加速傳遞及耗散. 通過對方程(1)求碰撞不變量φ=的矩, 可得模型方程對應的質量、動量及能量守恒的宏觀守恒方程如下:

CV為定容熱容.

宏觀守恒方程組(10)并不封閉. 為了實現輸運方程的理論表達形式封閉并具有良好的數學性質, 已有研究者提出了一系列的理論形式, 如正則化Burnett方程、R-13方程以及NCCR模型[32?35].盡管做了相應的簡化, 這些理論形式依舊相當復雜, 數值穩定性仍然需要進一步優化. 近年來數值求解Boltzmann模型方程的氣體動理論方法同樣體現出巨大的潛力. 利用氣體動理論方法數值求解非線性應力、熱流, 并結合經典CFD在計算效率、數值穩定性方面的優勢, 從數值層面描述非線性本構關系并封閉宏觀守恒方程組(10)無疑是新的途徑.

在Chapman-Enskog漸近分析方法中, 氣體速度分布函數f和非守恒量σij及qi具有多尺度特性. 它們的多尺度展開形式為

對Shakhov模型方程進行Chapman-Enskog漸近分析可得:

式中,μ為黏性系數,κ為熱傳導系數. (14) 式及(15)式即為N-S方程的線性本構關系式. 為了保持宏觀守恒方程與模型方程的一致性, 應力及熱流的高階項必須在計算中予以考慮. 值得一提的是,在數值計算過程中, 來自壁面的擾動是通過邊界層的剪切作用向主流區及遠場傳遞并在沿途進行耗散. 耗散(包括物理耗散及數值耗散兩個部分)對流場收斂及數值穩定性都有積極作用. 在宏觀方程的求解過程中, 保持一階應力、熱流項表達式具有的耗散性質可更好地利用宏觀方程擾動傳遞快的優勢, 加速收斂. 對于超聲速繞流, 增加宏觀方程的耗散對數值穩定性更為重要.

基于上述考慮, 將宏觀方程的應力及熱流項分解為兩部分, 一階線性項及高階項[37,39]:

并且, 高階項 H oTσ及 H oTq獲取方法如下:

采用迭代求解宏觀方程(10)獲得的宏觀量對方程(6)中的當地平衡態速度分布函數及碰撞頻率進行預測, 并重新構造隱式格式, 可將方程(9)改進為如下形式:

3 數值方法

對模型方程進行無量綱化, 引入特征長度L,特征溫度T0, 特征黏性系數μ0, 特征速度Cm∞=特征數密度n0, 特征壓力特征熱流及特征速度分布函數Shakhov模型方程(1)的黏性系數采用如下指數律:

本研究中ω=0.81 . 因此, 無量綱碰撞頻率等于:

為了簡化, 省去無量綱量的上標“?”. 若無特別說明, 后續量都為無量綱量.

基于第2節介紹的理論模型, 耦合宏觀方程本構關系的加速收斂方法使用開源代碼OpenCFDEC2 D[42]的核心代碼進行宏觀方程數值格式構建;使用本課題組GKUA核心代碼進行氣體動理論相關算法構建; 將二階NND格式用于Shakhov模型方程及宏觀方程的界面物理量重構; 界面通量采用Steger-Warming矢通量分裂方法計算; 使用完全漫反射邊界條件確定壁面反射速度分布函數. 其求解流程如下.

1)確定離散速度空間及數值積分方法. 為保證離散速度空間能體現整個流場的速度分布形態,對速度空間進行截斷通常需要滿足[43]

2) 確定第k時間步的速度分布函數fk及宏觀流動物理量Wk.

3) 通過 (18)式及 (19)式獲得高階應力項HoTσ及高階熱流項 H oTq.

4) 將第3步獲得的高階應力項 H oTσ及熱流項HoTq數值解代入(16)式及(17)式, 求解由方程(10)式、(16)式及(17()式構成)的方程組, 獲得預測的宏觀流動物理量Wk+1?. 即將高階應力項 HoTσ及熱流項 H oTq視作宏觀方程的固定源項, 迭代求解直至宏觀方程收斂到特定精度或達到指定的迭代步數. ()

5) 采用第4步獲得的宏觀流動量Wk+1?演化更新當地平衡態速度分布函數及碰撞頻率, 并采用LU-SGS方法求解構造的隱式方程(20), 獲得第k+1步的離散速度分布函數fk+1及宏觀流動物理量Wk+1.

6)判斷是否達到收斂標準. 本文采用全流場υk與υk+1的相對偏差均值作為收斂判斷標準:

式中,I為物理空間網格索引, N umcell為物理空間總網格數, ?t為時間步長. 若還未達到收斂標準,繼續進行第2—5步迭代.

4 算例驗證

4.1 方腔流動

對于方腔流動, 主流區黏性特性能否準確描述對流場計算結果, 尤其是溫度場的影響尤為明顯.因此, 被廣泛用于算法驗證. 為了與已有研究成果對比, 本文以方腔邊長為參考長度分別計算了努森數Kn= 1.000, 0.075, 0.010 的方腔稀薄流動. 努森數Kn定義為[44,45]:

計算區域為邊長L= 1的方塊區域, 所有壁面的無量綱溫度均為Tw=1.0 , 上壁面無量綱速度Uw=0.15, 其余壁面無量綱速度為0. 物理空間網格為 65 × 65, 在 [–5.67, 5.67]× [–5.67, 5.67]的速度相空間采用32 × 32個Gauss-Hermite離散速度坐標點. 圖1繪出了各狀態的溫度場對比情況(“GKUA”表示常規 GKUA 結果; “Coupled”表示本文提出的耦合加速收斂方法結果). 通過對比可以看出兩者良好符合, 反映了方腔內的流動從左側的低溫膨脹過渡到右側的高溫壓縮的變化過程, 這樣的過程對稀薄效應越強的流動狀態表現越突出.為了定量分析, 圖2繪出了中心線的速度分布計算結果與常規GKUA及文獻[45]結果的對比情況.計算結果準確反映了壁面速度滑移隨稀薄程度的增加而增大以及整個流域的非線性分布規律.

圖 1 方腔流溫度分布計算結果 (a) Kn = 1 ; (b) Kn =0.075; (c) Kn = 0.01 (GKUA: 彩 色 背 景 及 黑 色 實 線 ;Coupled: 紅色虛線)Fig. 1. Temperature distribution in cavity flow: (a) Kn = 1;(b) Kn = 0.075; (c) Kn = 0.01 (GKUA: coloured background and black solid lines; Coupled: red dashed lines).

圖 2 方腔中心線上 的速度場 (a) Kn = 1; (b) Kn =0.075; (c) Kn = 0.01Fig. 2. Velocity profiles at the central lines of the cavity:(a) Kn = 1; (b) Kn = 0.075; (c) Kn = 0.01.

圖3給出了耦合加速收斂方法計算收斂曲線與常規GKUA對比情況. 可知, 在小努森數Kn情況下本方法能大幅度減少計算收斂的迭代步數. 隨著努森數Kn增大, 加速收斂效果逐漸減弱. 表1列出了本文提出的耦合加速收斂算法的迭代步數與常規 GKUA 的迭代步數對比情況. 其中,Kn=0.01的加速比達到47倍. 分析認為, 隨著努森數Kn增大, 擾動傳遞由宏觀對流主導逐漸轉換為由分子碰撞擴散效應主導, 這削弱了借助宏觀方程加速收斂的優勢. 需要指出的是, 由于該算例采用的宏觀方程求解器為適于高速可壓縮流動問題的數值格式, 其在低速弱可壓的方腔流模擬中收斂較慢. 因此, 從計算耗時角度看, 耦合加速收斂方法的優勢受到嚴重影響. 但是, 從減少氣體動理方法的迭代步數來看, 該方法的效果可觀, 驗證了耦合本構關系思路的有效性.

圖 3 耦合加速收斂方法與常規GKUA的計算收斂歷史Fig. 3. The convergence history between coupled acceleration method and the conventional GKUA.

表 1 耦合加速收斂方法與常規 GKUA 的收斂情況對比Table 1. Convergence comparison between the conventional GKUA and the coupled acceleration method.

4.2 超聲速圓柱繞流

為了進一步探討本文發展的耦合加速收斂方法對超聲速稀薄繞流問題的模擬能力, 本文計算了來流馬赫數Ma= 5的圓柱超聲速稀薄繞流. 相較于低速流動, 超聲速繞流存在著強剪切邊界層、激波以及流動分離, 將導致速度分布函數嚴重偏離Maxwell平衡態分布. 激波結構、壁面物理量的結果對比可以有效反映耦合加速收斂方法描述強非平衡物理特征的準確性. 為了與已有研究成果對比, 以圓柱半徑為參考長度, 本文計算了來流努森數Kn= 0.01, 0.10 的超聲速稀薄繞流. 努森數Kn定義如下[31,44]:

圖 4 Kn = 0.01 圓柱繞流流場對比 (a) 壓力; (b) 溫度;(c) 馬赫數 (GKUA: 彩色背景及白實線; Coupled: 紅色虛線)Fig. 4. (a) Pressure, (b) temperature, (c) Mach number distribution around cylinder for Kn = 0.01 (GKUA: coloured background and white solid lines; Coupled: red dash lines).

無量綱壁面溫度等于來流靜溫, 即Tw=1.0 .本文的計算區域為: 圓柱半徑R= 1; 計算域遠場半徑Rf= 11; 物理空間網格: 周向 82 × 徑向 80,第一層網格高度 0.0002, 并按雙曲正切函數tanh 增長. 在 [–15, 15]× [–15, 15]的速度相空間采用100 × 100 個Gauss- Legendre離散速度坐標點.圖4繪出了Kn= 0.01案例的壓力場、溫度場和馬赫數場的結果及其與常規GKUA果的對比. 結果顯示, 采用本文提出的耦合加速收斂方法計算獲得的流場與常規GKUA計算獲得的流場符合良好,包括激波厚度、尾流區大小等細節. 圖5繪出了Kn= 0.1案例的壓力場、溫度場及馬赫數場的結果及其與常規GKUA果的對比. 同樣, 兩種方法獲得的結果符合良好. 并且, 對比Kn= 0.01 及Kn= 0.1案例的流場結果可知, 增大稀薄程度可使激波層、邊界層變厚, 降低物理量的梯度, 有利于減小不同格式的數值黏性差異導致的結果差異.圖6給出了本文的耦合加速收斂方法獲得的壁面物理量與常規GKUA及DSMC獲得的壁面物理量對比情況. 其中,Kn= 0.01 案例的 DSMC 計算結果由 DS2 V 軟件[46]計算得到,Kn= 0.1 案例的計算結果源自參考文獻 [31]. 橫坐標q= 180°對應于圓柱迎風面駐點. 可以看出, 三種方法獲得的壓力符合良好, 并且氣體稀薄程度增加將導致激波強度變弱, 進而增大壓力系數; 熱流結果總體符合,但頭部駐點的差異略有增大; 三種方法獲得的壁面剪切應力符合良好. 圖7給出了耦合加速收斂方法的收斂歷史及其與GKUA的對比情況. 可以看出,本方法能大幅度加快收斂速度. 對于Kn= 0.01的案例, 迭代步數加速約7.5倍, 計算耗時提升約7.3 倍 (GKUA 耗時 116 h, Coupled 耗時 15.9 h);對于Kn= 0.1的案例, 迭代步數加速約2.85倍,計算耗時提升約2.7倍(GKUA耗時77 h, Coupled耗時 28.6 h). 另外, 數值實驗發現, 采用更稀疏的網格有利于提升加速收斂效果, 但是, 對壁面物理量計算結果有負面影響. 分析認為, 稀疏網格有利于宏觀方程的內迭代收斂, 并利于發揮宏觀方程在擾動傳遞方面的優勢, 促進整體收斂. 另一方面,相較于只存在一階空間導數的氣體動理論方法,宏觀方程存在的二階導數項導致數值求解宏觀方程在網格無關性方面的要求高于氣體動理論方法.

圖 5 Kn = 0.1 圓柱繞流流場對比 (a) 壓力; (b) 溫度;(c) 馬赫數 (GKUA: 彩色背景及白實線; Coupled: 紅色虛線)Fig. 5. (a) Pressure, (b) temperature, (c) Mach number distribution around cylinder for Kn = 0.1 (GKUA: coloured background and white solid lines; Coupled: red dash lines).

圖 6 圓柱壁面的 (a) 壓力, (b) 熱流, (c) 剪切應力Fig. 6. (a) Pressure, (b) heat flux, and (c) shear stress profile along the wall surface of cylinder.

圖 7 耦合加速收斂方法與常規GKUA的超聲速圓柱繞流計算收斂情況對比Fig. 7. Comparison of the convergence history of supersonic flow around the cylinder between the coupled acceleration method and the conventional GKUA.

4.3 并列雙圓柱干擾繞流

圖 8 并列雙圓柱多塊網格布局Fig. 8. The multi-blocks mesh layout for two side-by-side cylinders.

圖 9 并列雙圓柱繞流流場對比 (a) 壓力; (b) 溫度; (c) 馬赫數 (GKUA: 彩色背景及白實線; Coupled: 紅色虛線)Fig. 9. (a) Pressure, (b) temperature, (c) Mach number field for two side-by-side cylinders (GKUA: coloured background and white solid lines; Coupled: red dash lines).

圖 10 上圓柱壁面的 (a) 壓力, (b) 熱流, (c) 剪切應力Fig. 10. (a) Pressure, (b) heat flux, and (c) shear stress profile along the surface of upper cylinder.

為了進一步驗證耦合加速收斂方法對復雜流場的處理能力, 本文對并列雙圓柱干擾繞流進行了模擬. 相較于單體繞流, 并列雙圓柱干擾繞流存在著激波融合、干擾及邊界層相互影響. 為了與已有研究成果對比, 以圓柱半徑為參考長度, 本文計算了來流努森數Kn= 0.1, 馬赫數Ma= 2 的超聲速雙圓柱干擾繞流. 努森數Kn定義與上一案例相同. 無量綱壁面溫度等于來流靜溫, 即Tw=1.0 . 利用多塊對接網格技術, 并列雙圓柱物理空間網格布局如圖8所示. 兩個半徑為1的圓柱的圓心距離為 4, 上下布置于計算域中. 在 [–6, 6]× [–6, 6]的速度相空間里布置了 40 × 40 個 Gauss- Legendre離散速度點. 圖9繪出了壓力場、溫度場及馬赫數分布的耦合加速收斂方法計算結果及其與常規GKUA計算結果的對比情況. 結果顯示, 采用本文提出的方法計算獲得的流場與常規GKUA計算獲得的流場符合良好, 包括激波、尾流區及雙圓柱干擾區等特殊區域. 為了定量對比, 圖10給出了本文方法獲得的上圓柱壁面壓力、熱流及剪切應力結果及其與常規GKUA結果、文獻[47]給出的DSMC 結果的對比情況. 橫坐標q= –180°—180°對應于從圓柱迎風面頂點沿逆時針分布的壁面位置. 結果顯示, 本文提出的方法獲得的壁面壓力、剪切應力與常規GKUA獲得的結果以及DSMC結果符合良好. 在熱流方面, 本文方法和常規GKUA,DSMC結果整體符合良好, 在頭部略有差異. 由此可以看出, 本文方法能夠準確反映常規GKUA在描述跨流域復雜流動問題的特性. 圖11給出了耦合加速收斂方法的收斂歷史及其與GKUA的對比情況. 其迭代步數加速約9.3倍, 計算耗時提升約2.5 倍 (GKUA 耗時 42.7 h, Coupled 耗時 16.8 h).對比上一案例可知, 本案例離散速度空間大小縮小約6倍, 導致宏觀方程內迭代耗時的占比增大, 影響了計算耗時的優化效果.

圖 11 耦合加速收斂方法與常規GKUA的并列雙圓柱超聲速繞流計算收斂情況對比Fig. 11. Comparison of the convergence history of supersonic flow around two side-by-side cylinders between the coupled acceleration method and the conventional GKUA.

5 結 論

為了提高氣體動理論數值方法的計算效率, 針對近空間飛行環境的稀薄氣體非平衡流動特點, 本文利用宏觀方程在擾動傳播、演化方面的優勢, 構造了基于耦合宏觀方程本構關系的氣體動理論加速收斂方法, 以加速氣體動理論方法的收斂速度.通過理論建模及數值驗證得到如下結論.

1)通過對Shakhov模型方程求碰撞不變量的矩, 建立了與模型方程相容的宏觀方程及本構關系, 并通過Chapman-Enskog漸近展開將宏觀方程的應力、熱流的一階項與高階項剝離. 以氣體動理論方法提供的應力、熱流高階項數值解作為紐帶, 實現了宏觀方程的封閉.

2)在氣體動理論統一算法框架下, 改進了具有加速收斂能力的有限體積LU-SGS全隱格式. 在格式中, 宏觀方程的收斂解被用于模型方程的當地平衡態速度分布函數及碰撞頻率的更新、演化.

3)通過對不同流域方腔流動的數值模擬, 證明了本研究提出的理論方法能正確描述不同稀薄條件下的本構關系. 獲得的溫度場、速度場與常規GKUA及參考文獻的結果符合良好. 并且, 本方法對近連續流區低速流動問題的加速收斂作用明顯,如對Kn= 0.01的方腔流動, 本方法相較常規GKUA的加速比為47. 隨著稀薄程度增加, 氣體對流輸運效應減弱, 宏觀方程對迭代的加速效果降低.

4)通過對圓柱繞流、并列雙圓柱干擾繞流的數值模擬, 驗證了方法在處理激波、壁面強剪切、流動分離等物理特征的能力. 獲得的流場及壁面物理量與常規GKUA及DSMC結果符合良好. 并且, 實現了數倍的加速收斂效果. 類似于方腔流動,隨著稀薄效應的增加, 加速效果逐漸降低. 結合各算例的數值試驗經驗, 分析認為, 加速收斂方法在Kn小于0.5的案例中都能實現加速收斂效果.

5)如何加速宏觀方程內迭代收斂速度及優化迭代策略, 減少宏觀方程的迭代耗時, 是進一步提升本方法效率必須考慮的重要問題, 有待深入研究.

感謝中國科學院力學研究所李新亮研究員及團隊提供的有益支持.