DSMC方法模擬強旋流場的初步研究

何 亮,姜東君,應純同

(清華大學工程物理系,北京 100084)

離心法分離同位素的主要原理是把同位素氣體混合物置于高速旋轉的離心機轉子中,在強大的離心力場中,使相對分子質量不同的組分分離開來。離心機流場的壓強在徑向上呈指數分布,壓強梯度隨半徑和旋轉速度的增大而增大。在離心機的中心區域,氣體非常稀薄,不能使用Navier-Stokes來描述其流場[1]。因此在傳統計算中,離心機內流場計算都限制在靠近壁面附近,滿足Navier-Stokes方程的區域,通過引入粘性邊界,在粘性邊界上設定自然邊界條件,并人為給定供料條件。由于該假設不能反映真實的流動狀況,所以需要對離心機中心區域流場進行模擬計算。

在稀薄氣體的模擬方法中,Bird[2]提出的直接蒙特卡羅方法(Direct Simulation Monte-Carlo,DSMC)是應用最廣泛的數值模擬方法。Roblin[3]等使用直接蒙特卡羅方法模擬了離心機供料流的流場,得到了供料射流存在的離心機流場分布圖,但是并沒有給出流動的更多細節,也沒有進一步討論稀薄區氣體運動對離心機環流的影響,因此,有必要對離心機中心區域進行深入研究。

本工作擬采用直接蒙特卡羅方法(DSMC)方法對強旋條件下的離心機流場進行數值模擬,以探討DSMC方法對復雜結構的離心機流場的求解能力。

1 DSMC方法概述

DSMC是一種通過跟蹤有限個模擬分子的運動來研究真實氣體流動的數值模擬方法。模擬分子的速度、位置信息和內能等物理信息都存儲在計算機中,并且在模擬計算區域中,根據分子運動、分子間的碰撞以及分子與邊界的作用而隨時間不斷更改這些微觀物理信息。

DSMC方法的重要思想是將分子的運動和碰撞解耦。在分子運動階段,模擬分子根據其本身的速度和計算時間 Δt更新位置信息,這些信息包括和邊界及物面的相互作用,根據不同的模型計算模擬分子作用后的新位置。在分子碰撞階段,分子碰撞取樣采用非時間計數器(No-Time-Counter)方法,在每一個網格單元中,恰當選擇碰撞對,并實現一定數目的碰撞,使Δt時間內碰撞與運動匹配,模擬流動與真實流動一致。

分子的運動和碰撞解耦的思想,即在 Δt內所有模擬分子按照速度運動一段時間,然后計算Δt時間內有代表性的分子間發生碰撞的過程,使得DSMC程序有隱含限制,時間步長Δt必須遠小于局部平均碰撞時間、網格長度必須小于局部平均分子自由程、每個網格單元必須包括一定數量的模擬分子以保證計算結果在統計上處于一定的概率范圍。

2 計算程序

本工作在Bird[2]提供的G2程序基礎上進行了發展和改進,使其適合于更復雜的流場特點。在程序修改過程中,增加了可設置邊界的總數,將計算程序從單精度運算改成雙精度運算,同時注意Fortran程序相關優化建議,注重提高程序的效率。

為了減少分子平均碰撞距離,在G2程序中采用固定子單元格(Fixed Sub-Cell)技術,將計算網格劃分為單元格和子單元格,將碰撞分子的選擇限制在相同或相鄰子單元格中(都屬于同一單元格)。這樣有利于減少碰撞分子對間的距離,但其代價是增加了計算時間和內存。根據文獻[3-5]的建議,在分子碰撞過程中,對即將要計算的單元格進行子單元格劃分,并將此單元格中的分子按子單元格順序進行排序,這就是動態子單元格(Transient-Adaptive Sub-Cell)技術。這樣可大幅減少計算時間和內存。文獻[6]建議采用虛擬單元格(Virtual Sub-Cell)技術,在網格中隨機選擇一個分子,計算與其他分子的距離,選擇距離最近的分子作為碰撞分子。本工作擬基于動態子單元格技術對程序進行修改,以減少內存需求并加快運行速度。

在軸對稱流場中,假如網格是均勻劃分,那么靠近軸的單元格體積會遠小于靠近壁面的單元格的體積,因此分子數會很少。為了解決這個問題,在G2程序中引入了徑向加權系數W F,其表達式如下:

(1)式中,r為分子的徑向半徑,r max為流場區域的最大半徑,R wf為加權常數。r處每個模擬分子所代表的真實分子數NUM用(2)式表示:

(2)式中,FNUM為每個模擬分子所代表的平均真實分子數,由于加權系數不能為0,因此上述加權系數計算公式就存在截斷半徑,這會影響靠近軸線處的宏觀量統計,進而導致G2程序無法準確求解強旋流場。文獻[7]建議使用如下的徑向加權系數公式:

