氣體運動論統一算法在跨流域轉動非平衡效應模擬中的應用

蔣新宇， 李志輝， 吳俊林

氣體運動論統一算法在跨流域轉動非平衡效應模擬中的應用

蔣新宇1，2， 李志輝1，2， 吳俊林1

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

從考慮轉動松弛變化特性的Rykov模型出發，基于Boltzmann模型方程求解跨流域氣體運動論統一算法原理與計算規則，采用轉動慣量描述氣體分子自旋運動，研究考慮轉動非平衡影響的Boltzmann模型方程數值求解方法.通過對氮氣激波結構、二維鈍頭體和三維尖雙錐跨流域繞流的模擬分析，驗證該算法的跨流域一致適用性.

轉動非平衡；統一算法；激波結構；鈍頭體；尖雙錐

0 引言

飛行器往返大氣層經歷連續流、近連續流、滑移過渡流以至高稀薄自由分子流［1］等，不同流域氣體熱力學性質及繞流狀態互不相同.在稀薄過渡流區，由于連續介質假設失效，N-S方程已不能有效描述飛行器繞流流動現象；而基于微觀分子運動與碰撞隨機統計模擬的DSMC方法因受網格劃分、時間步長限制，難以對低Knudsen數高空近連續稀薄流進行數值仿真［2－3］.如何準確模擬近連續滑移過渡流區高超聲速氣動力／熱環境，一直是人們所關心的問題［4－5］.

氣體分子運動論的基本方程—Boltzmann方程本身可描述各個流域的氣體分子輸運現象，通過求解Boltzmann方程可以研究各流域氣體動力學問題.但該方程是一個高度非線性積分－微分方程，幾乎不可能求出精確解.為避免求解復雜碰撞項，眾多學者基于守恒不變量，由數學上較簡單的統計和碰撞松弛模型代替Boltzmann方程碰撞項，提出許多氣體運動論模型方程.最簡單的BGK模型方程于1954年提出，隨后，人們發展了多種形式的Boltzmann模型方程，諸如Holway的橢球模型［6］、Shakhov基于BGK方程的修正而得到的高階推廣模型方程等［7］.因此，求解簡化的Boltzmann模型方程，是一個較為經濟有效的解決途徑.李志輝［8－11］從研究建立描述各流域均適用的氣體分子速度分布函數方程出發，吸收計算數學指數型積分求解原理，提出氣體運動論離散速度坐標法，發展經改進的Gauss-Hermite積分方法等系列離散速度數值積分技術，消除原分布函數對速度空間的連續依賴性，將Boltzmann模型方程化為在各個離散速度坐標點處基于時間、位置空間具有非線性源項的雙曲型守恒方程，利用二階Runge-Kutta方法數值求解源項方程，引入NND格式數值求解對流運動項.基于對氣體分子與物面相互作用原理的研究，發展可用于氣體分子速度分布函數有限差分求解的氣體運動論邊界條件數學模型.由此建立求解一維、二維、三維Boltzmann模型方程的氣體運動論統一算法，開展跨流域氣體繞流數值模擬應用［12－13］.為了研究近空間飛行器跨流域非平衡繞流問題，近年來，基于對轉動自由度松弛變化［14］、氣體分子運動論Rykov模型［15－16］的研究，在氣體分子速度分布函數演化更新求解中考慮轉動自由度影響，把氣體分子轉動能作為分布函數的自變量，在統一算法理論與計算框架下，提出考慮轉動非平衡效應Boltzmann模型方程數值算法［17］.本文在此基礎上，針對高、低不同馬赫數（1.53≤Ms≤25）氮氣一維激波結構，不同Knudsen數（10－3≤Kn∞≤10）二維鈍頭體繞流和三維雙錐外形繞流算例，研究考慮轉動非平衡效應Boltzmann模型方程在跨流域繞流問題中的模擬應用，計算剖析自由分子流到連續流跨流域中間過渡帶非平衡流動規律與繞流現象.

1 控制方程與計算方法

1.1 求解分子速度分布函數方程數值格式

