李 凡，李志輝，*，李中華，羅萬清，陳愛國

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

0 引言

軸對稱流動廣泛存在于自然界，無論是在軌航天器還是微納米器件中，都有很強的工程應用背景。針對軸對稱流動Boltzmann模型方程數(shù)值模擬，如果采用全三維模擬，存在計算量巨大和所需內(nèi)存量十分龐大的難題。相比于求解三維方程，求解針對軸對稱流動的簡化準二維Boltzmann模型方程，在計算量和內(nèi)存量上少了一個維度，計算開銷大大降低[1-2]，并且在軸對稱特性得到滿足的同時，也避免了繁瑣的網(wǎng)格處理和全三維模擬網(wǎng)格生成過程中的奇異性問題[3]。因此相對于全三維模擬，發(fā)展一種基于準二維軸對稱Boltzmann模型方程計算模擬的工程應用方法是十分必要的。

近幾十年來，基于準二維方程來模擬軸對稱流動的算法取得了飛速發(fā)展[2-15]。其中一部分算法[4-7]是基于準二維軸對稱納維-斯托克斯（Navier-Stokes,N-S）方程發(fā)展而來的。由于N-S方程是在連續(xù)性介質假設基礎上建立的，其控制方程本身在稀薄流區(qū)與實際情況存在著至少是定量上的差距，因此該類算法無法完成全流域軸對稱流動的模擬。另一部分算法[3,8-15]是基于Boltzmann方程發(fā)展而來。Boltzmann方程的建立基于如下兩個重要假設[16]：（1）分子相互碰撞時只考慮二體碰撞，即認為3個分子或3個以上分子同時碰撞在一起的幾率很小；（2）分子混沌假設。實踐表明，不僅稀薄氣體流動滿足這些假設，而且處于海平面標準條件下的大氣也滿足Boltzmann條件，因此Boltzmann方程本身適用于全流域。為了模擬全流域軸對稱流動，文獻[17-21]通過研究確立描述各流域微觀分子輸運現(xiàn)象的Boltzmann模型方程，發(fā)展了氣體分子運動論離散速度坐標法對速度分布函數(shù)進行數(shù)值離散，建立了可用于高、低不同馬赫數(shù)宏觀流動取矩的離散速度數(shù)值積分法。將計算流體力學有限差分無波動無自由參數(shù)（nonoscillatory nonfree dissipative, NND）格式[22]推廣拓展到基于時間、位置空間和速度空間的Boltzmann模型方程數(shù)值求解，提出能可靠模擬從稀薄氣體自由分子流到連續(xù)流的氣體動理論統(tǒng)一算法[17-21]與大規(guī)模并行計算方法[23]，并開展了微機電系統(tǒng)（micro electro mechanical systems,MEMS）[19-20]、近地空間航天飛行器跨流域氣體流動問題應用研究[24-25]，對本文進一步研究直接求解軸對稱Boltzmann方程可計算建模提供了較強的研究基礎與可行性。

目前，最常用的稀薄流計算方法是直接模擬蒙特卡洛（direct simulation Monte Carlo, DSMC）方法[8-11]。DSMC通過跟蹤大量分子的自由運動和分子之間的相互碰撞，經(jīng)過統(tǒng)計平均處理得到氣體的宏觀物理量，如密度、速度、能量、熱流等，與直接求解Boltzmann方程是等價的[10]。由于DSMC方法本身涉及到統(tǒng)計平均過程，其中的統(tǒng)計誤差造成了在低速流動中的統(tǒng)計漲落問題。另一方面，DSMC方法的物理空間網(wǎng)格步長受到分子平均自由程的限制，時間步長受到分子平均碰撞時間的限制[11]。當模擬連續(xù)流區(qū)的軸對稱流動時，會存在計算效率低下、內(nèi)存量十分龐大的問題，以致計算困難。因此基于DSMC模擬全流域軸對稱流動的途徑仍不是合適的選擇。另一種方法是基于準二維Bhatnagar-Gross-Krook （BGK）類模型方程[3,12-15]。BGK類模型方程是基于Boltzmann方程的簡化模型方程。研究發(fā)現(xiàn)，發(fā)展基于準二維BGK類模型方程的算法相比于基于笛卡爾坐標的算法來說，難度更大，其困難之處在于需要處理準二維模型方程中的軸對稱源項[1]。軸對稱源項來源于模型方程對流項的坐標變換，其中還包含速度空間的偏導。并且，離散求解軸對稱源項對算法格式的精度、穩(wěn)定性以及守恒性方面有很大影響[1]。