3 計算模型和邊界條件

文獻[8]認為,在旋轉圓柱體內,氣體運動的熱力學平衡解是等溫剛體運動,這種無能無運動可以滿足Boltzmann方程。其計算模型示于圖1。由圖1可以看出,邊界設置為:上邊界為側壁,左右兩端邊界為端蓋。在等溫剛體算例中,壁面條件是8 000 rad/s,溫度為300 K;在側壁溫度驅動的兩個算例中,左右端面溫差為20 K,側壁溫度為線性分布,壁面處角速度分別為7 000和8 000 rad/s。壁面反射模型均為完全漫反射模型,工作氣體為UF6。本程序中氣體分子模型為變徑硬球模型,碰撞采樣采用Bird提出的NTC方法。

圖1 剛體平衡解計算模型

4 結果與討論

4.1 等溫剛體的模擬分析

在旋轉角速度為8 000 rad/s下,離心機流場的分子數密度分布示于圖2。從圖2可以看出,分子數密度沿徑向存在變化,沿軸向基本沒有變化,并且在靠近壁面處,變化較劇烈。

分子數密度和旋轉速度在某一軸向位置沿徑向的分布變化分別示于圖3和圖4。由圖3和圖4可以看出,速度沿徑向線性變化,分子數密度沿徑向呈指數分布,且和理論結果基本一致。因此,可以認為,流場分布滿足等溫剛體分布。

圖2 分子數密度的流場分布圖分子密度數水平:1——5×1019;2——1×1020;3——2×1020;4——3×1020;5——4×1020;6——5×1020;7——6×1020;8——7×1020;9——8×1020;10——9×1020

圖3 分子數密度沿半徑的變化比較圓點為計算數據;實線為理論數據

圖4 旋轉速度沿半徑的變化比較圓點為計算數據;實線為理論數據

4.2 側壁熱驅動的模擬分析

文獻[9]認為,當側壁溫度為線性分布時,離心機流場會存在環流。本工作將側壁兩端溫差設定為20 K,通過DSMC模擬計算,可以得到離心機的軸向環流流線圖。圖5和圖6為不同旋轉速度下離心機流場的軸向速度分布和流線圖。從圖5和圖6可以看出,軸向速度分布沿軸向基本沒有變化;離心機旋轉速度越高,環流越靠近壁面。

圖5 7 000 rad/s旋轉速度下離心機流場的軸向速度分布和流線圖直線為軸向速度分布;曲線為流線

圖6 8 000 rad/s旋轉速度下離心機流場的軸向速度分布和流線圖直線為軸向速度分布;曲線為流線

5 結 論

本工作對G2程序進行了改進,并利用改進后的程序實現了強旋流場的DSMC模擬,通過對等溫剛體運動和側壁溫度為線性分布等兩種情況的模擬,可以看出改進后的G2程序能模擬離心機流場分布,并且能得到與以往文獻較吻合的結果。因此,可以使用DSMC方法進一步研究較為復雜的高速旋轉氣體離心機內部流場。

[1] Kai T.Basic characteristics of centrifuges(Ⅲ):analysis of fluid flow in centrifuges[J].J Nucl Sci Technol,1977,14(4):267-281.

[2] Bird GA.Molecular gas dynamics and the Direct Simulation of gas flows[M].New York:Oxford University Press,1994.

[3] Roblin P,Doneddu F.Direct Monte-Carlo simulations in a gas centrifuge[C]//The 22nd International Symposium on Rarefied Gas Dynamics.Sydney,Australia:AMER INST Physics,2000:169-173.

[4] Gallis MA,Torczynski JR,Rader DJ.Accuracy and convergence of a new DSMC algorithm[R].AIAA 2008-3913.Washington:AIAA,2008.

[5] Bird GA.The DS2V/3V program suite for DSMC calculations,in rarefied gas dynamics[C]//24th International Symposium,AIP Conference Proceedings 762.New York,2005:541-546.

[6] LeBeau GJ.Virtual Sub-Cells for the Direct Simulation Monte Carlo Method[R].AIAA 2003-1031.Washington:AIAA,2003.

[7] Liechty Derek S.Modifications to axially symmetric simulations using new DSMC(2007)algorithms[J].Rarefied Gas Dynamics,2009,1 084:251-256.

[8] Soga T,Ooue K.On the numerical simulation of rotating rarefied flow in the cylinder with smooth surface[C]//23rd International Symposium American Institute of Physics Conference Proceedings,Whistler:AIP,2003:210-217.

[9] 張存鎮.離心分離理論[M].北京:原子能出版社,1987.