盧肖勇 袁程 高陽

1) (核工業(yè)理化工程研究院, 天津 300180)

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

離子引出是激光法分離同位素技術中的關鍵環(huán)節(jié)之一, 對產品豐度和產率有重要的影響, 而離子引出中離子和原子間的共振電荷轉移過程則會對產品豐度造成污染, 因此在研究離子引出過程時應考慮共振電荷轉移的影響. 本文利用粒子模擬(PIC)法以及基于PIC法和雜化PIC法的混合算法研究了考慮共振電荷轉移的電場法離子引出過程, 通過對一維平行板法的引出方式進行數(shù)值計算, 獲得了離子引出過程中共振電荷轉移的基本性質和關鍵影響因素—共振電荷轉移截面、引出時間、背景原子密度和蒸氣寬度, 并據(jù)此得到了離子發(fā)生共振電荷轉移比例的經(jīng)驗公式; 二維工況下的平行板法、交替偏壓法、Π型電場法和M型電場法四種引出方式的計算結果則表明, 在其他條件相同的情況下, M型電場法具有最小的離子引出時間和共振電荷轉移損失比例. 本文的研究結論對于指導激光法分離同位素技術中離子引出裝置設計和實驗工藝設計具有比較重要的參考意義.

1 引 言

離子引出過程是激光法分離同位素技術(laser isotope separation, LIS)的核心環(huán)節(jié)之一, 由光電離產生的離子的引出率直接關系到最終產品的產率, 同時離子引出過程中的各種碰撞過程則會對產品的產率和豐度組成造成影響[1]. 其中, 關于離子引出率的研究結論已經(jīng)比較明確, 之前的學者通過理論和實驗的手段獲得了諸多影響離子引出率的關鍵因素[2-9]: 引出電壓、等離子體密度、等離子體尺寸、電子溫度、磁場[10-14]等等. 但是關于其中碰撞過程的研究則相對較少, 研究結論也有較大的差異[15,16], 因此考慮碰撞的離子引出過程還有待進一步的研究.

離子引出過程中主要存在以下碰撞過程: 光電離生成的離子和未電離原子間的碰撞(共振電荷轉移過程)、電子與未電離原子間的碰撞(二次電離過程)、光電離生成的電子與離子的碰撞(復合過程).其中復合過程會使精料量下降, 共振電荷轉移過程和二次電離過程則會影響產品的同位素豐度組成,從而影響分離效果. 由于二次電離過程和復合過程對離子總量的影響遠小于共振電荷轉移過程的影響[15], 因此在本文的研究中僅考慮共振電荷轉移過程對離子引出的影響.

本文采用PIC (particle in cell)法、以及基于PIC法和雜化PIC法的混合算法研究了一維和二維情況下考慮共振電荷轉移的離子引出過程. 在一維模擬中重點關注了共振電荷轉移過程的特征和主要影響因素, 并根據(jù)模擬結果和理論分析給出了離子發(fā)生共振電荷轉移比例的經(jīng)驗公式; 在二維模擬中分析了平行板電場法、交替偏壓法、Π型電場法和M型電場法等四種常見構型下的離子引出過程, 獲得了各構型下離子發(fā)生共振電荷轉移的比例, 并基于離子引出效率和共振電荷轉移比例對上述構型進行了對比評價.

2 數(shù)學模型

PIC方法是一種直接從第一性原理出發(fā)的等離子體粒子模擬方法[17,18], 其模擬的流程圖如圖1所示. 在完成模擬空間的參數(shù)設置、網(wǎng)格劃分和粒子信息初始化之后, 計算基本流程如下:

圖1 PIC法計算流程圖Fig. 1. Schematic diagram of PIC method.

1) 將網(wǎng)格內宏粒子的電荷分配到網(wǎng)格點上獲取網(wǎng)格點的電荷和電流;