基于轉動松弛特性Rykov模型［15－16］，采用轉動慣量來描述氣體分子的自旋運動，利用分子總角動量守恒作為一個新的碰撞不變量.基于氣體分子速度分布函數f＝f（r，ξ，t，e），這里r、ξ分別是位置空間和速度空間的坐標、e為內能，在求解Boltzmann模型方程的統一算法理論框架［8－13］下，引入權因子1和e對速度分布函數所依賴的速度空間進行無窮積分，構建經約化速度分布函數處理［17－18］，考慮轉動非平衡影響的Boltzmann模型方程，其無量綱形式為

這里，T∞和T*均為有量綱的值，但B＝T∞／T*是無量綱量．對于氮氣N2而言，T*＝91.5 K，1／δ＝1.55．

速度分布函數f0、f1是依賴于時間t、位置空間r和速度空間ξ的函數，需將速度空間離散降維，去掉其對速度空間的連續依賴性.氣體運動論離散速度坐標法［8－10，18］通過一套與積分規則相一致的離散速度坐標點分布，用具有N個元素的函數簇代替原分布函數對積分變量的連續依賴性，把模型方程轉換為在各離散速度坐標點處彼此獨立的關于位置空間和時間的雙曲型守恒方程，宏觀流動參數則可通過相應于離散速度坐標法的積分求解規則確定.如經速度空間離散處理的速度分布函數方程在計算平面內矩陣形式為

采用文獻［9－10，18］所建立二階有限差分格式，將時間分裂數值計算方法應用于方程（4）式中位置空間的三個對流運動方程和碰撞松弛源項方程進行耦合數值求解，得到在每個離散速度坐標點處直接求解分子速度分布函數的氣體運動論數值格式

1.2 邊界條件數值處理方法

對于流場邊界條件數值處理，按照統一算法基于分子速度分布函數邊界處理數學模型［11，13，18］，對于無窮遠處來流條件，我們以無量綱的來流宏觀參數n∞＝1、U∞、V∞、W∞和T∞＝1所表征的平衡態分布函數確定

第二種，適應系數通過分析氮氣中熱傳導的實驗數據得到［15］

2 數值計算結果與討論

2.1 高、低不同馬赫數定常正激波內流動

為檢驗考慮轉動非平衡影響的Boltzmann模型方程統一算法對高、低不同馬赫數激波內流動問題模擬能力，計算氮氣中激波內流動.為便于比較，引入無量綱參考值［18］：

圖1分別繪出馬赫數Ms＝1.53、3.2、10、25激波結構無量綱密度ρ*分布，圖中GKUA表示統一算法結果，Exp.為來自文獻［20］實驗數據，GBE為文獻［14，21］中的廣義退化Botzmann方程解算器計算結果.可看出，對于高、低不同馬赫數1.53～25激波內流動密度分布統一算法結果均與實驗數據、典型文獻模擬值吻合較好.圖2繪出考慮轉動非平衡影響的Ms＝25強激波內流動密度與溫度變化輪廓線，圖中實線n表示數密度，虛線T為全局溫度，點劃線Tt表示平動溫度，雙點劃線Tr表示轉動溫度.可看出，激波內流動從波前開始，溫度先于密度激發，不同模態的溫度中平動溫度又先于轉動溫度激發.激發的不同時正好說明了激波流動是一個非平衡變化過程.圖中還顯示出整個非平衡流動區間中，平動溫度總是遠高于轉動溫度，這也表明平動溫度的主導地位，造成平動溫度輪廓線在臨近激波下游過渡區出現嚴重的非平衡凸起“超出”現象.圖1、圖2關于弱激波Ms＝1.53到強激波Ms＝25的計算，證實統一算法在計算分析雙原子氮氣高、低不同馬赫數激波內流動問題方面的準確可靠性.

圖1 高低不同馬赫數下激波內流動密度分布Fig.1 Density distributions in shock wave at various Mach numbers

2.2 二維鈍頭體全飛行流域繞流

為了研究從自由分子流到連續流跨流域二維繞流流場中的內能非平衡影響，擬定圖3所示的一個鈍頭體［21］外形，其中方形柱段長度是圓頭半徑的兩倍，計算來流馬赫數Ma＝3、克努森數Kn＝10、1、0.1、0.001的繞流流場.圖4繪出了三種不同Kn數下，鈍體繞流流場的密度分布等值線圖.可以看出：在Kn＝10時，流場擾動區域很大，沒有出現參數強間斷形式的激波，流場為稀薄流；在Kn＝0.001時，流場存在明顯的激波，擾動無法傳播至上游，流場為連續流、近連續流；而在Kn＝1時，可以看到物體前面厚厚的激波強擾動，流場體現為過渡區繞流.通過對自由分子流到連續流全流域的鈍頭體繞流數值模擬，證明了統一算法的跨流域一致適用性.

