稀薄空氣中圓柱腔體內系統電磁脈沖的混合模擬*

張含天 周前紅 周海京 孫強 宋萌萌 董燁 楊薇 姚建生

(北京應用物理與計算數學研究所,北京 100094)

系統電磁脈沖廣泛存在于強電離輻射環境中,且難以有效屏蔽.為了評估稀薄空氣對系統電磁脈沖的影響,本文基于粒子-流體混合模擬方法,建立了三維非穩態模型,計算并分析了稀薄空氣等離子體的特性以及其與電磁場響應的相互作用.結果表明,壓力越高,光電子發射面附近的次級電子數密度越高,軸向分布的梯度越大,腔體中部的電子數密度在20 Torr (1 Torr=133 Pa)下出現峰值,而電子溫度隨壓力升高單調遞減.腔體內的稀薄空氣等離子體阻礙了空間電荷層的產生,電場響應峰值比真空條件下的低了一個數量級,電場脈沖寬度也顯著降低.光電子運動特性決定了電流響應的峰值,壓力升高,到達腔體末端的電流先增加再減小.而等離子體電流會抑制總電流的上升速率,并使電流響應出現拖尾.最后,將數值模擬結果與電子束模擬系統電磁脈沖的實驗結果進行比較,驗證了本文混合模擬模型的可靠性.本研究所采用的混合模擬方法相比于粒子云網格-蒙特卡羅碰撞方法,大幅減小了計算消耗.

1 引言

電離輻射(γ和X 射線)輻照下的金屬表面會通過光電效應等過程發射高能電子,高速運動的電子流在系統(如衛星、導彈以及各類電子儀器)中激勵出的瞬態電磁場響應被稱為系統電磁脈沖(system generated electromagnetic pulse,SGEMP).

特別地,又將在系統內部激勵出的SGEMP 稱為內電磁脈沖(internal electromagnetic pulse,IEMP)[1]或者腔體SGEMP.系統電磁脈沖難以有效屏蔽,會對電子設備產生干擾和損傷[2],有關SGEMP 損傷機制和防護方法的研究最早源自于高空核試驗,在其他強電離輻射環境,如高功率激光設備和加速器中,SGEMP 也是十分重要的干擾源[3].

為了降低實驗復雜性與成本,SGEMP 的實驗研究往往通過直接注入高能電子束來模擬光電子,避免使用大型射線源,即便如此,相關實驗的難度依然很高,公開的實驗數據也非常稀缺[4-6].因此,數值模擬已經成為研究SGEMP 的重要手段,代表性成果是美國Woods 等[7]編制的二維模擬程序ABORC.以往的研究中,真空環境下的SGEMP一直是關注的重點[8].近五年左右,西北核技術研究所的Wang 等[9-12]基于自研程序UNIPIC-3D對SGEMP 進行模擬,該程序支持辛算法以及共形網格,可以有效地處理系統電磁脈沖中的電子發射邊界問題.國內研究人員對其中存在的空間電荷限制效應[13,14]、壁面二次電子發射[15]等問題也進行了廣泛的研究[16].

然而,實際系統往往工作在氣體環境或含有充氣元器件,光電子與中性氣體間會通過碰撞電離產生大量的次級電子、離子對,次級電子又會從電磁場中獲得能量,繼續與中性氣體發生電離、激發等一系列反應,從而形成等離子體.實驗結果表明,等離子體與電磁場的相互作用會使SGEMP 的特性會發生顯著變化[17].Woods 等[18]最早通過數值模擬研究了高空稀薄空氣和預電離空氣對SGEMP的影響,結果表明,隨著氣壓升高,空間電荷層被破壞,系統表面電流增加[19].Chan 和Woods[20]以及Strasburg 等[21]使用電導率模型描述稀薄空氣等離子體.Pusateri 等[22]基于swarm 模型對高空核爆電磁脈沖中的二次電子進行建模.Angus 等[23]建立了零維反應動力學模型,研究電子束誘導稀薄空氣等離子體中的主要反應.Ribière 等[24]基于粒子云網格方法(particle in cell,PIC),研究0.9 MeV脈沖 X 射線在充氣(3—570 Pa)金屬腔內誘發等離子體的特性,在該模型中,二次電子和離子并不自恰生成,而是在空間中均勻加載,在時間上的變化服從預定義分布函數.在之前的研究中[25],作者已經使用PIC-MCC(粒子云網格-蒙特卡羅碰撞)方法研究了壓力0—500 mTorr (1 Torr=133 Pa)空氣等離子體與光電子-等離子體的相互作用,但對于更高壓力條件,該方法計算耗時過長,不具備實用性.為了更好地評估稀薄空氣條件下X 射線引起的系統電磁脈沖,本文基于PIC-流體混合模擬方法,建立了三維非穩態SGEMP 數值模擬模型,研究稀薄空氣等離子體的特性以及其對電磁場響應的影響.