2) 根據(jù)網(wǎng)格點上的電荷和電流以及邊界條件求解電磁場方程獲取網(wǎng)格點上的電磁場;

3) 根據(jù)網(wǎng)格點上的電磁場插值得到宏粒子的受力;

4) 根據(jù)經(jīng)典運動方程或相對論計算得到宏粒子在下一時刻的速度和位置;

5) 處理共振電荷轉移.

雜化PIC方法則是在PIC方法的基礎上假設電子處于平衡態(tài), 僅計算離子的運動過程. 由于不需要跟蹤電子的運動, 因此雜化PIC法的計算量遠小于PIC法. 混合算法是作者在處理電場法離子引出問題時提出的一種綜合PIC法和雜化PIC法的方法, 相較于PIC法, 混合算法的計算效率更高, 同時比雜化PIC法的適用范圍更大, 該算法的詳細介紹可參見文獻[19].

離子引出過程中, 在一個模擬時間步長Δt內,離子與背景原子的共振電荷轉移幾率Pia為:

其中,σT,ia(via)為相對速度為via的離子和原子的共振電荷轉移截面,na(r)為位置r處原子的數(shù)密度. 為了提高數(shù)值計算效率, 帶電粒子和中性粒子的共振電荷轉移采用零碰撞(null collision)技術以避免遍歷每個宏粒子.

3 模擬結果及分析

3.1 一維離子引出過程

本節(jié)以一維平行板電場法離子引出為例, 其示意圖如圖2所示. 在未做額外說明的情況下, 均采用如下的計算條件: 極板間距L= 5 cm; 在未電離前, 金屬蒸氣在極板間均勻分布, 其原子數(shù)密度為4.0 × 1011cm—3; 金屬蒸氣中同時存在A, B兩種同位素(其中同位素A為目標同位素), 其豐度均為50%, 兩種同位素的相對原子質量均取為200(這種取法不會影響最后的結論). 等離子體寬度為Lp= 3 cm, 位置位于極板中心, 其初始密度分布為均勻分布. 電子溫度為0.5 eV, 背景氣體中的原子和電離產生的離子溫度均為0.1 eV. 本節(jié)共計算5個典型算例, 其中同位素B的電離率均為0, 同位素A的電離率、初始離子密度和相應的引出電壓等參數(shù)詳見表1. 共振電荷轉移截面取2.0× 10—14cm2, 雖然實際的共振電荷轉移截面與離子和原子相對速度有關[20], 但是這種取法不影響本文的結論. 計算時考慮了光電離導致的光作用區(qū)內同位素A原子數(shù)密度的降低.

圖2 一維平行板電場法離子引出示意圖Fig. 2. Schematic diagram of one dimensional electric ion extraction of parallel type.

表1 一維算例的計算條件Table 1. Simulation parameters in onedimensional cases.

計算時采用的網(wǎng)格數(shù)為1000, 網(wǎng)格寬度為5 × 10—5m, PIC法和雜化PIC法的時間步長分別為100 ps和1 ns; 該算法已在作者之前工作中通過檢驗空間步長、時間步長和粒子數(shù)的方式完成收斂性檢驗, 詳見文獻[19].

圖3所示為考慮共振電荷轉移條件下, 算例1計算區(qū)域內的剩余離子/電子比例隨引出時間的變化曲線, 為了方便比較, 圖3中同位素B顯示的結果為實際結果的100倍(下同). 圖4為離子沉積至兩側收集板時的能量分布, 其中圖4(a)的顯示范圍為0—2.5 keV, 圖4(b)為0—20.0 eV. 作為對比, 圖4中也給出了未考慮共振電荷轉移時的計算結果. 圖4(a)中的計算結果表明, 當考慮共振電荷轉移時, 同位素A, B的離子沉積能量分布可分為三部分: 第一部分離子的沉積能量集中在2.0 keV附近; 第二部分在0—2.0 keV之間近似均勻分布;第三部分集中在0 keV附近. 當不考慮共振電荷轉移時, 同位素A離子的沉積能量只有集中在0和2.0 keV附近的兩部分. 這種現(xiàn)象解釋如下.

