一种计算非平衡等离子体中粒子能级布居的简化方法*

何新 江涛 高城† 张振福 杨俊波

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

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

获得粒子能级布居是研究非平衡等离子体辐射性质的一个重要方面. 对于复杂三维等离子体, 采用细致碰撞辐射模型虽然精确, 但计算耗费大. 本文提出了一种束缚态特征温度法, 能够快速计算得到非平衡等离子体中的粒子能级布居. 对非平衡氖等离子体算例的研究表明, 本文方法是有效的, 在等离子体非平衡程度不太高时与碰撞辐射模型符合较好. 在计算效率上, 本文方法比碰撞辐射模型至少提高了3000倍, 可极大节约计算资源和成本, 在工程计算中有重要实际意义.

1 引 言

在天体物理、X射线激光物理和约束聚变等领域中[1], 常常需要掌握等离子体的辐射特性, 进而研究其中辐射的输运与分配[2]. 为了获得辐射参数,必须知道等离子体中粒子的能级布居.

对局域热动平衡(local thermodynamic equilibrium, LTE)等离子体, 通过求解Saha方程可方便地得到粒子能级布居[3]. 然而, 很多情况下的等离子体处在非局域热动平衡状态(non-LTE, NLTE)[4].计算NLTE等离子体中的粒子能级布居是辐射特性研究必须面对的问题.

碰撞辐射(collisional-radiative, CR)模型是常用的NLTE等离子体的计算模型, 该模型考虑等离子体中所有碰撞和辐射等微观原子过程, 建立能级布居速率方程组并进行求解[5]. 根据原子参数的精密程度, 研究者发展了基于不同层次的CR模型, 如平均原子模型[6-9]、超组态模型[10-14]、细致组态模型[15-21]和细致能级模型[22,23]等. 一般说来,细致能级模型因为考虑了不同电荷态离子的能级结构, 计算精度最高, 但是计算耗费大; 而超组态模型把大量能量相近的精细能级近似处理成一个“能级”, 虽节省了计算时间, 但降低了计算精度[24-26]. 因此研究者们提出了一些兼顾计算精度和计算时间的CR模型和方法[27-29]. 例如Hansen等[27]使用混合了细致能级模型和超组态模型的方法获得原子参数和能级布居; Bauche等[29]在超组态模型中引入“有效温度”的概念, 能够方便获得相应能级上的粒子占据数. 虽然CR模型的种类很多[26], 但很多时候为了准确模拟NLTE等离子体辐射性质, 需要选取足够多数量的量子态, 使得速率方程组规模很大. 大规模的速率方程组求解复杂, 计算量大, 通常要借助高性能计算平台. 在三维尺度、参数梯度大的等离子体工程计算中, CR模型虽然精度高, 但难于实际应用. 因此, 研究能够保证一定精度同时低成本、高效率的计算方法,一直是研究者不断追求的目标, 也具有重要且实际的意义.

本文提出了一种束缚态特征温度简化方法, 用于NLTE等离子体中粒子能级布居的快速、简便计算. 以5种条件下的NLTE氖等离子体为例, 与CR模型的计算结果进行对比, 分析讨论该方法的准确度和效率.

2 束缚态特征温度法

本文不考虑处于完全热非平衡态的等离子体.在一般情况下的NLTE等离子体中, 可认为自由电子服从某一自由电子温度Te的麦克斯韦速度分布[2]. 考察NLTE等离子体中原子序号为A的z价粒子Az+( 0 ≤z≤A,z=0 表示原子), 设其电离能为Iz, 如图1所示.

图1 A z+ 的束缚态及连续态Fig. 1. Bound and continuum states of A z+.

一方面,Iz能级可看作Az+的束缚态. 引入Az+的束缚态特征温度, 即假设Az+各束缚态上的占据数服从温度为的Boltzmann分布, 则Iz能级上的“粒子数”满足:

其中Nz为Az+粒子数,为Iz能级的简并度,k为Boltzmann常数,为Az+在特征温度下的配分函数(和分别为Az+第j束缚态的能量和简并度).

另一方面,Iz能级上的粒子可看作由处于基态的A(z+1)+与速度为零的自由电子构成. 由于自由电子服从Te下的麦克斯韦速度分布, 则还满足[30]:

其中分别为A(z+1)+基态的占据数和简并度,Ne为自由电子数,me为自由电子质量,h为Planck常数. 类似地, 引入A(z+1)+的束缚态特征温度, (2)式可转化为

其中Nz+1为A(z+1)+粒子数,为A(z+1)+

联立(1)式和(3)式, 有

由此可计算出. 进而, 根据Boltzmann公式可得到Az+粒子的各能级占据数.

在应用上述方法时, 除需要粒子的能级参数外, 还需要各价粒子数密度、自由电子数密度和自由电子温度作为输入参数(来自于实验或计算). 具体计算过程中, 可直接令(因为z=A对应于粒子完全电离), 然后按照“倒序”依次求解出详细计算流程如图2所示.实际上, 总存在某个电离度范围p≤z≤q, 不在该电离度范围的粒子丰度相对很小. 此时, 可近似认为然后依次求解出进而计算主要粒子的能级布居.

图2 计算流程图Fig. 2. Calculation flowchart.

值得一提的是, 对于LTE等离子体, 必然有那么(4)式自然过渡为LTE条件下的Saha方程[30].

3 算例与讨论

表1列出了5种条件的NLTE氖(Ne)等离子体, 其中e—代表自由电子. 为了便于将本文计算的能级布居与CR模型结果进行直接对比, 表中Ne粒子百分数、自由电子(e—)数密度均采用了与CR模型相同的数据[22]. 另外, 表中未列出丰度小于10—6的粒子. 算例1—算例3主要研究Ne核总数密度相同、Te不同的情况; 算例4和算例5分别与算例2和算例3对应, 主要研究Te相同、Ne核总数密度不同的情况.

表1 算例参数Table 1. Cell parameters.

根据(4)式, 计算能级布居还用到各价Ne粒子的电离能(见表2). 等离子体环境中, 由于屏蔽效应, 原子(或离子)的电离能相比于孤立原子会下降. 表2中的数据由计算而得(等离子体条件:kTe为40 eV, Ne核数密度为1020cm—63). 在本文所研究的等离子体条件下, 电离能下降值变化不大, 因此忽略不同条件(自由电子温度、核总数密度)对Ne粒子电离能的影响.

表2 Ne原子及离子电离能Table 2. Ionization energy of neon atom and ions.

在下文各图中, 能级布居均是非简并的, 标有“Boltzmann”和“Saha”的数据分别由Te下的Boltzmann公式和Saha方程计算而得:

为反映等离子体的非平衡程度, 引入(5)式和(6)式的比值α(对同种粒子的任意能级均相等).显然, 对平衡等离子体,α=1 ; 反之,α偏离1. 为评估本文计算结果与CR模型的符合程度, 对各能级计算出二者的比值δ(δ=1 表示完全符合).

图3为算例1条件下本文计算的能级布居与CR模型结果的对比. 图中, 1.00＜α≤2.13 , 可知等离子体是弱非平衡的. 如图3(a), (b), (d), (e),对于Ne, Ne+, Ne3+和Ne4+,δ与1的平均偏差不大于10%, 说明本文计算的能级布居与CR模型符合很好. 如图3(c)所示, 对于Ne2+,δ与1的平均偏差稍大(—31%), 这可能主要是由于Ne2+的电离能数据存在误差.

图3 算例1的非简并能级布居Fig. 3. Non-degenerate electronic level populations for Case 1.