2 數值模型

本文建立的模型關注高能光電子本身的運動、高能光電子與中性氣體的相互作用以及在該過程中激勵出的電磁場響應.主要計算流程如圖1 所示,分別使用PIC 模型描述高能光電子的運動,使用流體模型描述中性氣體電離產生的次級電子特性.高能光電子與次級電子共同貢獻電流密度,影響電磁場,電磁場又通過電磁力改變光電子和次級電子的運動狀態,從而實現耦合模擬.

圖1 粒子-流體混合模擬流程圖Fig.1.Flow chart of the hybrid particle-fluid model.

需要指出的是,在稀薄氣體環境中(次級電子的碰撞自由程與腔體尺寸相當),蒙特卡羅碰撞(Monte Carlo collision,MCC)方法也是一種有效地描述光電子-中性氣體碰撞電離產生次級電子-離子對,以及隨后發生雪崩電離過程的建模方法[26,27],該方法對次級電子、離子也使用粒子進行描述(即全粒子方法),并通過碰撞截面獲得各類碰撞發生的概率.然而,在使用PIC-MCC 方法進行計算時,設置的網格尺寸Δx和時間步Δt需要同時滿足如下條件:

其中(1)式為Courant-Friedrich-Lewy(CFL)條件,c為光速.(2)式中的ωp為等離子體頻率.(3)式表示網格尺寸需要滿足精細網格不穩定性條件(fine grid instability),ξ為數量級為1 的參數,受插值函數以及電子能量分布影響[28],為德拜長度,Ne,Te為電子數密度和電子溫度,ε0,kB,e分別為真空介電常數、玻爾茲曼常數以及電子電量.若網格尺寸不滿足(3)式,則會造成非物理的“自加熱”現象,并導致計算逐漸發散[29,30].然而,較高氣壓下,產生的等離子體電子溫度低且密度高,如當電子數密度1020m—3,PIC-MCC 方法中就需要令Δx≈10—6m,Δt≈10—15s,模擬的計算開銷很大.在作者之前的研究中[25],使用PIC-MCC方法計算壓力數百毫托條件下的腔體SGEMP 過程,四核并行需要超過20 h.對于更高壓力下的SGEMP 問題,即便是使用粒子合并等加速手段,其三維建模在現有計算能力下也難以實現.因此,擬采用流體方法模擬等離子體,相比于全粒子的PIC-MCC 方法,可以大幅減小計算量,如本文涉及算例的典型計算耗時約為1 h.

2.1 光電子的電磁粒子模擬方法

電磁場部分的建模,采用時域有限差分法求解,利用(4)式所示的Maxwell 方程的兩個旋度方程.并使用Esirkepov 的密度分解法[31]以滿足電荷連續性方程?·J=-?ρ/?t.

其中E,B,H,D和J分別為電場強度、磁感應強度、磁場強度、電位移矢量和電流密度.電磁場在空間上基于Yee 元胞離散化,即將E定義在網格的棱邊上,B定義在網格的面心上[32].在時間上,上述電磁場的旋度方程使用如(5)式的蛙跳格式推進,上標為推進的時間步[29]:

高能光電子與其誘導產生的等離子體共同貢獻電流密度J=Jp+Js,從而實現與電磁場的耦合.其中,Jp為高能光電子的電流密度,該電流密度通過在PIC 方法中統計光電子的運動獲得,而Js為次級電子電流項,求解過程見2.2 節.

在采用PIC 方法描述高能光電子的運動時,通常不會跟蹤所有實際電子,而是用“宏粒子”代表若干相鄰電子的運動,本文中一個“宏電子”代表106個實際電子.高能光電子的位移x和速度v隨時間t的變化服從(6)式和(7)式.

其中Ztar*與J*分別為修正的靶原子序數、等效平均激發能量,修正依據為Bohr 電子剝離判據[35,36],入射粒子能量越低,靶的Ztar*與J*都越小,數值上可用多項式擬合的經驗公式表達.zinc*為入射電子的等效電荷數,F(τ) 為相對論修正因子,詳見文獻[37].(8)式為修正形式的Bethe 公式,可用于計算介質對中等能量電子(數10 eV—10 keV)的阻止本領.圖2 為由(8)式得到的N2和O2對不同能量電子的質量約化阻止本領[38],算得結果的高能部分(＞104eV)與NIST 數據庫[39]一致.

圖2 N2 和O2 對電子的約化阻止本領(左)[38];電子-N2 的電離碰撞截面(右)Fig.2.Mass stopping power for incident electron in N2 and O2 (left)[38];ionization cross sections between electrons and N2 (right).

穿過介質時,光電子的運動方向也發生變化.本文用多重散射模型處理光電子運動方向的偏轉[40,41],在更新粒子動量方程(7)式后,將粒子運動方向由p=pep變為p′=pcosθep+psinθ(cosφe1+sinφe2),如圖2 右上所示,散射在方位角方向均勻分布φ=2πRφ,Rφ為0 到1 之間均勻的隨機數.散射的極角〈θ〉=σRθ,Rθ為滿足正態分布的隨機數,,均方散射角 〈θ2〉的變化率[40]滿足(9)式,

其中Z=7.2 為空氣平均原子序數;Natom為背景氣體中原子的數密度;p=meγv為光電子動量[42];a為原子屏蔽半徑,根據Thomas-Fermi 模型,可取a=a0/Z1/3,a0為Bohr 半徑[43];? 為約化普朗克常數.

2.2 等離子體的模擬方法

光電子電離中性氣體以及次級電子雪崩電離共同導致了系統中等離子體的產生.對于等離子體部分的建模[44],本模型使用流體方法,求解次級電子的數密度Ne、漂移速度vd、電子能量U=kBTe,以獲得次級電子的電流密度Js=eNevd.

(10)式為次級電子數密度Ne守恒方程.Spe=Npe(vpeNairσi)=JpNairσi/e為光電子與中性氣體碰撞產生次級電子的速率.不同于冷態氣體的電擊穿過程,在X 射線整個脈沖時間范圍內,高能光電子會持續注入系統.SGEMP 發展的最初階段,光電子電離產生次級電子的過程占主導,其中Nair,Npe分別為空氣數密度以及入射光電子數密度,vpe為光電子的運動速率,σi為電離碰撞的碰撞截面,該截面與靶分子種類相關,通過電子能量查表插值獲得[45-47].圖2(右側坐標)給出了電子-N2的電離碰撞截面,電子能量小于2 keV 的數據來自Itikawa[45],高能部分數據通過RBEB(relativistic binary-encounter-Bethe)模型[48]計算.右側第二項為次級電子電離反應導致的電子數增加,ki與katt分別為電離、吸附率.因為程序存在除Ne的運算,因此,初始次級電子密度設置為105m—3,經過測試,該初始值不會影響計算結果.