由于軸對稱源項具有強非線性的特點，許多研究正在尋找簡化軸對稱源項的方法。Bergers[26]假設氣體分子速度分布函數(shù)在軸向上處于平衡態(tài)，將非線性的軸對稱源項簡化為線性的，使得求解更加容易。He等[27]將軸對稱源項中的氣體分子速度分布函數(shù)基于Chapmann-Enskog （CE）進行展開，構造了低速軸對稱流動的格子Boltzmann方法。Li等[28]基于CE展開，發(fā)展了求解連續(xù)流區(qū)軸對稱流動的氣體動理學格式（gas-kinetic scheme for Navier-Stokes equations,GKS-NS）。但是這些通過簡化軸對稱源項而構造出的算法，由于其特殊性，不能夠適用于全流域模擬。只有基于原始的軸對稱源項進行離散求解，才能對全流域軸對稱流動進行模擬。根據(jù)這一思路，李詩一等[3]基于原始笛卡爾坐標下的時間演化解構造出了柱坐標下的軸對稱源項的時間演化解，提出了基于有限體積法的全流域軸對稱流動統(tǒng)一氣體動理學格式（unified gas-kinetic scheme for axisymmetric flow,UGKS-AS）。

研究發(fā)現(xiàn)，若能適當修正BGK模型碰撞松弛參數(shù)和當?shù)仄胶鈶B(tài)分布，經(jīng)修正的BGK方程將能用于描述遠離平衡態(tài)以至全流域復雜的氣體流動問題[29]?；贐oltzmann方程可計算建模，李志輝等建立了一套全流域下的氣體動理論統(tǒng)一算法[17-21,24-25]（gaskinetic unified algorithm, GKUA）。本文將基于氣體動理論統(tǒng)一算法GKUA，采用新方法直接求解軸對稱源項，并離散求解準二維Boltzmann模型方程，發(fā)展基于有限差分法的全流域軸對稱流動氣體動理論統(tǒng)一算法，并用于空天飛行環(huán)境噴管流動模擬[30-32]。

1 氣體動理論統(tǒng)一算法Boltzmann模型方程

1.1 軸對稱Boltzmann模型方程推導

將宏觀流體力學與微觀分子動力學連接起來的介觀Boltzmann速度分布函數(shù)方程，可描述氣體分子速度分布函數(shù)基于位置空間、速度空間任一時刻由非平衡態(tài)向平衡態(tài)的演化：

其中，f、f1和f′、f′1分別為兩個分子碰撞前后的分子速度分布函數(shù)，V、V1為碰撞前的分子速度，g為氣體分子相對運動速率，b為碰撞因子，?為方位角。

氣體動理論統(tǒng)一算法將該方程化為各流域不同尺度間分子輸運現(xiàn)象統(tǒng)一模型方程，其在三維笛卡爾坐標系下的形式如下所示：

該方程通過描述氣體流動過程中分子速度分布函數(shù)f對位置空間、速度空間和時間的變化率關系，基于氣體分子碰撞松弛參數(shù)v和當?shù)仄胶鈶B(tài)速度分布函數(shù)fN，將不同流域流態(tài)控制參數(shù)、宏觀流動物理量、氣體黏性輸運特性、熱力學效應及氣體分子相互作用規(guī)則與分子模型用統(tǒng)一表達式聯(lián)系起來，由非平衡態(tài)向平衡態(tài)的演化，將各個流域不同尺度間分子輸運現(xiàn)象統(tǒng)一[17-21,24-25]。本文以此為基礎，定義柱坐標系下的徑向速度Vr和周向速度Vφ為：

