束缚态特征温度方法及应用*

何新 江涛 张振福 杨俊波

1) (国防科技大学文理学院,长沙 410073)

2) (中国空气动力研究与发展中心计算空气动力研究所,绵阳 621000)

随着高超声速飞行器速度增大,激波层空气等离子体中的原子发射谱线成为辐射加热主要来源,因此研究原子激发非常重要.考虑到处于热非平衡态的空气等离子体,平衡态统计理论不适用.精细物理模型(如碰撞辐射模型)虽然可以处理热非平衡问题且准确度高,但计算量太大,难于工程应用.本文采用束缚态特征温度法,结合FIRE II 激波管实验中的非平衡空气等离子体,对原子激发进行了分析.计算得到的原子能级布居与碰撞辐射模型符合,说明简化计算是合理的,计算效率提高了2000 倍以上,且能够保证一定的精度.

1 引言

在载人航天、太空探测等应用中,高超声速飞行器返回时被大强度激波包裹,激波层中的空气将变为等离子体[1].当飞行器速度很高时,其所承受热负荷的很大一部分源于激波层的辐射加热[2].为了评估辐射加热、指导热防护设计,必须研究空气等离子体的辐射特性[3,4].

分析原子激发是这其中的一个十分重要的方面,因为在高超声速激波层空气等离子体中,原子的辐射将占据主导地位[5,6].例如阿波罗飞行器返回时,约90%的激波层辐射来自于原子发射谱线[7].

在高超声速激波层中,部分空气等离子体处于热平衡态,部分处于热非平衡态[8].对于热平衡空气等离子体可以方便地利用Boltzmann 分布或Saha方程得到原子能级布居[9].然而,对于热非平衡空气等离子体,计算其中的原子能级布居是一项挑战.

特别是对于高超声速航天器大尺度三维激波层计算,其中常常包含超大量的热力学状态点(有限元),在这种情况下,已发展的、较常用的碰撞辐射(collisional-radiative,CR)模型[3,4,8−11],虽然能够处理热非平衡问题,且准确度高,但计算耗费超大,甚至无法实现[12].

为了既能处理热非平衡问题,又在误差可接受的前提下降低计算成本,研究者开发了一些简化计算方法.例如,采用准稳态(quasi-steady-state,QSS)近似、多温度Boltzmann 分布等来计算得到空气等离子体的原子激发数据[13−17],以便与大型飞行器流场计算耦合应用.然而,目前这些简化方法普遍基于精细CR 模型的思想,需要用到许多微观粒子相互作用速率系数,这些系数难以保证准确,从而导致这些简化方法之间、它们与CR 模型之间存在较明显的偏差[18,19],直接影响到飞行器辐射加热评估和热防护策略选择.可以说,针对高超声速大尺度三维等离子体计算需求,能够保证一定精度且快速计算分析原子能级布居,是研究者非常关注的课题和不断追求的目标.

在前期工作中,所提出的束缚态特征温度法无需用到微观粒子相互作用系数,在激光等离子体相互作用类似情形的计算中得到验证.本文在该方法基础上,对高超声速空气等离子中的原子激发进行研究,尝试发展适用于高超声速激波层等离子体的解决途径.选取高超声速实验中典型的热非平衡和热平衡空气等离子体作为研究对象,计算其中氮和氧原子能级布居,并与CR 模型、其他简化模型的结果进行对比,分析计算准确度和计算效率.

2 计算方法

本文关注广泛存在于载人航天器返回、太空探测器再入等航天应用领域的空气等离子体状态,其自由电子温度一般不高(30000 K 以下),且为弱电离.通常在这样的空气等离子体中,只需考虑原子(N,O)及其一价离子(N+,O+)的存在[15].

根据束缚态特征温度法[20],若设某种原子的电离能为I,则表征该种原子能级布居的特征参数Tb可由下式计算:

其中Tb为该种原子的束缚态特征温度,Te为自由电子温度,n,n+和ne分别为原子、一价离子和自由电子的数密度,me为自由电子质量,k和h分别为Boltzmann 常数和Planck 常数,Q(Tb)和Q+(Te)分别为原子和一价离子的配分函数:

其中Ei和gi分别是原子第i能级的能量和简并度,分别是一价离子第j能级的能量和简并度.

若已知n,n+,ne和Te,可根据(1)式计算出Tb,从而原子第i能级的非简并布居为：

实际上,(1)式可称为修正的Saha 方程.这样写的好处是对热平衡和热非平衡空气等离子体都适用.对于热平衡空气等离子体,必然存在Tb=Te,则(1)式就是众所周知的Saha 方程[1];对于热非平衡空气等离子体,所求解出的Tb将与Te不同,二者之间的差别反映了能级布居偏离热平衡分布的程度.

3 典型空气等离子体算例及结果

选取针对Fire II 工程的地面激波管实验空气等离子体为研究对象[21],其中的辐射主要来源于原子谱线,对原子能级布居(尤其是较高能级布居)计算非常重要.表1 给出了所选取的空气等离子体状态参数(包括自由电子温度、粒子数密度),它们分别对应Fire II 飞行过程中的1634,1636 和1643 s 时间点.

表1 空气等离子体参数Table 1.Parameters of air plasmas.

由(1)—(3)式可知,计算能级布居还需要原子及其离子的能级参数.本文中,N 和O 原子能级数据来自文献[22],N 和O 原子电离能采用美国国家标准与技术研究院(NIST)数据[23],N+和 O+离子能级数据来自文献[1].