圖2 Ms＝25氮氣激波結構內流動密度分布與平動、轉動和全場溫度分布Fig.2 Distribution of density，translational temperature，rotational temperature and global temperature in nitrogen shock wave at Ms＝25

圖3 鈍頭體外形尺寸示意圖Fig.3 Blunt body geometries

圖4 不同Kn數下鈍體繞流的密度分布Fig.4 Density contours around a blunt body

為了剖析過渡區流動的非平衡效應，圖5、圖6分別繪出了Kn＝1與Kn＝0.001兩種狀態對應的繞流流場溫度分布，其中圖（a）表示平動溫度等值線、圖（b）表示轉動溫度等值線.可以看出：在Kn＝0.001對應的近連續流區，平動溫度與轉動溫度分布趨于一致，流場等值線結構基本無差別，流動可以認為是平衡流；而Kn＝1對應的高稀薄流區，平動溫度與轉動溫度彼此差別較大，計算得到的物體前部高溫區平動溫度最大值約為3，而轉動溫度最大僅1.6左右，兩者相差近一倍，流場表現出強烈的稀薄非平衡特性.

2.3 雙錐外形體跨流域繞流

為了檢驗考慮轉動非平衡影響的Boltzmann模型方程統一算法對三維飛行器跨流域繞流模擬能力，擬

圖5 鈍頭體繞流（Ma＝3，Kn＝1）流場溫度分布等值線Fig.5 Temperature contours around a blunt body at Ma＝3，Kn＝1

圖6 鈍頭體繞流Ma＝3，Kn＝0.001流場溫度分布等值線Fig.6 Temperature contours around a blunt body at Ma＝3，Kn＝0.001

定圖7所示雙錐外形.為了與文獻［22］使用GBE解算器計算高稀薄流區來流馬赫數Ma＝3、克努森數Kn＝1的上述尖頭雙錐外形繞流狀態分析比較，本文計算了該外形以Ma＝3通過稀薄過渡區Kn＝1、0.1的繞流流場.圖8繪出了統一算法計算上述尖雙錐外形繞流Ma＝3、Kn＝1駐點線密度分布與文獻［22］Wilson計算結果（符號Δ表示）的對比，可看出兩種方法有較好的一致，證實考慮轉動非平衡效應的統一算法用于計算三維飛行器跨流域繞流問題的可行性.圖9繪出Ma＝3、Kn＝1時的稀薄過渡流區雙錐繞流流場的

圖7 雙錐外形示意圖Fig.7 Bicone geometries

圖8 Ma＝3，Kn＝1駐點線密度分布Fig.8 Density distribution along the stagnation line at Ma＝3，Kn＝1

圖9（a）平動溫度和圖9（b）轉動溫度分布，可以看出兩者有較大的差距，這體現了高Kn數稀薄過渡流區存在嚴重的流動非平衡，說明了在稀薄流區考慮轉動平衡效應的必要性.圖10繪出Kn＝1和Kn＝0.1兩種狀態下尖錐繞流流場全局溫度分布.從圖中可以看出，Kn＝1時，流動呈現較強的稀薄流特征，激波強擾動不明顯，流動參數沒有強間斷，物體前方擾動區域范圍很大；而Kn＝0.1時，流動趨近物體前端出現激波強擾動現象.

圖9 Ma＝3，Kn＝1尖雙錐繞流溫度分布Fig.9 Temperature contours around bicone at Ma＝3，Kn＝1

圖10 不同努森數尖雙錐繞流全局溫度分布Fig.10 Temperature distributions around bicone at different Kn