經(jīng)過坐標變換后，得到柱坐標系下的準二維Boltzmann模型方程為：

根據(jù)文獻[1]，定義：

從而得到柱坐標系下的守恒型準二維Boltzmann模型方程為：

其中，

其中，x、r為柱坐標系下的空間位置，Vx、ζ、ω為分子速度，fM為麥克斯韋平衡態(tài)速度分布函數(shù)，p為氣體壓力，T為氣體溫度，Pr為普朗特數(shù)，n為氣體分子數(shù)密度，q為熱流，U、V、W為氣體宏觀流動速度，n∞為遠前方未擾來流氣體分子數(shù)密度，T∞為未擾來流溫度，R為氣體常數(shù)，χ為黏性系數(shù)的溫度依賴冪指數(shù)，λ∞為未擾來流氣體分子平均自由程。

1.2 軸對稱源項數(shù)值離散

使用三角算子來離散軸對稱源項：

該三角算子具有二階精度，同時滿足穩(wěn)定性，可用于粗網(wǎng)格，也可基于顯式Mac Cormack算子分裂思想，對軸對稱源項進行前差預測、后差校正離散：

1.3 Boltzmann模型方程有限差分數(shù)值格式求解

采用守恒型離散速度坐標法消除氣體分子速度分布函數(shù)對積分變量的連續(xù)依賴性，基于氣體分子速度分布函數(shù)方程對流運動和碰撞松弛的非定常特性，采用非定常時間分裂數(shù)值方法，將離散后的模型方程分裂為代表碰撞松弛變化的非線性源項方程、位置空間對流運動方程以及軸對稱源項方程，進行耦合迭代求解：

采用解析方法差分求解方程式（12），建立差分格式Ui+1=Ls(Δt)Ui，解除了數(shù)值格式穩(wěn)定性條件下分子碰撞時間對時間步長的限制；將所得結果作為初值帶入方程（13），采用NND-4（a）預測、校正兩步格式求解方程（13），建立差分格式Ui+1=Lr(Δt)Ui；將所得結果作為初值代入方程（14），采用NND-4（a）預測、校正兩步格式求解方程（14），建立差分格式Ui+1=Lx(Δt)Ui；將所得結果作為初值代入方程（15），采用二階龍格庫塔方法差分離散方程（15）的時間導數(shù)項，采用三角算子或者顯式Mac Cormack算子分裂差分離散軸對稱源項，進而求解方程（15），建立差分格式Ui+1=Lω(Δt)Ui。由此構造了求解軸對稱Boltzmann模型方程的有限差分數(shù)值格式如下：

考慮到實際流動中氣體分子速度分布函數(shù)方程的對流運動與碰撞松弛同時進行，為了避免差分離散引起計算順序的差別，將公式（16）前、后時間步交替耦合計算，得到求解跨流域統(tǒng)一軸對稱Boltzmann模型方程的有限差分格式為：

2 算法驗證與分析

2.1 同軸圓筒間的定常旋轉流動

考慮同軸圓筒間的氣體流動，外筒保持靜止，內(nèi)筒旋轉，流動僅與圓筒半徑r相關，可模型化為一維軸對稱模型。

圖1給出了圓筒的示意圖，外筒半徑R2=2 ，內(nèi)筒半徑R1=1。內(nèi)外筒壁面溫度Tw=1，內(nèi)筒旋轉速度，無量綱氣體常數(shù)R=0.5，氣體平均密度為ρav=1，壁面采用漫反射邊界，模擬氣體為單原子氣體。

圖1 圓筒間旋轉流動示意圖[1]Fig. 1 Schematic of rotational flow between two cylinders[1]

對于近連續(xù)過渡流0.02≤Kn≤10.0時，內(nèi)筒旋轉速度，氣體動理論統(tǒng)一算法GKUA計算所得無量綱密度和速度與統(tǒng)一氣體動理學格式（unified gas-kinetic scheme, UGKS）[3]結果對比如圖2所示，可以看出二者計算結果吻合較好。