圖3 剩余離子比例隨引出時間的變化曲線 (算例1)Fig. 3. Plots of remaining ion ratio versus extraction time(Case 1).

不考慮共振電荷轉移時, 同位素A離子在向負極板運動時被電場加速, 由于鞘層兩側的電勢差約為2.0 kV, 因此從負極板引出的離子能量集中在2.0 keV附近. 同時, 正極板附近也會形成微小的鞘層, 小部分被正極板引出的離子的能量為數(shù)電子伏, 在圖4(b)中則表現(xiàn)為沉積能量集中在0 keV附近. 當考慮共振電荷轉移時, 由于引出開始時等離子體中沒有同位素B離子, 因此引出過程中的同位素B離子均由共振電荷轉移產生. 由同位素B離子的沉積能量分布可以看出, 一部分共振電荷轉移發(fā)生在屏蔽層中, 因此被引出至收集板時離子的沉積能量集中在2.0和0 keV附近(原理同上);另一部分共振電荷轉移發(fā)生在鞘層中, 因此被引出至收集板時同位素B離子能量在0—2.0 keV之間. 同位素A離子的能量分布與同位素B類似, 與不考慮共振電荷轉移時的計算結果相比, 一部分同位素A離子的沉積能量分布在0—2.0 keV之間,這是由同位素A離子在鞘層中與同位素A原子發(fā)生共振電荷轉移作用所致.

圖4 兩側收集板上的離子沉積能量分布 (算例1) (a) 整體圖; (b) 局部放大圖Fig. 4. Energy distribution of ions deposit on the collection plates (Case 1): (a) General distribution; (b) local distribution.

由于發(fā)生共振電荷轉移的離子和原子只交換電荷, 不交換速度, 因此發(fā)生在鞘層中的共振電荷轉移相當于使鞘層中離子的平均速度減小了, 從而使離子引出速率略微變慢. 只是由于在鞘層中發(fā)生共振電荷轉移的離子比例較低, 因此引出速度的下降幅度十分微弱, 如圖3所示. 屏蔽層中的共振電荷轉移則無此影響. 同位素B剩余離子比例則受到其產生速率和引出速率的共同影響, 因此呈現(xiàn)先增加后減少的趨勢.

圖4中還有兩個現(xiàn)象需要解釋: 一是同位素A離子的沉積能量分布曲線中2.0 keV位置處的能量峰左側均有一個小峰, 這是由于離子引出過程中鞘層形成階段鞘層內電勢波動造成的[19], 同位素B離子的沉積能量曲線對應能量處未出現(xiàn)明顯峰型, 其原因是同位素B離子是由共振電荷轉移過程產生, 鞘層形成階段產生的同位素B離子量少, 因此其沉積能量曲線在2.0 keV附近無其他明顯的能量峰; 二是兩種同位素2.0 keV位置處的能量峰不是左右對稱的, 左側偏大, 且2.0 keV能量峰峰值處的離子能量小于2.0 keV, 這是由于鞘層剝離和離子在其中的加速運動同時進行, 離子在加速運動中感受到的電場會逐漸減小, 因此電場傳遞給離子的能量小于2.0 keV. 圖5所示為算例1改變不同初始等離子體密度時, 離子運動至兩側收集板時的能量分布. 可以看出, 隨著初始等離子體密度的增加, 鞘層剝離速度變慢, 因此能量峰的位置逐漸接近2.0 keV.

圖5 不同初始等離子體密度條件下, 兩側收集板上的離子沉積能量分布(算例1)Fig. 5. Energy distribution of ions deposit on the collection plates with several initial plasma densities (Case 1).

