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

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

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

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

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

1 DSMC方法概述

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

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

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

2 計(jì)算程序

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

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

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

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

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

3 計(jì)算模型和邊界條件

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

圖1 剛體平衡解計(jì)算模型

4 結(jié)果與討論

4.1 等溫剛體的模擬分析

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

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

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

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

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

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

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

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

5 結(jié) 論

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

[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.