我行還計算分析了該外形在Ma＝3、Kn＝0.000 1連續流區的繞流流場.如圖11所示.從密度、溫度、馬赫數等參數的等值線圖上均可以明顯看出斜激波和弓形激波的相互干擾.從熱流云圖上可以看出，脫體弓形激波強度遠強于第一個錐產生的附體斜激波，熱流峰值出現在弓形激波處和斜激波后的附面層內，說明飛行器在第一個錐面存在較大熱流分布，需要考慮熱防護，該圖直觀再現了雙錐外形連續流區繞流存在斜激波和弓形激波相互干擾的復雜流動結構，計算結果與理論分析相吻合.從馬赫數等值線圖上也可以看出該繞流的簡要流動特征：遠前方氣流繞流第一個錐時會產生一道較弱的斜激波，經斜激波壓縮后氣流速度降低，但仍然是超音速流動，氣流繼續向前繞流第二個錐時會產生一個脫體的弓形激波，兩道激波在兩個錐面結合處會聚相遇相互擠壓影響，產生一道更強的正激波，波后氣流變為亞音速.當氣流完全繞過錐體，在第二個錐面結束處因流動急劇膨脹而產生一束膨脹波系，重新將繞流氣體加速到超音速往后流去.

3 結論

基于對氣體分子轉動自由度和Rykov模型，應用稀薄流到連續流氣體運動論統一算法數值求解考慮轉動非平衡影響的Boltzmann模型方程.針對不同馬赫數定常正激波結構、二維鈍頭體全飛行流域繞流、尖雙錐外形體跨流域繞流問題計算分析，得到如下結論：

圖11 Ma＝3，Kn＝0.000 1尖雙錐繞流流場Fig.11 Flow around bicone at Ma＝3，Kn＝0.000 1

1）通過對寬廣馬赫數范圍內（1.53≤Ms≤25）考慮轉動非平衡影響的氮氣激波結構計算結果與實驗數據、典型文獻模擬值比較分析，驗證了考慮轉動非平衡影響的統一算法對求解強、弱激波一維流動問題的準確可靠性，并分析顯示了高超聲速情況下激波結構中的非平衡流動特征與變化規律.

2）通過對二維鈍頭體外形從自由分子流到連續流跨流域繞流流場的模擬研究，顯示了跨流域不同流動特征，驗證了考慮轉動非平衡影響的統一算法對全飛行流域模擬計算的一致適用性.通過對比分析繞流流場平動溫度與轉動溫度分布，證實來流Kn數越高，流動非平衡效應越明顯.

3）通過對尖雙錐外形跨流域繞流狀態的計算分析，描述了稀薄氣體轉動非平衡效應和連續流區的激波干擾現象.驗證了考慮轉動非平衡效應的統一算法對三維復雜飛行器跨流域繞流問題的模擬能力.

本文工作僅考慮轉動非平衡影響的Boltzmann模型方程統一算法在跨流域繞流問題方面的初步應用與計算檢驗，將其應用于跨流域復雜高超聲速氣動力／熱繞流問題以及考慮轉動能與振動能、混合氣體非平衡效應的適應性研究等，還有待進一步深入.

致謝：本文工作得到國家自然科學基金重大研究計劃‘近空間飛行器的關鍵基礎科學問題'重點支持項目（91016027）資助，大部分并行計算在國家超級計算天津中心和國家超級計算濟南中心完成，特此表示感謝.

［1］Qian X S.Superaerodynamics，mechanics of rarefied gases［J］.Journal of Aeronautics Science，1946，13（12）：653－664.

［2］Koppenwallner G，Legge H.Drag of bodies in rarefied hypersonic flow［C］∥Progress in astronautics and aeronautics：Thermophysical aspects of reentry flows，New York，AIAA Paper 85－0998，1985，103：44－59.

［3］李志輝，吳振宇.阿波羅指令艙稀薄氣體動力學特征的蒙特卡羅數值模擬 ［J］.空氣動力學學報，1996，14（2）：230－233.

［4］Whitehead A J.NASP aerodynamics［J］.AIAA Paper 89－5013，1989.

［5］莊逢甘，崔爾杰，張涵信.未來空間飛行器的某些發展和空氣動力學的任務 ［C］∥中國第一屆近代空氣動力學與氣動熱力學會議論文集，2006，8：1－12.

［6］Holway JLH.New statisticalmodels for kinetic theory，Methods of construction［J］.Phys Fluids，1966，9（9）：1658－1673.