由于電流密度較小,次級電子動量守恒方程(11)式中忽略了磁場力evd×B,νm為動量轉移碰撞頻率.能量守恒方程(12)式右側第一、第二項分別為次級電子從電場中獲得的能量、通過與環境氣體碰撞損失的能量(U0為環境氣體溫度),νw為能量轉移碰撞頻率.右側第三項表示次級電子與背景氣體電離(即雪崩電離)造成的能量損失,wc為空氣的平均電離能,約為14 eV[22].右側第四項pe相當于光電子電離反應給等離子體注入的能量.為次級電子的平均初始能量,確定該值準確大小較為困難,且不同文獻里的取法也存在差異[44],這里簡要地給出估算過程.對于電離反應,光電子的能量損耗等于εinc-εsc=εiz+εse,εinc,εsc為碰撞前后光電子的能量,εiz為電離能,εse為次級電子能量.(8)式所述Bethe 公式給出了總的能量損失(包括電離、激發等各類碰撞過程),又由于電離為主要的耗能反應,所以可利用阻止本領估計次級電子平均能量的上限,即,算得的次級電子平均能量如圖3 所示.

圖3 次級電子平均能量Fig.3.Mean energy of the secondary electrons.

Longmire 和Longley[40]用類似的方法估算了高能電子與空氣相互作用的能量損失(圖3 中黑色離散點),與本文的估算結果比較接近,差別是因為使用了不同的電離碰撞截面數據.需要指出的是,一般計算中認為高能電子每消耗34 eV 能量產生一對電子-離子對,該值實際上是對整個電離過程平均后的結果,與本文的估算值并不矛盾.由圖3可知,能量在MeV 量級的電子約消耗86 eV(加上電離能)產生第一對電子-離子對,第一次級電子的能量較高,會繼續電離中性氣體,逐漸產生其它的次級電子-離子對,綜合下來,相當于高能電子每消耗34 eV 能量電離出一個次級電子[49].

(10)—(12)式的計算需要額外使用三個系數,即動量轉移碰撞頻率νm、能量轉移碰撞頻率νw和電離率ki,統稱為swarm 參數,所以也稱本節用于描述等離子體的流體模型為swarm 模型[50].本文使用BOLSIG+軟件[51]求解二階玻爾茲曼方程來獲得上述三個參數,并將其擬合為電子溫度U的函數用于后續計算[22].計算考慮了彈性碰撞、電離、振動、轉動激發等多類碰撞過程(含e-N2的25 個反應和e-O2的17 個反應),詳見文獻[52].計算得到的密度約化swarm 參數結果如圖4 所示,該結果與Higgins 等[44]總結的實驗結果具有較好的一致性.

圖4 約化swarm 參數隨電子能量的變化Fig.4.Normalized swarm parameters as a function of electron energy.

2.3 計算區域和邊界條件

圖5 為本文涉及的計算模型,所模擬的金屬腔體為空心圓柱,內直徑2R=100 mm,高度H=50 mm.總的來說,SGEMP 的模擬可以分為兩個階段.首先,X 射線光子與金屬靶材相互作用使得金屬表面發射出一系列光電子,光電子的產額、能譜和角分布等特性目前主要借助蒙特卡羅方法計算光子-電子輸運獲得[53,54].第二階段,將光電子作為激勵源,計算系統中的電場和磁場響應.數值模擬中,這兩個階段之間為弱耦合,即第二階段的總光電子數目由第一階段計算,且假設發射光電子的電流密度與X 射線注量具有相同的時間分布[8].本模型主要關注由光電子激勵出電磁場響應這一過程,即不包含第一階段X 射線與金屬壁面的相互作用,直接假設左側金屬端面(x=0,半徑r≤25 mm區域)發射單向(沿x軸向右)、單能(20 keV)的光電子[55],光電子電流密度滿足如下脈沖函數:

圖5 計算模型Fig.5.Schematic of the calculation domain.

其中脈沖的半高寬τ=2 ns.由發射電流I=πr2Jemi=eΔnemi/Δt,可獲得每一時間步注入的光電子數量Δnemi.本文模型用于研究氣體環境下腔體SGEMP 的基本特性,因此對光電子的發射特性進行了簡化.實際情況中光電子具有更為復雜的時空分布以及能譜分布,相應地只需要在建模時修改描述電子發射的參數,如在3.2 節中就對一電子束模擬SGEMP 實驗中的電磁場響應進行了建模.