圖2 同軸圓筒間定常旋轉流動GKUA與UGKS方法結果比較Fig. 2 Comparison of results between GKUA and UGKS methods for steady rotational flow between two coaxial cylinders

圖3給出了GKUA在連續(xù)流區(qū)和自由分子流區(qū)與其他方法計算的周向速度分布結果比較，其中當處于連續(xù)流區(qū)時，Kn=0.0001，uφin=0.1，處于自由分子流區(qū)時，Kn=100.0，uφin/=0.5。從圖3中可以看出，在連續(xù)流區(qū)時，GKUA所得的周向速度分布與不可壓的N-S解析解吻合較好；在自由分子流區(qū)時，GKUA所得結果與自由流解吻合較好。從圖2、圖3可以看出氣體動理論統(tǒng)一算法GKUA可以準確可靠模擬全局克努森數(shù)0.0001≤Kn≤100各流域的軸對稱旋轉Couette流動。

圖3 同軸圓筒間定常旋轉流動GKUA與N-S、自由流方法結果比較Fig. 3 Comparison of results between GKUA method and other methods for steady rotational flow between two coaxial cylinders

2.2 同軸圓筒間的非定常旋轉流動

在同軸圓筒旋轉速度突然改變的非定常旋轉流動模擬中，內(nèi)筒靜止，外筒旋轉。當t≤50時，外筒轉速uφout/=1.0；當t＞50時，uφout/=-1.0。圖4給出了Kn=0.02下不同時刻GKUA計算所得無量綱速度、溫度、壓力與DSMC結果[33]對比，二者吻合一致。

圖4 同軸圓筒間非定常旋轉流動GKUA與DSMC方法結果比較Fig. 4 Comparison of GKUA and DSMC methods for unsteady rotational flow between two coaxial cylinders

由于GKUA的求解過程是通過長時間的非定常模擬過渡到最終穩(wěn)態(tài)，因此可以捕捉到同軸圓筒間氣體流動的整個非定常演化過程，這也是區(qū)別于傳統(tǒng)離散速度法 （discrete velocity method, DVM）或離散坐標法（discrete ordinate method, DOM）的獨特優(yōu)勢。

2.3 軸對稱噴管道內(nèi)流動

選用型面如圖5所示的噴管進行軸對稱噴管內(nèi)流動模擬。噴管喉道處直徑5.1 mm，圓角半徑1.3 mm，入口截面距喉道30 mm，出口截面距喉道50.7 mm，喉道前后壁面與軸線夾角分別為30°和20°。噴管入口壓力pin=474Pa，溫度Tw=300K，外部壓力pout=1.529Pa，溫度Tout=203.03K。

圖5 軸對稱噴管幾何尺寸[34]Fig. 5 Geometrical size of the nozzle[34]

圖6給出了軸對稱噴管內(nèi)流動物理量沿噴管軸線分布的GKUA和N-S方法結果的比較。橫軸為軸線，取喉道截面處為坐標原點，從圖中可以看出，兩種方法計算結果變化趨勢相符合，這也證實了本文GKUA方法應用于噴管內(nèi)流動連續(xù)介質流動精細求解的正確性。

圖6 軸對稱噴管內(nèi)流動物理量沿噴管軸線分布的GKUA和N-S方法結果比較Fig. 6 Comparison of the results of GKUA and N-S methods for the profiles of flow physical quantity in axisymmetric nozzle along the nozzle axis

圖7給出了噴管軸線處GKUA方法計算的密度分布結果與Rothe[34]實驗值的比較。橫軸為距喉道的距離，Rt為喉道半徑。從圖中可看出模擬值與實驗值符合較好。通過與實驗數(shù)據(jù)進行對比，進一步確認了本文算法的可靠性。

圖7 噴管軸線密度分布模擬值與實驗值比較Fig. 7 Comparison of the simulated and experimental values of the density profile of the nozzle axis

