稀薄氣體中高超速圓柱繞流仿真研究

黃 海 龍

(1.哈爾濱工業(yè)大學(xué) 土木學(xué)院，哈爾濱 150090;2.北京宇航系統(tǒng)工程研究所，北京 100076)

隨著飛行高度逐漸增加，當(dāng)飛行器高度達(dá)到或超過70 km 時，氣體的稀薄效應(yīng)開始明顯.該區(qū)域內(nèi)，Kn較大，氣體處于非連續(xù)區(qū)，此時，流體連續(xù)介質(zhì)假設(shè)失效，傳統(tǒng)的連續(xù)流體運(yùn)動方程(NS 方程)和相應(yīng)的CFD 方法已經(jīng)失效.針對這種稀薄氣體效應(yīng)，當(dāng)采用適當(dāng)?shù)难芯糠椒?，例如采用分子動力學(xué)模型和相應(yīng)的數(shù)值方法，其中直接模擬蒙特卡洛方法(DSMC)等方法就是該領(lǐng)域內(nèi)典型的研究方法.

直接模擬Monte carlo(DSMC)方法是依賴物理的概率模擬方法，在求解過渡領(lǐng)域流動的眾多解析、數(shù)值和模擬方法中，目前只有DSMC 方法可以模擬三維復(fù)雜真實稀薄氣體流動［1－6］.直接模擬蒙特卡洛方法可以模擬飛行器穿越高度為70～128 km 范圍內(nèi)的過渡區(qū)以及飛行器在70 km 以上高度巡航時的流場稀薄效應(yīng)，其所涉及的可研究領(lǐng)域包括高超聲速巡航飛行器、空間站、衛(wèi)星、返回艙等多種飛行器.

本文將用連續(xù)流體運(yùn)動方程(NS 方程)和DSMC 方法對四種典型Kn下的圓柱繞流進(jìn)行計算，分析比較不同Kn時，流場中的間斷結(jié)構(gòu)變化情況以及圓柱體后方的分離區(qū)變化情況.其中，Kn=0.001 時，流動環(huán)境處于連續(xù)介質(zhì)區(qū)，可采用連續(xù)流體運(yùn)動方程(NS 方程)求解，其余三種繞流流動采用DSMC 方法研究.

1 DSMC 方法概述

Bird［2］提出的DSMC(Direct Simulation Monte－Carlo)方法是依賴物理的概率模擬方法，它來源于分子動力學(xué)方法，采用幾率論方法判斷分子間是否發(fā)生碰撞.DSMC 方法從微觀角度出發(fā)，利用較少量的模擬分子，代表真實流體的大量分子，用計算機(jī)模擬由于氣體分子運(yùn)動、碰撞而引起動量和能量的輸運(yùn)、交換，產(chǎn)生氣動力和氣動熱這一宏觀物理過程.在求解過渡領(lǐng)域流動的眾多解析、數(shù)值和模擬方法中，DSMC 方法是惟一取得巨大成功的方法，從宏觀參量到細(xì)觀速度分布函數(shù)的水平上，該方法均能得到實驗的支持［3－7］.

DSMC 的模擬過程主要包括兩個部分:分子的微觀運(yùn)動和分子碰撞.一般，在實際情況中，分子的運(yùn)動和碰撞是耦合在一起的，而數(shù)值模擬方法中，很難將二者解耦，這給計算過程造成了一定困難.因此，在使用DSMC 方法對稀薄氣體流動進(jìn)行模擬計算時，必須對分子運(yùn)動和分子碰撞加以特殊對待.在微觀領(lǐng)域下，考慮限制DSMC 時間步長小于平均碰撞時間，這時，在平均意義上，單個分子在一個時間步長內(nèi)最多只有一次碰撞的機(jī)會;同時，再要求網(wǎng)格尺寸小于分子平均自由程，則碰撞只會在同一網(wǎng)格中的分子間發(fā)生.這樣就可在數(shù)值計算中，將分子的運(yùn)動和碰撞實現(xiàn)解耦.在模擬初始，根據(jù)物理條件，在計算域內(nèi)布置模擬分子并給定其位置和速度，在每一個時間步內(nèi)，分別計算運(yùn)動和碰撞.處理運(yùn)動時還包括分子進(jìn)入或逸出計算域(對應(yīng)于進(jìn)出口邊界條件)、與壁面的相互作用;碰撞的處理主要是碰撞對的選擇、碰撞截面積的計算.流動的宏觀量則由網(wǎng)格內(nèi)模擬分子的運(yùn)動參數(shù)統(tǒng)計平均得到.

2 計算條件和軟件設(shè)置

2.1 計算條件

本文計算了來流馬赫數(shù)為5，以圓柱半徑為特征長度的Kn分別為0.001、0.01、0.1、1 四種條件下的圓柱繞流.

圖1 流動分區(qū)示意圖