由于同位素A離子與A原子之間的共振電荷轉移不會影響最終產品的豐度, 因此以下部分只研究同位素A離子與B原子之間的共振電荷轉移過程. 圖6所示為引出電壓分別取1.0和2.0 kV(分別對應算例1和算例2)時, 剩余離子比例隨引出時間的變化曲線, 圖7給出了相應條件下極板上離子沉積能量分布曲線. 當引出電壓分別為1.0和2.0 kV時, 同位素A離子中發(fā)生共振電荷轉移的比例分別為1.278%和1.151%, 其中, 在屏蔽層中和鞘層中發(fā)生的比例分別為0.663%和0.615%(1.0 kV), 0.371%和0.780% (2.0 kV), 計算結果如表2所列. 結合圖6中的剩余離子比例隨引出時間的變化曲線和(1)式, 可以得到如下結論: 在其他條件不變的情況下, 隨著引出電壓的升高, 引出時間減小, 離子在屏蔽層中運動的路程減小, 從而使離子屏蔽層中發(fā)生共振電荷轉移的幾率下降(0.663% → 0.371%). 當引出電壓分別為1.0和2.0 kV時, 從負極板引出的離子的比例分別為70.24%和84.62%, 由此可得由負極板引出的離子在通過鞘層時發(fā)生共振電荷轉移的幾率從0.876%升至0.921%, 這是因為引出電壓的升高使更多的離子由負極板引出, 導致從負極板引出的離子在鞘層中的平均運動路程增長, 發(fā)生共振電荷轉移的幾率升高. 因此隨著引出電壓的增加, 通過鞘層的離子比例增加, 且離子在鞘層中的平均運動路程略有增長, 這兩個原因導致離子在鞘層中發(fā)生共振電荷轉移的幾率增大.

表2 各算例的數(shù)值計算結果和公式評估結果Table 2. Simulation and empirical formula results of several simulation cases.

圖6 不同引出電壓條件下, 剩余離子比例隨引出時間的變化曲線(算例1和算例2)Fig. 6. Plots of remaining ion ratio versus extraction time with several applied voltages (Case 1 and Case 2).

圖7 不同引出電壓條件下, 兩側收集板上的離子沉積能量分布(算例1和算例2)Fig. 7. Energy distribution of ions deposit on the collection plates with several applied voltages (Case 1 and Case 2).

圖8所示為背景金屬蒸氣原子數(shù)密度分別取2.0 × 1011, 4.0 × 1011和8.0 × 1011cm—3(分別對應算例1、算例3、算例4)時, 剩余離子比例隨引出時間的變化曲線. 圖9給出了相應條件下極板上離子沉積能量分布曲線. 在背景金屬蒸氣原子數(shù)密度改變時, 同位素豐度和初始離子密度均不變. 在上述三種背景原子蒸氣密度下, 同位素A離子中與B原子發(fā)生共振電荷轉移的比例分別為0.576%,1.151%和2.250%, 與背景原子蒸氣密度基本上成正比, 且背景原子蒸氣密度的變化不影響共振電荷轉移發(fā)生在鞘層中和屏蔽層中的比例.

圖8 不同背景原子蒸氣密度條件下, 剩余離子比例隨引出時間的變化曲線(算例1、算例3、算例4)Fig. 8. Plots of remaining ion ratio versus extraction time with several background atomic densities (Case 1, Case 3,Case 4).

圖9 不同背景原子蒸氣密度條件下, 兩側收集板上的離子沉積能量分布(算例1、算例3、算例4) (a) 整體圖; (b) 局部放大圖Fig. 9. Energy distribution of ions deposit on the collection plates with several background atomic densities (Case 1,Case 3, Case 4): (a) General distribution; (b) local distribution.