3 結果與討論

3.1 腔體內壓力對系統電磁脈沖的影響

圖6 給出了不同壓力(0—20 Torr)條件下,金屬腔體內部光電子以及次級電子的空間分布隨時間的演化過程.圖6 第一行為真空條件下的光電子分布,圖中每個點代表一個宏粒子,如2.3 節所述,光電子初始運動方向自左往右.1 ns 時,最右側光電子區域所占半徑相比于最左側基本不發生變化.脈沖開始后2 ns,光電子在距發射面約20 mm處發生聚集,形成了高密度的光電子層,右側光電子區域的空間尺寸變大,半徑由x=0 處的25 mm增加至x=50 mm 處的約42 mm.3 ns 時,腔體內部空間大致被分為兩個部分,右側光電子的密度要顯著低于左側.靠近金屬腔體側壁的光電子具有更快的運動速度,高密度光電子層在空間上呈“漏斗狀”,厚度約為5 mm,隨后,高密度光電子層逐漸整體向右運動.

圖6 光電子與次級電子的分布(0—20 Torr,1—4 ns)Fig.6.Distributions of photoelectrons (red dot) and secondary electrons (0—20 Torr,1—4 ns).

相比于真空條件,腔體內空氣壓力5 Torr 條件下,光電子的演變出現了顯著的區別.1—4 ns,光電子區域在徑向方向上沒有出現擴張,大部分光電子在到達右側端面后仍在r＜ 25 mm 的范圍內,且沒有出現光電子高密度薄層.氣壓升高至10 Torr,光電子在右側區域再次出現了擴張,20 Torr 條件下,徑向方向上的擴張更為顯著,但均沒有出現真空條件下的光電子薄層.

圖6 還給出了不同時刻,10 Torr 和20 Torr條件下次級電子的空間分布.由于最初的次級電子是由光電子電離產生,所以2—3 ns 次級電子在空間上的分布大致與光電子相同,即腔體左側電子數密度高、空間半徑小,右側次級電子較分散,電子密度也低.次級電子高密度區域逐漸向右側發展,其變化滯后于光電子的變化,4 ns 時,光電子脈沖結束,但左側區域次級電子的密度仍高于右側.次級電子的數密度隨腔體內壓力增加而增大,10 Torr、4 ns 時,次級電子的最大數密度約為2.5 × 1018m—3,20 Torr 條件下,最大數密度約為5 × 1018m—3.光電子的峰值數密度約為2 × 1015m—3,在光電子入射后很短時間內,系統中次級電子數量就占主導.

圖7(a)給出了真空條件下,不同軸向位置上的電場強度Ex隨時間的變化.x=2 mm 處的電場約在2.5 ns 達到峰值4.8 × 106V/m,如此高電場的產生主要是因為空間電荷效應(space-charge limited effect).由于電荷守恒,金屬在X 射線的輻照下發射出光電子,金屬器壁本身隨即帶等量的正電荷,光電子發射面與光電子之間產生電場,且該電場的方向與光電子初始運動方向相同.隨著發射出光電子電流密度的升高,電場強度快速增強,在電場力作用下,光電子速度顯著降低并在發射面附近聚集,形成了所謂的空間電荷層(space charge barrier layer),也就是圖6 中的次級電子高密度區域.在空間電荷層內,光電子被減速,甚至被反向加速返回到發射面,發射面同時又不斷向空間電荷層內補充電子,維持空間電荷層的動態平衡.該空間電荷層的存在也會阻礙后續光電子向右運動,導致圖6 中真空條件下、3 ns 時右側區域具有較低的電子密度.2 ns 后,由于注入光電子的電流密度逐漸減小,空間電荷層內的電場強度達到峰值后無法繼續維持,造成了圖7(a)中x=2 mm 處電場的快速降低,空間電荷層也逐漸向右移動.圖6 中金屬側壁附近空間電荷層的運動更快,主要是由于側壁附近較低的光電子密度,造成了相對更弱的空間電荷限制效應.圖7(a)中x=10 mm 處的電場在2.5 ns 附近出現了負值,這是因為此時光電子的電流密度接近峰值,空間電荷效應增強,高負電量的空間電荷層與光電子發射面之間的距離縮短,空間電荷層運動到了x=10 mm 的左側位置.

