何 亮,姜东君,应纯同

(清华大学工程物理系,北京 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.