圖1 是流動分區(qū)示意圖.隨著Kn的增加，氣體流動可分為連續(xù)流、滑移流、過渡流和自由分子流四個流區(qū)，許多情況下，將滑移流區(qū)和過渡流區(qū)統(tǒng)稱為過渡區(qū)，除了連續(xù)流區(qū)外的其他區(qū)域統(tǒng)稱為稀薄流區(qū).從圖1 可以看出，Kn=0.001 處于連續(xù)流區(qū)和滑移流區(qū)的邊界，這時氣體比較稠密，連續(xù)介質(zhì)假設(shè)仍然成立，因此，此區(qū)域內(nèi)的流體流動常用NS 方程進(jìn)行計算.而如果用DSMC 計算的話，計算量太大，對CPU和內(nèi)存的需求都很大，計算時間也很長，計算結(jié)果也與求解NS 方程的結(jié)果一致，因此，本文采用連續(xù)流方程(NS 方程)計算Kn=0.001 時的圓柱繞流流場，可利用商軟Fastran 進(jìn)行仿真.

對于Kn為0.01，0.1、1 的三種情況，本文將采用基于DSMC 方法的程序進(jìn)行計算.下面將以Kn=0.01為例進(jìn)行詳細(xì)介紹.

下面的算例中，取圓柱取直徑D=100 mm，來流溫度，壁面溫度.

根據(jù)Kn的定義有:

其中:λ為分子運(yùn)動的平均自由程;L為流動特征長度，本文中取為圓柱半徑D/2.

由可以求得:

DSMC 算法要求，劃分網(wǎng)格時，每個網(wǎng)格尺度Δx、Δy、Δz均要小于當(dāng)?shù)氐姆肿悠骄杂沙苔?的1/3，實際計算中可以放寬到小于λ.

分子運(yùn)動的平均自由程還有如下關(guān)系:

其中:Mref為分子參考質(zhì)量，對空氣而言取Mref=48.1×10－27kg;Dref為分子參考直徑，流動介質(zhì)為空氣時，取Dref=4.19×10－10m;ρ為空氣密度.

再結(jié)合氣體分子數(shù)密度n 與氣體密度ρ 的關(guān)系式:

可以得到對應(yīng)Kn=0.01 下的氣體分子數(shù)密度n:

分子熱運(yùn)動的平均速度c為:

其中，k=1.3806×10－23J/K為波爾茲曼常數(shù);T為當(dāng)?shù)販囟?代入來流溫度T∞，并結(jié)合分子平均碰撞時間的定義，可以得到:

DSMC 算法要求，時間步長要小于分子平均碰撞時間Dt 的1/3.

來流溫度選為T∞=300K，來流馬赫數(shù)Ma∞=5，可以得到來流速度為:

2.2 網(wǎng)格劃分

本文使用CGNS 格式的結(jié)構(gòu)網(wǎng)格.采用單區(qū)或者多區(qū)的CGNS 格式的結(jié)構(gòu)網(wǎng)格，同時實現(xiàn)分區(qū)并行計算.在計算前，由于流場中分子平均自由程λ的分布未知，因此，需要對流場進(jìn)行預(yù)估，比如在來流部分滿足小于來流分子平均自由程λ 的1/3，在激波后，網(wǎng)格要適度加密，在尾流區(qū)網(wǎng)格可以適度放寬.

下面給出Kn=0.01 的算例所用的網(wǎng)格在圓柱壁面附近的網(wǎng)格大約尺寸見表1.

表1 網(wǎng)格尺寸(Kn=0.01)

該網(wǎng)格的計算域體積Vcell約為5.5×10－5m3，總網(wǎng)格數(shù)Ncell為540 000，因此，每個模擬分子代表的真實分子數(shù)Nn為:

