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

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

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

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

1 引 言

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

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

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

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

2 理論模型

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

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

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

采用離散速度坐標(biāo)法對約化速度分布函數(shù)方程進(jìn)行速度空間數(shù)值離散, 可得

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

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

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

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

的控制方程:

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

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

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

CV為定容熱容.

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

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

對Shakhov模型方程進(jìn)行Chapman-Enskog漸近分析可得:

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

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

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

采用迭代求解宏觀方程(10)獲得的宏觀量對方程(6)中的當(dāng)?shù)仄胶鈶B(tài)速度分布函數(shù)及碰撞頻率進(jìn)行預(yù)測, 并重新構(gòu)造隱式格式, 可將方程(9)改進(jìn)為如下形式:

3 數(shù)值方法

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

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

為了簡化, 省去無量綱量的上標(biāo)“?”. 若無特別說明, 后續(xù)量都為無量綱量.

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

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

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

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

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

5) 采用第4步獲得的宏觀流動量Wk+1?演化更新當(dāng)?shù)仄胶鈶B(tài)速度分布函數(shù)及碰撞頻率, 并采用LU-SGS方法求解構(gòu)造的隱式方程(20), 獲得第k+1步的離散速度分布函數(shù)fk+1及宏觀流動物理量Wk+1.

6)判斷是否達(dá)到收斂標(biāo)準(zhǔn). 本文采用全流場υk與υk+1的相對偏差均值作為收斂判斷標(biāo)準(zhǔn):

式中,I為物理空間網(wǎng)格索引, N umcell為物理空間總網(wǎng)格數(shù), ?t為時間步長. 若還未達(dá)到收斂標(biāo)準(zhǔn),繼續(xù)進(jìn)行第2—5步迭代.

4 算例驗證

4.1 方腔流動

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

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

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

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

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

4.2 超聲速圓柱繞流

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

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

圖 5 Kn = 0.1 圓柱繞流流場對比 (a) 壓力; (b) 溫度;(c) 馬赫數(shù) (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) 剪切應(yīng)力Fig. 6. (a) Pressure, (b) heat flux, and (c) shear stress profile along the wall surface of cylinder.

圖 7 耦合加速收斂方法與常規(guī)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 并列雙圓柱多塊網(wǎng)格布局Fig. 8. The multi-blocks mesh layout for two side-by-side cylinders.

圖 9 并列雙圓柱繞流流場對比 (a) 壓力; (b) 溫度; (c) 馬赫數(shù) (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) 剪切應(yīng)力Fig. 10. (a) Pressure, (b) heat flux, and (c) shear stress profile along the surface of upper cylinder.

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

圖 11 耦合加速收斂方法與常規(guī)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 結(jié) 論

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

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

2)在氣體動理論統(tǒng)一算法框架下, 改進(jìn)了具有加速收斂能力的有限體積LU-SGS全隱格式. 在格式中, 宏觀方程的收斂解被用于模型方程的當(dāng)?shù)仄胶鈶B(tài)速度分布函數(shù)及碰撞頻率的更新、演化.

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

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

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

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