圖7(b)給出了10 Torr 壓力下,軸向電場強度隨時間的變化.x=2 mm 處的電場約在0.8 ns達到峰值2.8×105V/m,僅為真空條件下電場峰值的1/17,電場脈寬(以下均為半高全寬)也由真空下的1.7 ns 減小到約0.8 ns.x=10 mm 處電場峰值稍滯后于x=2 mm 處,與40 mm 處電場的波形相似,但極性相反.可以發現,在10 Torr 壓力下,氣體的存在顯著削弱了腔體SGEMP 中電場響應的峰值以及脈沖寬度,這主要是由于次級電子的產生以及增殖,削弱了空間電荷效應.值得注意的是,x=2 mm 與x=10 mm 處的軸向電場在約2.5 ns 后出負值,并在3 ns 左右達到負峰值,約3 × 104V/m,隨后又逐漸衰減到0.由圖6 可知,次級電子的空間分布并不均勻,靠腔體左側的密度更高,軸向方向次級電子數密度的梯度導致了沿x軸負方向的電場.

圖7 真空(a)以及10 Torr 壓力(b)下,軸向電場Ex 隨時間的變化(x=2,10,40 mm;r=0 mm)Fig.7.Time-dependent electric field Ex (x=2,10,40 mm;r=0 mm) for (a) vacuum and (b) 10 Torr.

圖8 給出了靠近腔體右側端面(x=48 mm,r=48 mm)的切向磁感應強度Bφ隨時間的變化.真空條件下,Bφ約在1.9 ns 到達負峰值—5 × 10—4T,脈寬為1.1 ns.10 Torr 壓力下,磁感應強度峰值(—7.9×10—4T)約為真空條件下的1.6 倍,脈沖寬度增加至2.4 ns.但進一步提高腔體壓力,磁感應強度的幅值會發生降低,20 Torr 下Bφ的峰值已經與真空條件下Bφ相當,50 Torr 對應的峰值Bφ僅為真空下的18%,但脈沖持續時間均大于4 ns.真空條件下,切向磁感應強度主要是由光電子的定向運動產生的,2 ns 后Bφ出現快速下降,是因為近發射面空間電荷層的存在阻礙了光電子向右側端面運動,如圖6 所示.

圖8 不同壓力下,切向磁感應強度Bφ 隨時間的變化(x=48 mm;r=48 mm)Fig.8.Time-dependent magnetic flux density Bφ at different pressures (x=48 mm;r=48 mm).

非真空環境下,光電子電流以及次級電子電流均對磁感應強度有貢獻,圖9 分別給出了軸向不同位置上的光電子電流以及次級電子電流,該電流是通過對垂直于x軸的截面上的電流密度Jp以及Js求面積分獲得.光電子電流在模擬時間范圍內均為負.0—2.5 ns,通過x=10 mm 截面的次級電子電流與光電子電流的方向相反,這是由于初始階段光電子占主導,腔體中的光電子與發射面之間會建立較強的電場,在這一正向電場力的作用下,次級電子漂移速度的方向與光電子運動方向相反,次級電子電流反過來又抑制總電流的快速上升,總電流出現峰值的時刻相比光電子電流峰值時刻滯后了約0.5 ns.在次級電子增殖過程中,軸向方向上存在電子密度梯度,產生了負電場(如圖7(b)所示),次級電子電流也隨之反向.在入射光電子脈沖的靠后階段,約3.5 ns 后,次級電子電流開始占主導,并達到負峰值,約—85 A.同時,次級電子通過(11)式中的動量交換項減速,次級電子電流逐漸減小為零.次級電子電流的存在導致了總電流在入射光電子脈沖結束后出現“拖尾”,通過x=10 mm 截面的總電流在8 ns 后才接近0.在光電子入射后的最初2 ns,由于電場為負,圖9(b)中x=40 mm 處的次級電子電流與光電子電流同向,隨后此處的總電流也同樣出現了拖尾.