在計算過程中，可以根據(jù)具體條件對Nn進(jìn)行適當(dāng)?shù)恼{(diào)整，比如，當(dāng)計算條件不足時，需要增大Nn，流場中的模擬分子數(shù)就會減少，但此時的計算精度會下降，因此，對Nn的調(diào)整需謹(jǐn)慎.

該網(wǎng)格分了四個Block，可以四個線程并行計算.利用MPI 進(jìn)行分區(qū)并行計算，為了提高計算效率，可以對計算域進(jìn)行合理分區(qū)，使每個區(qū)的計算量大致相同，這樣才不會因為某個區(qū)內(nèi)的計算量過大，導(dǎo)致整體計算速度降低，見圖2.

圖2 網(wǎng)格

3 結(jié)果分析

下面給出Kn分別為0.001、0.01、0.1、1 四種條件下的圓柱繞流的計算結(jié)果.其中連續(xù)介質(zhì)區(qū)(Kn=0.001)的工況由Fastran 完成，非連續(xù)介質(zhì)區(qū)(Kn=0.01、Kn=0.1、Kn=1)的三個工況由DSMC方法計算得到.

3.1 Kn=0.001

圖3 給出了Kn=0.001 工況下的馬赫數(shù)云圖和等值線.從圖3 可以明顯的看到，圓柱前方有清晰的脫體弓形激波，圓柱后緣遇到內(nèi)折再次產(chǎn)生激波，圓柱正后方有一條狹長的條狀低速區(qū)，低速區(qū)兩側(cè)有兩條滑移線.圖4 給出了Kn=0.001 工況下的壓力云圖和等值線，弓形激波后有明顯的高壓區(qū).圖5 是圓柱附近的流線圖，可以看出，圓柱后方出現(xiàn)了分離，有很大的漩渦區(qū).

圖3 馬赫數(shù)云圖(Kn=0.001)

圖4 壓力云圖(Kn=0.001)

圖5 圓柱附近流線(Kn=0.001)

3.2 Kn=0.01

圖6 給出了Kn=0.01 工況下的馬赫數(shù)云圖和等值線.從圖6 中還可以明顯的看到圓柱前方有清晰的脫體弓形激波，圓柱后緣遇到內(nèi)折再次產(chǎn)生激波，圓柱正后方有一條狹長的條狀低速區(qū)，低速區(qū)兩側(cè)有兩條滑移線.但相較Kn=0.001 的結(jié)果，激波和滑移線都明顯變粗，弓形激波的位置前移.圖7 給出了Kn=0.01 工況下的壓力云圖和等值線，弓形激波后的高壓區(qū)較Kn=0.001 的結(jié)果減弱.圖8 是圓柱附近的流線圖，圓柱后方仍存在分離區(qū)，但分離區(qū)范圍減小.

圖6 馬赫數(shù)云圖(Kn=0.01)

圖7 壓力云圖(Kn=0.01)

圖8 圓柱附近流線(Kn=0.01)

3.3 Kn=0.1

圖9 給出了Kn=0.1 工況下的馬赫數(shù)云圖和等值線.從圖9 中可以看出，弓形的脫體激波變得更厚，圓柱后緣的激波已經(jīng)不明顯，低速區(qū)兩側(cè)的兩條滑移線已經(jīng)消失.圖10 給出了Kn=0.1 工況下的壓力云圖和等值線，弓形激波后的高壓區(qū)進(jìn)一步減弱.圖11 是圓柱附近的流線圖，圓柱后已沒有分離區(qū)，分離區(qū)消失.

圖9 馬赫數(shù)云圖(Kn=0.1)

圖10 壓力云圖(Kn=0.1)

圖11 圓柱附近流線(Kn=0.1)

圖12 馬赫數(shù)云圖(Kn=1)

3.4 Kn=1

圖12 給出了Kn=1 工況下的馬赫數(shù)云圖和等值線.從圖12 中已經(jīng)看不出明顯的激波間斷，弓形激波變得很厚.圖13 給出了Kn=1 工況下的壓力云圖和等值線，弓形激波后的高壓區(qū)進(jìn)一步減弱.圖14 是圓柱附近的流線圖，圓柱后也已沒有分離區(qū)，分離區(qū)消失.

圖13 壓力云圖(Kn=1)

圖14 圓柱附近流線(Kn=1)

3.5 不同Kn下結(jié)果對比

綜合上面的計算結(jié)果，可以看出，隨著Kn的增加，繞圓柱的流場主要有以下變化:流場中的間斷，包括激波和滑移線都隨著Kn的增加而變厚和減弱，到了Kn=1 時，已經(jīng)沒有原來意義的激波;隨著Kn的增大，圓柱后方的分離區(qū)逐漸減小，直到消失.

4 結(jié)語

本文利用CFD－FASTRAN和ESI－RARE 軟件對四個典型Kn(Kn分別為0.001、0.01、0.1和1)，來流馬赫數(shù)為5 條件下的圓柱繞流進(jìn)行了計算，經(jīng)過分析比較發(fā)現(xiàn)，隨著Kn的增大，流場中的激波間斷結(jié)構(gòu)變?nèi)踝兇?，甚至消失，圓柱后方的分離區(qū)也逐漸變小直到消失.

［1］沈 青.DSMC 方法與稀薄氣流計算的發(fā)展［J］.力學(xué)進(jìn)展，1996，26(1):1－13.

［2］BIRD，GRAEME A.Molecular Gas Dynamics and the Direct Simulation of Gas Flows［M］.Clarendon Press:Oxford，1994.

［3］陳偉芳，吳其芬.Boltzmann 方程求解方法綜述［J］.國防科技大學(xué)學(xué)報，1999，21(1):4－7.

［4］樊 菁，劉宏立，沈 青，等.直接統(tǒng)計模擬位置元方法中的分子表面反射確定論判據(jù)［J］.空氣動力學(xué)學(xué)報，2000，18(2):180－187.

［5］李志輝，張涵信.稀薄流到連續(xù)流的氣體運(yùn)動論統(tǒng)一算法研究［J］.空氣動力學(xué)學(xué)報，2003，21(3):256－266.

［6］王 艷，胡煥林，姚達(dá)毛.TPMC和DSMC 方法在真空技術(shù)領(lǐng)域的應(yīng)用［J］.真空，2004，41(4):102－105.

［7］JIN D Q，WANG C，WEI Y J，et al.Axial cavitating flow study of underwater circular cylinder［J］.哈爾濱商業(yè)大學(xué)學(xué)報:自然科學(xué)版，2010，26(5):586－591.