［7］Shakov E M.Generalization of the Krook kinetic relaxation equation［J］.Fluid Dynamics，1968，3（1）：95－96.

［8］李志輝，張涵信.稀薄流到連續流的氣體運動論統一數值算法初步研究 ［J］.空氣動力學學報，2000，18（3）：251－259.

［9］李志輝，張涵信.稀薄流到連續流的氣體運動論模型方程算法研究［J］.力學學報，2000，34（2）：145－155.

［10］Li Z H，Zhang H X.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum［J］.J of Comput Phys，2004，193：708－738.

［11］Li Z H.Gas-kinetic unified algorithm for re-entering complex problems covering various flow regimes by solving Boltzmann model equation［M］∥Jason H.Advances in spacecraft technologies，InTech Publisher，2011：273－332.

［12］Li Z H，Zhang H X.Gas-kinetic numerical studies of three-dimensional complex flows on spacecraft re-entry［J］.JComput Phys，2009，228：1116－1138.

［13］李志輝，張涵信.基于Boltzmann模型方程各流域三維復雜繞流問題統一算法研究 ［J］.中國科學，2009，39（3）：414－427.

［14］Agrwal R K，Chen R，Cheremisin F G.Computation of hypersonic shock wave flows of a diatomic using the generalized Boltzmann equation［C］.AIAA 2007－4541，39th AIAA Thermophysics Conference，Miami，FL，June 25－28，2007.

［15］Rykov V A.Model kinetic equation of a gaswith rotational degrees of freedom［J］.Fluid Dynamics，1975，（10）：959－966.

［16］Rykov V A，Titarev V A，Shakhov E M.Shock wave structure in a diatomic gas based on a kinetic model［J］.Fluid Dynamics，2008，43（2）：316－326.

［17］李志輝，吳俊林，彭傲平，張涵信.考慮轉動非平衡效應Boltzmann模型方程統一算法與跨流域繞流問題模擬研究［C］∥第七屆全國流體力學學術會議論文集，CSTAM2012-B03-0195，11月12～14日，廣西桂林，2012.

［18］李志輝.從稀薄流到連續流的氣體運動論統一數值算法研究［D］.中國空氣動力研究與發展中心研究生部，2001.

［19］吳俊林，李志輝，蔣新宇.考慮轉動能影響的一維／二維Boltzmann-Rykov模型方程數值算法應用研究 ［J］.計算物理，2013，30（2）：329－339.

［20］Alsmeyer H.Density profiles in argon and nitrogen shock wavesmeasured by the absorption of an electron beam.［J］.JFluid Mech，1976，74：497－513.

［21］Chen R，Agrwal R K，Cheremisin F G.Computation of supersonic flow past2-D bodies using a Boltzmann solver［C］.AIAA 2005－1102.

［22］Wilson C D，Agarwal R K，Tcheremissine F G.Computation of hyersonic flow of a diatomic gas in rotationa nonequilibrium past3D blunt bodies using the generalized Boltzmann equation［C］.AIAA 2009－3836.

Application of Gas-K inetic Unified Algorithm Covering Various Flow Regimes for Rotational Non-equilibrium Effect

JIANG Xinyu1，2，LIZhihui1，2，WU Junlin1
（1.China Aerodynamics Research and Development Center，HAI，Mianyang 621000，China；2.National Laboratory for Computational Fluid Dynamics，Beijing 100191，China）

With gas-kinetic unified algorithm （GKUA）based on Boltzmann equation for various flows regimes，rotational nonequilibrium effect is investigated in Rykov model，in which spin movement of gas molecules is described with moment of inertia. Numericalmethod of Boltzmannmodel equation involving rotationalnon-equilibrium effect is developed.Nitrogen shock wave structure，flows around a 2D blunt-head body and flows around a 3D tine bicone are simulated to validate themethod in various flows regimes.

rotational non-equilibrium；GKUA；shock wave structure；blunt-head body；tine bicone

date： 2013－03－17；Revised date： 2013－07－17

V211.3

A

1001-246X（2014）04-0403-09

2013－03－17；

2013－07－17

國家自然科學基金（91016027）資助項目

蔣新宇（1987－），男，碩士，助理工程師，從事稀薄氣體動力學研究

李志輝，E-mail：zhli0097＠x263.net