圖9 截面上光電子、次級電子電流隨時間的變化(10 Torr) (a) x=10 mm;(b) x=40 mmFig.9.Time-dependent current of photoelectrons and secondary electrons (10 Torr) for the plane:(a) x=10 mm;(b) x=40 mm.

由圖9 可知,當前計算條件下,次級電子電流會影響總電流峰值出現的時刻并造成電流響應的脈寬增大,但總電流的峰值主要還是由入射光電子的運動決定.因此,圖8 中磁感應強度峰值在10 Torr 以上減小,主要是因為更高壓力條件下,背景中性氣體對光電子的散射作用增強,更多的光電子在運動過程中被散射,能到達右側端面的光電子數降低.

圖10 給出了2 ns 時對稱軸上的次級電子數密度Ne分布,從左往右次級電子數密度大致是遞減的.腔體左側(x＜ 15 mm),次級電子數密度隨壓力增加而增加,靠近左側端面的Ne由5 Torr 下的1.5 × 1018m—3增大為50 Torr 下的5.0 × 1018m—3,主要是因為(10)式中的光電子電離、次級電子雪崩電離源項均正比于背景中性氣體數密度,高背景氣體壓力有利于電離反應進行.然而,50 Torr 下,Ne在軸向方向上的梯度更大,沿x軸的衰減更快,在x＞ 15 mm 已經低于20 Torr 壓力下的Ne.

圖10 不同壓力下,對稱軸上次級電子數密度的分布(r=0 mm,t=2 ns)Fig.10.Distributions of secondary electrons along the axis at different pressures (r=10 mm,t=2 ns).

圖11 進一步給出了x=24 mm 處次級電子數密度沿徑向的線積分隨時間的變化.不同壓力下,腔體中部的電子數密度先快速上升,入射光電子脈沖結束后,Ne逐漸趨于飽和.Maulois 等[56]指出,在數十Torr 壓力條件下,O2的解離吸附等使電子密度衰減的反應約在數百納秒后才占主導,腔體內次級電子、離子數密度逐漸降低,并恢復成初始中性氣體的整個過程可能要持續數個微秒.5—20 Torr,Ne線積分的最大值隨壓力增大而增大.0—4 ns 內,次級電子密度的平均上升率由5 Torr 的1.98 × 1017m—3/ns 變為20 Torr下的5.84 × 1017m—3/ns.然而,50 Torr 下的Ne線積分略低于20 Torr 下的結果,這主要是因為高氣壓下光電子散射更嚴重,參與后續雪崩電離的種子電子少,次級電子增長發生了滯后.雖然圖10中x=24 mm,r=0 mm 處50 Torr 的次級電子數密度已經略低于10 Torr 下的Ne,但50 Torr下Ne徑向分布范圍更大,所以其徑向線積分仍高于10 Torr 下的結果.文獻[5,57]中也報道過類似的次級電子數密度隨氣壓上升先增加,再減少的變化關系.

圖11 不同壓力下,x=24 mm 處次級電子數密度的徑向線積分隨時間的變化Fig.11.Line-integrated secondary electrons at x=24 mm at different pressures.

圖12 為對稱軸上次級電子的溫度分布.脈沖開始后2 ns,5—50 Torr 條件下的電子溫度在2—8 eV,且隨著壓力增加,電子溫度單調遞減,這是因為等離子體能量守恒方程(12)式中的能量轉移碰撞頻率正比于背景氣體數密度,壓力越高,次級電子通過碰撞損失的能量也就越大.而同一壓力條件下,靠近右側端面的電子溫度更高,這主要是受雪崩電離過程影響.次級電子從電場中獲得的能量,除了通過與環境氣體碰撞產生損失外,還會通過雪崩電離反應被損耗.2 ns 時,腔體右側通過雪崩電離產生的次級電子較少(如圖10 所示),因此電子能量損失較低.圖12 還給出了5 Torr,不同時刻的電子溫度分布,2 ns 后,對稱軸上的電子溫度隨著時間增加而降低,又因為腔體右側的次級電子數量也逐漸提高(如圖6 所示),軸向的溫度梯度發生了減弱.