圖10 所示為算例5中剩余離子比例隨引出時間的變化曲線, 圖11給出了相應條件下極板上離子沉積能量分布曲線. 同位素A離子中與B原子發(fā)生共振電荷轉移的比例為1.203%, 其中發(fā)生在屏蔽層中和鞘層中的比例分別為0.549%和0.654%.考慮到離子從負極板引出的比例為75.78%, 即可得到由負極板引出的離子在通過鞘層時發(fā)生共振電荷轉移的幾率為0.864%. 另外, 由圖11中局部放大部分可以看出, 在5.0—8.0 keV之間, 同位素B的沉積速率略大于同位素A, 這是因為同位素A的電離率較大導致背景原子中同位素A的原子數(shù)密度低于同位素B, 在算例1—算例4中此種現(xiàn)象則不明顯.

圖10 兩側收集板上的離子沉積能量分布(算例5)Fig. 10. Plots of remaining ion ratio versus extraction time(Case 5).

圖11 兩側收集板上的離子沉積能量分布(算例5)Fig. 11. Energy distribution of ions deposit on the collection plates (Case 5).

通過以上計算結果可以發(fā)現(xiàn): 在電場法離子引出過程中, 屏蔽層和鞘層均會發(fā)生共振電荷轉移,根據(jù)(1)式, 除了共振電荷轉移截面, 離子發(fā)生共振電荷轉移的幾率只受到背景原子密度和離子相對原子運動路程的影響. 屏蔽層中離子和原子均以熱速度運動, 因此此區(qū)域內的離子發(fā)生共振電荷轉移的幾率主要取決于引出時間, 即可近似為離子停留在屏蔽層中的時間. 因此離子在屏蔽層中發(fā)生共振電荷轉移的比例Pshield可采用下式估算:

其中vT,ia為原子和離子間的平均相對熱運動速度,t為離子引出時間. 另外, 由負極板引出的離子在通過鞘層時發(fā)生共振電荷轉移的幾率Psheath可采用下式估算:

其中La為被負極板收集的離子在背景氣體中平均運動路程長度的2倍, 在估算時采用背景氣體寬度代替. 綜上, 離子引出過程中離子發(fā)生共振電荷轉移的比例Ptotal可用下式估算:

其中αcathode為離子從負極板引出的比例. 由(2)式—(4)式計算得到的共振電荷轉移的幾率列在表2的括號內, 其中通過鞘層的離子發(fā)生共振電荷轉移的幾率的計算結果均大于模擬結果, 這是因為靠近正極板的離子從正極板引出, 由負極板引出的離子在背景氣體中的平均路程長度小于背景氣體寬度的一半.

3.2 二維離子引出過程

本節(jié)計算平行板法、交替偏壓法、Π型電場法和M型電場法四種引出構型下的二維離子引出過程, 四種引出構型的示意圖如圖12所示. 計算參數(shù)如下: 左右收集板間距L= 5 cm, 高度H=6 cm; 尾料板距離左右收集板上沿為2 cm, 其長度為5 cm; 等離子體尺寸為Lp×Hp= 3 cm × 3 cm,位于左右收集板中心, 其初始密度分布為均勻分布;M型電場法中M極板的長度為5.0 cm, 即M極板下端恰好位于等離子體中心, 上端位于尾料板中心, M極板垂直于尾料板; 金屬蒸氣在極板間均勻分布, 其原子數(shù)密度為4.0 × 1011cm—3; 金屬蒸氣中A, B兩種同位素的豐度均為50%, 同位素A的電離率為2.5%, 同位素B的電離率為0, 即等離子體的初始離子密度為5.0 × 109cm—3. 其余計算參數(shù)同3.1節(jié). 平行板法、Π型法和M型法的引出電壓為2.0 kV, 交替偏壓法采用頻率為1.0 MHz、峰值為2.0 kV的正弦交流電壓. 計算時采用的網(wǎng)格數(shù)為250 × 400, 網(wǎng)格寬度為2 × 10—4m, PIC法和雜化PIC法的時間步長分別為100 ps和1 ns.