图1 为状态点1634-25 的氮原子非简并能级布居.标有“Johnston”的数据由文献[7]中方法计算而得;标有“Boltzmann”的数据由Te下的Boltzmann 分布,即(3)式中的Tb换成Te计算而得(下同);标有“Saha”的数据由另一形式的Saha 方程计算而得(下同):

众所周知,若等离子体处于热平衡态,则Te下的Boltzmann 分布与Saha 方程的结果必然是相等的.图1 中,标有“Boltzmann”和“Saha”的数据只存在轻微偏离,说明该状态点空气等离子体处于近平衡态.根据图1,对于N 的低能级和高能级布居,计算结果与CR 模型一致;对于N 的中间一些能级布居,计算结果介于CR 模型和Johnston 方法.

图1 状态点1634-25 的氮原子能级布居Fig.1.Energy level populations for N of Case 1634-25.

图2 给出了状态点1634-10 的氮原子和氧原子非简并能级布居.“Boltzmann”和“Saha”数据之间存在明显偏离,说明此状态点空气等离子体为非平衡.这是由于距离激波面较近,粒子之间能量松弛不充分造成的.

如图2(a)所示,除N 原子的第2,3 能级外,本文计算得到的其他能级布居与CR 模型一致;而此条件下采用Spradian 模块得到的结果明显大于CR 模型.如图2(b)所示,除O 原子的第2—4能级外,本文计算得到的其他能级布居与CR 模型符合;而Spradian 模块的结果明显大于CR模型.

CR 模型是更为精细的物理模型,准确度更高.图1 和图2 的结果表明,若以CR 模型为参照讨论计算精度,则本文计算精度与Johnston 方法接近,优于Spradian 模块.

图3 为状态点1634-7 的氮原子非简并能级布居.显然,此状态点为非平衡空气等离子体.如图3(a),计算得到的第2,3 能级占据数略小于CR 模型,但其他能级布居与CR 模型一致,这与图2 中的现象类似.图3(b)中还与其他简化模型的结果进行了对比.可以看出,计算结果与Johnston 方法符合,在一些中间能级与QSS Abba 方法符合,但整体上明显低于QSS Park 方法.如果仍然参照CR 模型讨论计算精度,图3 的结果表明,在此状态条件下,本文计算精度与Johnston 和QSS Abba 方法接近,优于QSS Park 方法.

图2 状态点1634-10 的氮原子和氧原子能级布居 (a) N;(b) OFig.2.Energy level populations for N and O of Case 1634-10:(a) N;(b) O.

图3 状态点1634-7 的氮原子能级布居Fig.3.Energy level populations for N of Case 1634-7.

图4 为状态点1636-5 的氮原子和氧原子非简并能级布居.此状态点空气等离子体虽然也是非平衡的,但由于粒子数密度较大,粒子之间能量松弛较图3 所示情形充分,因此“Boltzmann”和“Saha”数据之间的偏离减小.对于N 的中间一些能级,计算得到的占据数略大于CR 模型.但总体来说,在此状态点,计算结果与CR 模型符合.

图4 状态点1636-5 的氮原子和氧原子能级布居 (a) N;(b) OFig.4.Energy level populations for N and O of Case 1636-5:(a) N;(b) O.

图2—图4 中,本文所计算N 原子第2—3 能级、O 原子第2—4 能级的占据数比CR 模型偏小,对辐射加热评估的影响较小.这是因为N 和O 原子的发射谱线并不包含上述能级的自发辐射跃迁[21].另外,对于N 和O 原子几个中间能级的占据数,本文虽然与CR 模型存在一定偏差,但作为一种简化计算结果,较其他简化方法在计算精度上持平或有明显提高.

图5 为状态点1643-5 的氧原子非简并能级布居.此状态点对应Fire II 的飞行高度低,激波层中粒子能量松弛充分,因此空气等离子体处于近平衡态.由图5 可知,计算结果与CR 模型一致.

图5 状态点1643-5 的氧原子能级布居Fig.5.Energy level populations for O of Case 1643-5.

上述5 个不同状态点的结果表明,本文计算得到的空气等离子体原子能级布居与CR 模型基本一致,是合理有效的.在计算精度上,本文与QSS Abba和Johnston 简化方法接近,优于QSS Park 和Spradian 简化方法.

在工程应用中,计算效率是必须考虑的问题.借助普通笔记本电脑(CPU:2×2.60 GHz,Matlab程序),计算得到表1 所有状态点的能级布居数据约需7 s;而CR 模型借助IBM 服务器(CPU:6×2.53 GHz,Fortran 程序)约需4.5 h.即使忽略计算平台性能、程序语言效率、算法流程设计方面的差别,计算速度也比CR 模型提高了2000 倍以上.本文计算效率的提升主要源于不需要求解能级布居方程组.CR 模型尽可能多地考虑了影响粒子能级布居的大量微观过程,需要求解大规模的能级布居速率方程组(方程数量取决于所考虑能级数);而文中其他简化方法虽然减少了所考虑的微观过程,甚至进一步采用QSS 近似,但仍然需要求解能级布居方程组.因此,本文计算可与高超声速飞行器流场计算耦合,极大降低计算成本.

4 结论

采用束缚态特征温度法研究了高超声速激波层空气等离子体中的原子激发.以针对Fire II 工程的激波管实验空气等离子体为算例,对N 和O 原子能级布居进行了计算,并与CR 模型、其他简化模型的结果进行了对比.所研究的空气等离子体热力学状态包括近平衡态、非平衡态.结果表明,本文计算得到的原子能级布居与CR 模型基本一致,计算精度与其他简化方法接近甚至有一定提高.在计算效率上,本文比CR 模型提高了2000 倍以上,在工程应用(如高速飞行器辐射加热评估)中可大大节约计算成本.