圖12 對稱軸上的次級電子溫度分布(r=0 mm)Fig.12.Distributions of secondary electron temperature along the axis (r=0 mm).

3.2 與實驗比較

利用射線源對關注的系統開展模擬輻照實驗是研究SGEMP 最直接的方法,但此類實驗成本高昂、準備周期長,實施難度大,特別是高注量、大空間范圍的射線源難以獲得,所以往往直接用電子束模擬發射面光電子的特性進行實驗研究.20 世紀70 年代,美國陸軍哈瑞戴蒙德實驗室進行了一系列電子束模擬不同條件下腔體SGEMP 效應的實驗[17],本節參考其中短脈沖實驗條件開展數值模擬工作,并將模擬結果與已公開的實驗結果進行比較,以進一步驗證本文混合模擬模型的準確性.實驗腔體與圖5 類似,尺寸略大,為高150 mm、直徑300 mm 的圓柱腔體,電子束從左側端面射入.模擬中假設電流密度在整個發射面上均勻,圖13為電子束的電流以及發射電子的平均能量隨時間的變化,均來自于實測數據[17].受測量條件限制,實驗文獻無法給出發射電子的能譜以及角分布,這里參考Chan 和Woods[20]的處理方法,假設每一時刻發射電子的能量均等于平均能量,電子能量高于25 keV 單向發射(沿x軸向右),低于25 keV 的光電子以cos2θ的角分布概率發射,θ為與x軸的夾角.

圖13 實驗中電子束電流(左)與電子能量(右)隨時間的關系Fig.13.Electron beam current (left) and energy (right) as functions of time.

圖14 給出了腔體內壓力為3 Torr 時,Z=12.5 cm,r=13.8 cm 處的磁場強度隨時間的變化.模擬結果的趨勢與實驗結果基本一致,混合模擬模型得到的磁場強度峰值比實驗結果高了約10%.造成差別的原因可能是發射電子角分布的假設不夠準確,圖14 中也給出了假設所有電子均以cos2θ的角分布發射的模擬結果,該假設條件下,磁場強度峰值出現降低.也有可能是因為計算初始的數個納秒,次級電子密度較低,用基于流體的swarm 模型描述等離子體存在偏差.但總的來說,本文使用的混合模型已經可以用較低的計算資源獲得腔體SGEMP 的動態響應.

圖14 Z=12.5 cm,r=13.8 cm 處的磁場強度(腔體內壓力3 Torr)Fig.14.Magnetic field at Z=12.5 cm,r=13.8 cm for an air pressure of 3 Torr.

4 結論

為了研究背景氣體對系統電磁脈沖特性的影響,本文考慮了稀薄空氣對光電子的減速、散射以及等離子體與電磁場響應的相互作用,實現了基于PIC-流體混合模擬方法的三維非穩態SGEMP數值模擬.計算了0—50 Torr 壓力下的圓柱腔體內系統電磁脈沖,主要結論如下.

1) 稀薄空氣等離子體削弱了光電子發射面附近的空間電荷效應,阻礙了空間電荷層的產生,電場響應峰值比真空條件下的低了一個數量級,電場脈沖寬度也顯著降低.

2) 光電子運動特性決定了電流響應的峰值,壓力升高,腔體末端磁感應強度先增加再減小.但等離子體電流會抑制總電流的上升速率,并造成脈寬增大,電流響應出現拖尾.

3) 壓力越高,發射面附近的次級電子數密度越高,軸向分布的梯度越大,腔體中部的電子數密度在20 Torr 下出現峰值.電子溫度隨壓力的增高單調遞減.

4) 數值模擬結果與電子束模擬SGEMP 實驗結果進行比較,驗證了本文混合模擬模型的可靠性.