一種計算非平衡等離子體中粒子能級布居的簡化方法*

何新 江濤 高城? 張振福 楊俊波

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倍, 可大大節約計算成本, 對于復雜三維等離子體計算非常有用.