基于宏觀方程數(shù)值本構(gòu)關(guān)系的氣體動(dòng)理論加速收斂方法*

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

1) (中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速所, 綿陽(yáng) 621000)

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

3) (國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室, 北京 100191)

1 引 言

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

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

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

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

2 理論模型

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

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

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

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

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

式中I為物理空間網(wǎng)格單元索引, 符號(hào)“?”表示物理量在網(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格式, 也可采用考慮氣體分子碰撞遷移物理過(guò)程, 具有漸近保持特性的重構(gòu)方式, 如DUGKS的重構(gòu)方式[40,41]. 為了簡(jiǎn)化問(wèn)題, 本研究采用經(jīng)典的CFD重構(gòu)格式進(jìn)行界面重構(gòu). 界面數(shù)值通量的計(jì)算采用具有迎風(fēng)特性的Steger-Warming矢通量分裂方法:

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

在等式兩邊同時(shí)加上n時(shí)刻殘差:

的控制方程:

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

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

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

CV為定容熱容.

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

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

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

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

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

并且, 高階項(xiàng) H oTσ及 H oTq獲取方法如下:

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

3 數(shù)值方法

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

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

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

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

1)確定離散速度空間及數(shù)值積分方法. 為保證離散速度空間能體現(xiàn)整個(gè)流場(chǎng)的速度分布形態(tài),對(duì)速度空間進(jìn)行截?cái)嗤ǔＰ枰獫M足[43]

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

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

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

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

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

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

4 算例驗(yàn)證

4.1 方腔流動(dòng)

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

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

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

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

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

4.2 超聲速圓柱繞流

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

圖 4 Kn = 0.01 圓柱繞流流場(chǎng)對(duì)比 (a) 壓力; (b) 溫度;(c) 馬赫數(shù) (GKUA: 彩色背景及白實(shí)線; 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).

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

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

圖 11 耦合加速收斂方法與常規(guī)GKUA的并列雙圓柱超聲速繞流計(jì)算收斂情況對(duì)比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é) 論

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

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

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

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

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

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

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