圖12 二維電場法離子引出示意圖 (a) 平行板電場法; (b) 交替偏壓法; (c) Π型電場法; (d) M型電場法Fig. 12. Schematic diagram of two dimensional electric ion extraction: (a) Parallel type; (b) alternately biased parallel type;(c) Π-type; (d) M-type.

圖14 交替偏壓電場法中, 收集板上的離子沉積能量分布 (a) 同位素A; (b) 同位素BFig. 14. Energy distribution of ions deposit on the collection plates in alternately biased parallel type: (a) Isotope A; (b) isotope B.

圖15 Π型電場法中, 收集板上的離子沉積能量分布 (a) 同位素A; (b) 同位素BFig. 15. Energy distribution of ions deposit on the collection plates in Π-type: (a) Isotope A; (b) isotope B.

圖13 —圖16分別為四種引出構型下極板上同位素A和B離子沉積能量分布曲線. 由于共振電荷轉移對離子引出時間的影響比較微弱, 本節(jié)計算得到的剩余離子/電子比例隨引出時間的變化曲線與圖3中基本相同, 因此不再贅述. 如圖13所示, 平行板電場法中左右極板上同位素A和B離子沉積能量分布與一維工況的模擬結果類似; 在交替偏壓法的模擬結果中, 由于極板上的電壓是周期性的交變電壓, 交變電壓的幅值是2.0 kV, 因此沉積能量的最大值明顯小于2.0 keV. 且有相當一部分離子從尾料板引出, 這是因為加交變電壓的極板有一半的時間是處于正電位, 此時離子從另外一側的極板和尾料板引出. Π型電場法和M型電場法中左右極板上的離子沉積能量分布與平行板電場法中負極板上的結果類似, 不同點在于Π型電場法中離子沉積能量的峰值在1.2 keV左右, M型電場法則為1.9 keV, 這是因為Π型電場法中等離子體電位小于正極板電位[19], M型電場法則擁有較快的離子引出速率. M型電場法中沉積在左右極板上的離子能量的最大值接近2.5 keV, 這是由于M型電場法中電子振蕩階段造成等離子體電位高于零電位, 因此在電勢回落階段引出離子的能量大于2.0 keV. 另外, M型電場法中有一部分離子沉積在M極板上, 其沉積能量的分布類似于平行板電場法中正極板上的離子能量分布.

圖13 平行板電場法中, 收集板上的離子沉積能量分布 (a) 同位素A; (b) 同位素BFig. 13. Energy distribution of ions deposit on the collection plates in parallel type: (a) Isotope A; (b) isotope B.

圖16 M型電場法中, 收集板上的離子沉積能量分布 (a) 同位素A; (b) 同位素BFig. 16. Energy distribution of ions deposit on the collection plates in M-type: (a) Isotope A; (b) isotope B.

表3給出了四種電場法離子引出構型下同位素A離子中與同位素B原子發(fā)生共振電荷轉移的比例. 可以看出, 在其他條件相同的情況下, M型電場法離子引出過程中離子發(fā)生共振電荷轉移的比例最小.

表3 四種引出構型中同位素A離子與同位素B原子發(fā)生共振電荷轉移的比例Table 3. Resonant charge exchange ratio between A-ion and B-atom in four electrode configurations above.

4 總 結

本文利用混合算法和PIC法計算了考慮共振電荷轉移的電場法離子引出過程, 結果表明: 離子在屏蔽層和鞘層中都會發(fā)生共振電荷轉移, 共振電荷轉移過程會略微延長離子引出時間, 發(fā)生共振電荷轉移的幾率主要取決于共振電荷轉移截面、引出時間、背景原子密度和蒸氣寬度. 在其他條件相同的情況下, M型電場法中離子發(fā)生共振電荷轉移的比例最小. 結合之前論文中的研究結果, 從縮短引出時間和降低離子發(fā)生共振電荷轉移比例的角度, M型電場法在本文介紹的四種離子引出構型中表現(xiàn)最優(yōu).