2.4 噴管核心區(qū)羽流環(huán)境算法驗證與分析

為了準確分析空天飛行環(huán)境噴流/羽流特征，進行了低密度風洞羽流實驗與軸對稱Boltzmann模型方程統(tǒng)一算法及DSMC方法的驗證結合。高真空機組提供擴壓器實驗段長時間穩(wěn)定的真空環(huán)境，軸對稱姿控發(fā)動機噴管中的氣流在真空環(huán)境中急劇膨脹形成羽流。圖8、圖9分別給出了推力2 N發(fā)動機和推力10 N發(fā)動機羽流流場結構。發(fā)動機噴管根據(jù)姿控發(fā)動機的原型進行設計加工，其中2 N發(fā)動機噴管的收縮段和擴張段都是錐形型面，喉道直徑為1.2 mm，出口直徑為8.5 mm；10 N發(fā)動機噴管的擴張段型面是一圓弧，喉道直徑為3 mm，出口直徑為29.4 mm。圖9（a）為輝光放電顯示的流場，圖9（b）為數(shù)值模擬得到的流場，可看出幾種方法得到的流場結構一致。對于2 N發(fā)動機噴管，氣流沿錐形型面膨脹，邊界層厚度隨著膨脹逐漸增加。由于10 N發(fā)動機噴管擴張段的型面為一段圓弧，對流動有一定的壓縮作用，因而在噴管內(nèi)形成一道激波。激波沿弧形壁面發(fā)展，與邊界層相互干擾，形成一個復雜的波系結構，膨脹邊界范圍較大。比較分析圖8、圖9流動結構表明，推力2 N和10 N所形成的噴流/羽流流動狀態(tài)有較大差別，2 N推力較小，所形成的羽流呈完全膨脹擴散羽流狀；而10 N推力形成的噴管內(nèi)流動在擴張段因邊界層壓縮，形成上下壁面附近兩個壓縮激波從噴管出口噴出，在離噴口一定位置相交衍生出兩道膨脹波系，往后擴張形成馬赫盤魚尾狀羽流往后膨脹擴散。

圖8 推力2 N發(fā)動機羽流流場結構Fig. 8 Plume flow structure in the engine with 2 N thrust

圖9 推力10 N發(fā)動機羽流流場結構Fig. 9 Plume flow structure in the engine with 10 N thrust

圖10給出了推力10 N發(fā)動機羽流軸線上皮托管壓力分布比較。模擬計算的氣體介質為氮氣，計算狀態(tài)為：總壓p0=8.0×105Pa，總溫T0=773 K。從圖中可以看出數(shù)值模擬結果與風洞實驗結果比較接近，壓力分布符合較好。

圖10 在羽流中沿中心軸的皮托管壓力比較Fig. 10 Comparison of the pitot pressure along the central axis in the plume

3 結論

通過數(shù)值模擬同軸圓筒間的定常/非定常旋轉流動和軸對稱噴管內(nèi)流動，驗證確認了本文軸對稱Boltzmann模型方程統(tǒng)一算法處理全局克努森表征的全流域軸對稱噴管流動的準確可靠性，同時相對于傳統(tǒng)的離散速度坐標法，本文算法可靠捕捉軸對稱流動的非定常演化過程。

通過與低密度風洞實驗對比，計算得到的噴管出口核心區(qū)羽流結構一致、羽流軸線壓力分布一致。對于在軌航天器，本文算法可用于一體化計算其發(fā)動機內(nèi)部區(qū)域以及出口羽流軸對稱核心區(qū)域，相比于全三維模擬，其計算效率將得到很大提高。未來將深入研究考慮內(nèi)能激發(fā)多原子氣體以及混合氣體的軸對稱模型算法，并將該算法與DSMC方法和低密度風洞羽流實驗驗證相結合，3種方法互為補充完善，預期可發(fā)展揭示空天飛行環(huán)境噴流/羽流流場演變機理與流動特征的數(shù)值與實驗綜合分析模擬平臺。