图4为算例2条件下能级布居结果对比. 图中, 1.71≤α≤30.08 , 较算例1有所增大, 可知此条件下等离子体的非平衡程度提高. 对于各价Ne离子(除图4(d)中的Ne5+),δ与1的平均偏差都在17%以内, 说明此条件下计算出的能级布居与CR模型符合较好.

图4 算例2的非简并能级布居Fig. 4. Non-degenerate electronic level populations for Case 2.

图5 为算例3条件下能级布居结果对比. 此条件下, 7.63≤α≤40484 , 说明非平衡程度进一步增大.如图5(a)所示, 对于Ne5+,δ与1的平均偏差为38%. 如图5(b)—图5(d), 对于Ne6+, Ne7+和Ne8+,δ与1的平均偏差不大于20%. 结果表明, 此条件下本文计算得到的能级布居与CR模型偏离不大.

对比图3—图5可知, 在Ne核总数密度为1018cm—3条件下, 对于弱非平衡态, 本文计算的能级布居与CR模型符合很好; 随着非平衡程度增大, 本文计算结果与CR模型之间逐渐产生偏离.

图6为算例4条件下能级布居计算结果. 对于各价Ne离子,α均接近1, 说明等离子体处于近平衡态. 此条件较算例2的非平衡程度降低, 主要归因于Ne核总数密度增大. 图6中,δ与1的平均偏差不超过10%, 说明计算结果与CR模型符合很好. 对比图6与图4可知, 在kTe= 15 eV条件下,若非平衡程度提高, 则本文得到的能级布居与CR模型的偏离增大, 这与图3—图5反映的现象类似.

图5 算例3的非简并能级布居Fig. 5. Non-degenerate electronic level populations for Case 3.

图6 算例4的非简并能级布居Fig. 6. Non-degenerate electronic level populations for Case 4.

图7 为算例5条件下能级布居计算结果. 相比于算例3, 由于Ne核总数密度增加, 等离子体非平衡程度降低( 1.03≤α≤218 ). 对于图7(e)中的Ne7+,δ与1的平均偏差稍大(—13.5%); 对于图7中其他各价Ne离子,δ与1的平均偏差都不大于7.8%. 因此, 此条件下得到的能级布居与CR模型符合较好. 对比图7与图5发现, 现象与kTe=15 eV时类似.

图7 算例5的非简并能级布居Fig. 7. Non-degenerate electronic level populations for Case 5.

综上, 若以CR模型作为参照, 本文提出的束缚态特征温度法的计算准确度与等离子体非平衡程度有关. 当等离子体非平衡程度较弱(1.00＜α≤2.13)时, 本文方法与CR模型非常符合(平均偏差基本在10%以内); 随着非平衡程度的增大( 2.13＜α≤218 ), 与CR模型逐渐产生偏差(平均偏差基本在30%以内); 若非平衡程度不太高(α≤40484), 与CR模型符合较好(平均偏差不超过38%).

本文方法也适用于其他等离子体中的原子及离子. 如果原子及离子来源于不同种核, 可分别对每种核应用上述方法. 由于不需求解能级布居速率方程组, 可极大节约计算资源. 表3给出了5个算例总耗时对比, 忽略程序语言效率及计算平台性能的差别, 本文方法比CR模型至少提高了3000倍.在误差允许范围内, 这将在参数梯度大、复杂三维非平衡等离子体的工程计算中产生巨大效益.

表3 计算耗费对比Table 3. A comparison of calculation cost.

4 结 论

提出了一种用于计算非平衡等离子体中原子及离子能级布居的简化方法. 以几种条件(5 eV ≤kTe≤ 40 eV, 1018cm—3≤ 核数密度 ≤ 1020cm—3)的非平衡Ne等离子体为例, 计算研究了粒子的能级布居, 并与CR模型进行了对比. 结果表明:1)该方法是有效的, 当等离子体非平衡程度不太高时与CR模型符合较好; 2)计算速度比CR模型至少提高了3000倍, 可大大节约计算成本, 对于复杂三维等离子体计算非常有用.