何 亮,姜東君,應(yīng)純同

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

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

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

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

1 DSMC方法概述

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

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

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

2 計算程序

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

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

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

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

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

3 計算模型和邊界條件

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

圖1 剛體平衡解計算模型

4 結(jié)果與討論

4.1 等溫剛體的模擬分析

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

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

圖2 分子數(shù)密度的流場分布圖分子密度數(shù)水平: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 分子數(shù)密度沿半徑的變化比較圓點為計算數(shù)據(jù);實線為理論數(shù)據(jù)

圖4 旋轉(zhuǎn)速度沿半徑的變化比較圓點為計算數(shù)據(jù);實線為理論數(shù)據(jù)

4.2 側(cè)壁熱驅(qū)動的模擬分析

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

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

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

5 結(jié) 論

本工作對G2程序進行了改進,并利用改進后的程序?qū)崿F(xiàn)了強旋流場的DSMC模擬,通過對等溫剛體運動和側(cè)壁溫度為線性分布等兩種情況的模擬,可以看出改進后的G2程序能模擬離心機流場分布,并且能得到與以往文獻較吻合的結(jié)果。因此,可以使用DSMC方法進一步研究較為復(fù)雜的高速旋轉(zhuǎn)氣體離心機內(nèi)部流場。

[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] 張存鎮(zhèn).離心分離理論[M].北京:原子能出版社,1987.