陳再高

(強脈沖輻射環(huán)境模擬與效應全國重點實驗室;西北核技術(shù)研究所: 西安 710024)

全電磁粒子模擬方法通過求解帶電粒子與電磁場之間的耦合方程,能準確描述電磁場隨時間的演化過程及帶電粒子的運動參數(shù),在等離子體相關(guān)物理現(xiàn)象研究[1]、真空電子器件的設計與優(yōu)化[2-3]、介質(zhì)擊穿模擬[4]及核爆炸人工輻射帶模擬[5]等方面得到了廣泛應用。此外,高能射線作用于飛行器等電子設備,在電子系統(tǒng)內(nèi)外激發(fā)光電子或康普頓電子,由于電荷分布不均勻,飛行器表面會產(chǎn)生較大的置換電流,最終在系統(tǒng)周圍,甚至內(nèi)部激勵強電磁脈沖,全電磁粒子模擬方法也能用于模擬高能射線產(chǎn)生系統(tǒng)電磁脈沖等物理現(xiàn)象[6-7]。國內(nèi)外均開展了全電磁粒子模擬算法的研究及軟件的研制,其中,國外研究機構(gòu)和企業(yè)研制了MAGIC[8]、KARAT[9]和CST[10]等全電磁粒子模擬軟件,國內(nèi)西北核技術(shù)研究所和西安交通大學聯(lián)合研制了2.5維粒子模擬軟件UNIPIC[11]和3維并行粒子模擬軟件UNIPIC-3D[12],電子科技大學研制了3維粒子模擬軟件CHIPIC[13],中國工程物理研究院研制了粒子模擬軟件NEPTUNE[14],這些軟件在高功率微波源器件的設計及參數(shù)優(yōu)化方面發(fā)揮了重要作用。

在全電磁粒子模擬軟件中,電磁場方程一般采用時域有限差分 (finite difference time domain,FDTD)方 法,該方法的計算量與網(wǎng)格數(shù)呈線性關(guān)系,具有程序?qū)崿F(xiàn)相對簡單及容易并行等優(yōu)點。FDTD方法采用正交結(jié)構(gòu)網(wǎng)格。傳統(tǒng)的階梯網(wǎng)格處理復雜邊界會帶來偏差,而共形網(wǎng)格能精確描述復雜器件結(jié)構(gòu)的邊界,基于共形網(wǎng)格的FDTD方法能大幅提高復雜結(jié)構(gòu)的計算精度[15]。因此,在全電磁粒子模擬算法中引入計算機輔助設計技術(shù),實現(xiàn)共形網(wǎng)格的自動化生成,能大幅提高全電磁粒子模擬算法計算復雜結(jié)構(gòu)的精度。此外,將全局優(yōu)化算法與全電磁粒子模擬算法結(jié)合,能實現(xiàn)真空電子器件幾何參數(shù)及電學參數(shù)的全局優(yōu)化[16-17]。

最新發(fā)展的時域間斷有限元方法在Maxwell方程的弱解積分形式中引入數(shù)值通量,通過方程離散獲得非結(jié)構(gòu)網(wǎng)格單元內(nèi)電磁場的推進公式,每個單元中的電磁場值只與相鄰單元有關(guān),不用求解全域矩陣的逆。時域間斷有限元算法能采用非結(jié)構(gòu)網(wǎng)格描述復雜結(jié)構(gòu)邊界,同時具備FDTD方法顯式求解方程的優(yōu)點,將間斷有限元引入到全電磁粒子模擬算法,能大幅提高復雜幾何模型的處理能力,同時,在電磁場的時域求解方程中引入了電勢矯正方程,能減少時域間斷有限元全電磁粒子模擬方法中電荷不守恒帶來的偏差,這是粒子模擬方法研究的熱點。

本文介紹了作者所在團隊在全電磁粒子模擬方法研究及軟件研制方面取得的進展,通過20多年的技術(shù)研究及軟件研制,研究了基于計算機輔助設計的共形結(jié)構(gòu)網(wǎng)格生成技術(shù),研制了2.5維全電磁粒子模擬軟件UNIPIC及3維全電磁粒子模擬軟件UNIPIC-3D,基于并行遺傳算法及2.5維粒子模擬算法研制了真空電子器件優(yōu)化設計軟件,開發(fā)了基于時域間斷有限元的2.5維全電磁粒子模擬程序。

1 2.5維全電磁粒子模擬技術(shù)

2.5維粒子模擬軟件UNIPIC軟件主要用于模擬真空電子器件內(nèi)部的電磁波與相對論電子束之間的非線性相互作用,包括非相對論與相對論范疇的束-波互作用物理現(xiàn)象[11]。UNIPIC軟件可模擬的2維問題包括x-y,z-r,r-φ3種坐標系;能設置理想導體、有耗導體、注入波及吸收波邊界等邊界條件,可通過熱發(fā)射、場致發(fā)射、二次電子發(fā)射、爆炸發(fā)射及光電發(fā)射等方式產(chǎn)生帶電粒子。另外,UNIPIC軟件還能模擬帶電粒子與中性氣體的彈性、電離、激發(fā)、電荷交換及復合等碰撞現(xiàn)象[18-19]。本文作者所在團隊將UNIPIC軟件的核心計算模塊與并行遺傳算法相結(jié)合,研制了真空電子器件并行優(yōu)化軟件。

1.1 UNIPIC軟件的數(shù)值算法

圖1為全電磁粒子模擬核心算法流程框圖。UNIPIC軟件中電磁場由Maxwell方程描述,帶電粒子運動由Newton-Lorentz方程描述,Maxwell方程與Newton-Lorentz方程之間的耦合采用CIC(cloud-in-cell)插值算法。采用FDTD方法求解電磁場隨時間的演化過程,每隔一定時間步用Maxwell方程的散度方程進行偏差校正。FDTD方法的時間步長受穩(wěn)定性條件的限制,在算法上選用中心差分的顯式格式,通過加入阻尼因子實現(xiàn)時偏格式,進一步減少電磁場計算帶來的數(shù)值噪聲。初始電磁場可采用輸入?yún)?shù)或輸入文件的方式加載,也可通過求解Poisson方程獲得。UNIPIC軟件具有直角坐標系、柱坐標系及極坐標系下2維Poisson方程求解功能。

在2維模擬中,粒子由2維坐標和3維速度表示,采用Boris提出的“半加速-旋轉(zhuǎn)-半加速”方法求解。依據(jù)電磁場的穩(wěn)定性條件,粒子的推進距離不超過一個網(wǎng)格,網(wǎng)格邊界處粒子在一個時間步長內(nèi)可能會出現(xiàn)穿過多個網(wǎng)格的情況,這時需要根據(jù)粒子運動軌跡與網(wǎng)格線的交點,分段計算粒子運動產(chǎn)生的電流。當帶電粒子運動較慢時,采用多次電磁場推進/一次粒子推進的方式提高算法的計算效率。UNIPIC軟件還具有加載任意分布初始粒子的功能。

由于電磁場及電流分布在離散的網(wǎng)格點或線上,而粒子位置則位于計算區(qū)域內(nèi)任意位置,粒子推進時產(chǎn)生的電流需要分配到網(wǎng)格點或線上,推進粒子運動時粒子所在位置處的電磁場值需要通過網(wǎng)格點上的電磁場插值獲得。為避免出現(xiàn)粒子的自推進現(xiàn)象,同時控制插值所需計算量,UNIPIC軟件實現(xiàn)了一階最近網(wǎng)格點插值及二階雙線性插值2種方式。由于插值會造成一定的數(shù)值偏差或噪聲,軟件對數(shù)值噪聲進行濾波處理,減少長時間計算所帶來的累積偏差。

圖1 全電磁粒子模擬核心算法流程框圖Fig.1 Flow chart of fully electromagnetic particle simulation method

1.2 粒子模擬中碰撞電離數(shù)值模擬技術(shù)

在氣體放電器件和填充等離子體器件中,氣體分子、原子和帶電粒子之間會發(fā)生碰撞,一般采用蒙特卡羅碰撞(Monte Carlo collision,MCC)模型來處理[18]。MCC是采用碰撞時間隨機的方法,而PIC方法中粒子推進和場推進的時間步長是固定的,因此采用在一個時間步長內(nèi)隨機決定粒子間是否發(fā)生碰撞的方法來實現(xiàn)蒙特卡羅碰撞。

碰撞模型中的粒子分為靶粒子和源粒子。假設第i個源粒子和靶粒子之間可發(fā)生N種碰撞,設第j種(j∈[1,N])碰撞的碰撞截面為σj(εi),總碰撞截面為所有類型的碰撞截面之和,表示為

(1)

第i個源粒子與靶粒子發(fā)生碰撞(包括所有N種碰撞)的概率,即碰撞總概率可表示為

Pi=1-exp[-σT(εi)nt(xi)viΔt]

(2)

其中:xi,vi分別為第i個源粒子的位置和速度;nt(xi)為第i個源粒子位置處靶粒子的數(shù)密度;Δt為PIC推進的時間步長。

為避免在每次推進中對每個源粒子進行大量計算,在PIC-MCC方法中引入空碰撞概念。碰撞總概率可設定為一個常數(shù)Pnull,每次推進時只需隨機抽取部分源粒子進行碰撞處理。隨機抽取時,對每一個源粒子產(chǎn)生一個隨機數(shù)R1(R1∈[0,1]),如R1

表1 粒子碰撞類型的選擇Tab.1 Selection of particle collision types

其中:υj為源粒子與靶粒子可能發(fā)生的第j種碰撞類型的概率;υ′為總碰撞概率。程序能處理多種碰撞過程,包括電子與中性原子的彈性碰撞、電子與中性原子的激發(fā)碰撞、電子與中性原子的電離碰撞、離子與中性原子的彈性碰撞及離子與中性原子的電荷交換碰撞等過程。在各種碰撞類型中,需處理碰撞后各粒子的位置、速度和方向。以電子與中性原子的彈性碰撞為例,碰撞后電子的運動方向是隨機的,散射角χ和子午面角φ為可表示為

(3)

其中:R3,R4為隨機數(shù)(R3,R4∈[0,1]);E為入射電子能量。

碰撞后電子的速度可表示為

(4)

其中:me,mi分別為電子和離子的質(zhì)量;v為入射電子速度。

1.3 基于并行遺傳算法的全電磁粒子模擬算法

遺傳算法是一種全局優(yōu)化算法,將全電磁粒子模擬算法與遺傳算法相結(jié)合,能實現(xiàn)真空電子器件的全局優(yōu)化設計。在采用遺傳算法對相對論返波管(relativistic backward wave oscillator,RBWO)進行優(yōu)化時,需將全電磁粒子模擬程序模塊計算得到的RBWO輸出功率作為每一個個體的適應度。由于計算適應度值要模擬改變結(jié)構(gòu)參數(shù)后的RBWO模型,直到器件有穩(wěn)定的功率輸出,并通過輸出功率的平均值得到該個體的適應度值,導致每個個體適應度值的計算時間較長。研制的真空電子器件優(yōu)化軟件采用并行遺傳優(yōu)化算法[16-17, 20],該方法能在國產(chǎn)超算上同時計算每個個體的適應度值,節(jié)省優(yōu)化設計所需的時間。

圖2為RBWO并行優(yōu)化遺傳算法流程框圖。

圖2 RBWO并行優(yōu)化遺傳算法流程框圖Fig.2 Flow chart of parallel genetic algorithm for RBWO

在遺傳算法中需要生成大量的偽隨機數(shù),軟件采用單獨的一個進程(主進程)負責偽隨機數(shù)的生成,采用其他進程(子進程)完成適應度值的計算。在第一步時,主進程讀取優(yōu)化參數(shù)設置的取值范圍,并通過偽隨機數(shù)生成算法為每個優(yōu)化的參數(shù)在取值范圍內(nèi)產(chǎn)生一組新的參數(shù)值,同時對參數(shù)進行浮點數(shù)編碼;主進程將新生成的參數(shù)分組,每一組均包含了所有需要優(yōu)化的參數(shù),將每一組參數(shù)發(fā)送給相應的子進程;所有子進程均執(zhí)行第二步,每一個子進程通過接收的參數(shù)生成新的RBWO模型文件,調(diào)用全電磁粒子模擬計算模塊讀入模型文件進行數(shù)值計算,并將器件的輸出功率求平均值,得到該個體的適應度值,最后將計算的適應度值發(fā)送給主進程;主進程單獨執(zhí)行第三、四、五步,主進程在獲取每個個體的適應度值后,首先將適應度值最大的個體挑選出,并對剩下的個體進行選擇操作,如果某一個體的適應度值較大,則它被選擇存活到下一代的概率就變大,如圖2的第三步所示;主進程以一定概率將參數(shù)進行交叉以及變異操作產(chǎn)生新的參數(shù),如圖2的第四步、五步所示;最后主進程將新生成的優(yōu)化參數(shù)值再次分組,并將生成的參數(shù)值再次發(fā)送給子進程;每個子進程均需要執(zhí)行第六步,每個子進程根據(jù)接收到的參數(shù)值,生成新的RBWO模型文件,子進程調(diào)用全電磁粒子模擬計算模塊讀入模型文件進行模擬計算,獲得該個體的適應度值,并再次發(fā)送給主進程。在第七步,主進程在接收到每個子進程發(fā)送的適應度值后,判斷最佳優(yōu)化結(jié)果是否滿足指標要求,如果不滿足最優(yōu)指標則跳往第三步開始下一輪的參數(shù)優(yōu)化,直到完成目標優(yōu)化并輸出優(yōu)化結(jié)果。

1.4 2.5維全電磁粒子模擬軟件優(yōu)化RBWO

RBWO是一種能量轉(zhuǎn)換效率非常高的真空電子器件,在設計過程中需調(diào)節(jié)的參數(shù)非常多[20],可采用本文建立的優(yōu)化方法對RBWO的慢波結(jié)構(gòu)進行自動優(yōu)化。慢波結(jié)構(gòu)底部沿著r方向高度的取值范圍為0.011～0.015 m,頂部的取值范圍為0.015～0.02 m。對RBWO慢波結(jié)構(gòu)的9個底部位置和9個頂部位置所對應的高度進行浮點數(shù)編碼,模型的其他參數(shù)和驅(qū)動電壓保持不變,初始種群數(shù)設置為80,進化50代后,就可得到RBWO慢波的優(yōu)化結(jié)構(gòu)。采用UNIPIC軟件模擬優(yōu)化前均勻慢波結(jié)構(gòu)及優(yōu)化后非均勻慢波結(jié)構(gòu)返波管模型,優(yōu)化前后,返波管的平均輸出功率如圖3所示。均勻慢波結(jié)構(gòu)RBWO的模擬結(jié)果如圖3(a)所示,結(jié)果表明在12 ns后輸出的平均功率為760兆瓦;優(yōu)化慢波結(jié)構(gòu)參數(shù)后,在12 ns后輸出的平均功率為1.16 GW,如圖3(b)所示,輸出功率相比優(yōu)化前能夠提高52%。數(shù)值模擬結(jié)果表明,2.5維粒子模擬軟件UNIIPIC能夠準確模擬軸對稱真空電子器件,同時融合了并行遺傳算法的全電磁粒子模擬算法,能用于真空電子器件的優(yōu)化設計。

(a) Original

(b) Opitimized圖3 優(yōu)化前后,RBWO的平均輸出功率Fig.3 Average output power of original and optimized RBWO

2 3維全電磁粒子數(shù)值模擬技術(shù)

3維并行粒子模擬軟件UNIPIC-3D軟件可用于模擬3維真空電子器件內(nèi)電磁波與相對論電子束之間的非線性相互作用[8-9]。軟件集成的計算機輔助設計工具能用于構(gòu)建復雜的3維幾何模型,實現(xiàn)的結(jié)構(gòu)網(wǎng)格剖分模塊沿著x,y,z方向?qū)碗s3維模型進行共形結(jié)構(gòu)網(wǎng)格剖分,采用區(qū)域分解算法實現(xiàn)了核心計算模塊的并行計算。軟件能設置理想導體、有耗導體、端口模式加載及完全匹配層等電磁場邊界條件,帶電粒子可通過場致發(fā)射、束發(fā)射及爆炸發(fā)射等方式產(chǎn)生帶電粒子。

2.1 3維全電磁粒子模擬共形算法

隨著真空電子器件工作頻率的提高,特別是器件的工作頻率達到太赫茲波段,真空電子器件的高頻結(jié)構(gòu)也會變得更精細,通過縮小網(wǎng)格的尺寸能提高計算的精度,但會導致器件模擬計算量的大幅增加。因此,為更準確地描述真空電子器件中的細小結(jié)構(gòu),提高軟件模擬復雜器件結(jié)構(gòu)的計算精度,采用共形磁場計算技術(shù)提高復雜結(jié)構(gòu)器件的模擬計算能力。研制的3維粒子模擬軟件采用了基于ECT-CFDTD(enlarged cell technique conformal FDTD)原理的電磁場共形計算方法,通過基于計算機輔助設計方法的射線追蹤方法,獲得網(wǎng)格線與幾何體的交點,這些交點、網(wǎng)格頂點及幾何體之間的拓撲關(guān)系共同構(gòu)成了共形網(wǎng)格系統(tǒng)。在共形網(wǎng)格系統(tǒng)中定義了幾何量及物理量之間的拓撲關(guān)系,采用ECT-CFDTD方法將所有離散的電磁場方程轉(zhuǎn)換為局部有限積分方程,避免了求解大型矩陣,并保持不同區(qū)域的推進方程具有相同的表達方式。

ECT-CFDTD方法的基本原理如圖4所示。

(a) Conformal grid used in CFDTD method

(b) Conformal electron emitter圖4 ECT-CFDTD方法的基本原理Fig.4 Basic principles of the ECT-CFDTD

當對面積為Sp的小面元計算電磁場時,電場的推進方程保持與經(jīng)典的FDTD方法相同,只對磁場的推進方程進行修改,該方法能在有效保證器件模擬計算精度的同時,無須大幅減少時間步長。當網(wǎng)格中真空所占面積較小時磁場的計算方法可表示為

(5)

(6)

(7)

(8)

(9)

(10)

2.2 3維全電磁粒子模擬并行技術(shù)

采用3維全電磁PIC粒子模擬軟件模擬真空電子器件需耗費的計算資源非常多,采用并行計算的方法能同時使用多個進程對同一個器件模型進行模擬,大幅減少器件模擬所需的時間。研制的3維并行粒子模擬軟件采用規(guī)則分區(qū)的方式實現(xiàn)了粒子模擬軟件中電磁場計算的并行化。首先,在不同的方向上定義不同的分割位置;然后,根據(jù)分割位置,將全局網(wǎng)格區(qū)域分割為多個子區(qū)域,每個子區(qū)域?qū)粋€物理區(qū)域,同時也對應唯一的進程標號;同時,在物理區(qū)域定義的基礎上,定義了與物理區(qū)域?qū)臄U展區(qū)域,用于給出當前進程的局部網(wǎng)格區(qū)域。

以2維情況下的區(qū)域分解算法為例,圖5為并行計算的區(qū)域分解及擴展網(wǎng)格層。

(a) Computational domain decompostion

(b) Extented grid layer of physical domain圖5 并行計算的區(qū)域分解及擴展網(wǎng)格層Fig.5 Domain decomposition and extended grid layer by parallel computation

由圖5(a)可見,由Decomposition1和Decomposition2定義了2個分割位置,將全局網(wǎng)格分割為4個區(qū)域,同時也將全局物理區(qū)域分割為4個子物理區(qū)域,分別標示為PhysRgn1,PhysRgn2,PhysRgn3,PhysRgn4,分別對應進程1、進程2、進程3和進程4。將物理區(qū)域的上下邊界各擴展一層為物理區(qū)域添加守護網(wǎng)格層,圖5(b)所示為進程1所計算的PhysRgn1區(qū)域擴展得到的局部網(wǎng)格區(qū)域ExtendedRgn1。每一個進程的擴展網(wǎng)格區(qū)域與其他進程的物理區(qū)域進行求交運算,若當前進程的擴展網(wǎng)格區(qū)域與某個進程的物理區(qū)域求交運算結(jié)果不為零,則需進行電磁場數(shù)據(jù)的交換。

在上述規(guī)則分區(qū)的基礎上,實現(xiàn)了帶電粒子的并行推進。在每一個時間步長內(nèi),每個進程推進所計算物理區(qū)域內(nèi)部所有帶電粒子的運動方程,確定帶電粒子經(jīng)單個時間步后的速度和位置;由帶電粒子推進后所處的位置,判斷是否將該帶電粒子發(fā)送至其他相鄰進程,當帶電粒子位于當前進程的守護網(wǎng)格區(qū)域時,就需將粒子信息發(fā)送給相鄰網(wǎng)格區(qū)域的進程。

2.3 3維全電磁粒子軟件模擬真空電子器件

采用自行研制的3維共形全電磁粒子模擬軟件對雙電子注平板回旋管進行了數(shù)值模擬研究,圖6為雙電子束平板回旋管示意圖[21]。平板回旋管的工作原理如圖6(a)所示。圖6(b)中,黃色標示為帶狀回旋電子束,紅色為平板回旋管的金屬結(jié)構(gòu)。兩個帶狀電子回旋電子束回旋中心在y方向的位置分別位于y=b0/6和y=5b0/6,每一個帶狀回旋電子束的線電流密度均為18 A·cm-1。網(wǎng)格數(shù)為240×60×380,在國產(chǎn)超級計算機上采用228個進程進行并行計算,得到電場強度Ex的3維空間分布、平板回旋管輸出功率隨時間的變化關(guān)系及電子的實空間和相空間分布,分別如圖7-圖9所示。

(a) Projection of electron orbit in the planar gyrotron

(b) Planar gyrotron with double-beam圖6 雙電子束平板回旋管示意圖Fig.6 Schematic diagram of double-beam planar gyrotron

圖7 電場強度Ex的3維空間分布Fig.7 3D distribution of Ex

(a) Output power in x direction

(b) Output power in -x direction圖8 平板回旋管輸出功率隨時間的變化關(guān)系Fig.8 Output powers of the planar gyrotron vs. time

(a) In real space

(b) In phase space

由圖7可見,Ex在z方向的相位變化呈8π正弦穩(wěn)態(tài)分布,y方向的相位變化為3π正弦穩(wěn)態(tài)分布,電磁場的分布模式穩(wěn)定,沒有出現(xiàn)模式競爭問題。

由圖8可見,2個端口處輸出的功率大小相同,因功率輸出的方向相反,依據(jù)功率的計算公式所得到的幅值相反,輸出的總功率為24.2 kW。由圖9可見,隨著帶狀回旋電子束沿著z軸方向傳輸,經(jīng)與電磁波的相互作用,電子束在角向發(fā)生速度調(diào)制,將電子束的動能轉(zhuǎn)化為電磁波能量。

3 基于時域間斷有限元的全電磁粒子模擬技術(shù)

研制了基于時域間斷有限元空間的2.5維全電磁粒子模擬程序,通過引入輔助矯正勢抑制了電荷不守恒所帶來的偏差[22-23]。程序采用三角形網(wǎng)格對計算區(qū)域進行非結(jié)構(gòu)網(wǎng)格剖分,實現(xiàn)了基于標量基函數(shù)及矢量基函數(shù)的電磁場時域間斷有限元算法,采用總場/散射場邊界條件實現(xiàn)了電壓波加載邊界,采用復頻率完全匹配層實現(xiàn)了端口的截斷邊界條件,同時實現(xiàn)了軸對稱和金屬壁等電磁場邊界條件。針對非結(jié)構(gòu)網(wǎng)格中帶電粒子與電磁場之間的相互作用,實現(xiàn)了帶電粒子的快速定位及電荷和電流分配算法。

3.1 2維電磁場時域間斷有限元算法

(11)

(13)

再次應用分部積分,得到

(14)

(15)

其中,算子〈·〉,[·]的定義為

上標“+”表示在鄰接單元表面上電磁場值。

為了減少電荷分配不守恒帶來的偏差,在程序中采用了最近發(fā)展的修正勢Maxwell方程[24],表示為

(16)

其中:Φ,Ψ為修正勢,它們及相關(guān)輔助方程的引入是為減小數(shù)值計算中電荷不守恒帶來的偏差;κ和γ為修正系數(shù)。由式(16)可知,當Φ→0,Ψ→0時,方程退化為Maxwell方程。

從式(16)可推導得到

(17)

式(17)表明修正勢Φ和Ψ可用來描述電荷不守恒和磁荷不守恒的偏差,且這種偏差以波的形式傳遞。當參數(shù)κ>1,γ>1時,由式(17)可知,偏差的傳播速度遠高于電磁波的速度。因此,在求解帶有修正勢的Maxwell方程時,電荷不守恒及磁荷不守恒引起的偏差會以非?？斓乃俣葌鞑サ接嬎銋^(qū)域之外。含修正勢的Maxwell方程間斷伽遼金形式的離散與不含修正勢Maxwell方程間斷伽遼金形式的離散具有相同的步驟,各場量的空間近似仍采用Lagrange插值,在單元邊界上數(shù)值通量的計算會因引入修正勢而有所變化。

3.2 非結(jié)構(gòu)網(wǎng)格中粒子快速定位以及電流計算

全電磁粒子模擬算法在計算帶電粒子受到的洛倫茲力時,將帶電粒子看成半徑無限小質(zhì)點,需計算粒子所在位置處的電磁場值,并代入到洛倫茲力公式中,得到電磁場對帶電粒子的作用力。該電磁場值通過粒子所在單元及鄰近單元上的電磁場的空間插值得到。為得到粒子所在位置處的電磁場,需知道粒子所在物理單元的編號及粒子在該單元內(nèi)的局部坐標。為根據(jù)粒子的物理坐標快速查找粒子所在的單元編號,程序?qū)τ嬎阌騼?nèi)的所有非結(jié)構(gòu)網(wǎng)格單元采用了空間分群策略。將整個計算區(qū)域均勻劃分為由Nx×Ny個直角網(wǎng)格構(gòu)成的虛擬結(jié)構(gòu)網(wǎng)格,并將與同一個虛擬網(wǎng)格有空間重疊的非結(jié)構(gòu)三角單元組合為一個群。由于虛擬網(wǎng)格排列較為規(guī)則,通過粒子的物理坐標,可快速計算粒子所在虛擬結(jié)構(gòu)網(wǎng)格的索引,進而獲得粒子所在單元三角形網(wǎng)格群,提高在非結(jié)構(gòu)網(wǎng)格中定位粒子的效率。確定粒子所在的三角網(wǎng)格單元后,根據(jù)等參單元與物理單元的坐標變換公式計算得到粒子在等參單元上的局部坐標,表示為

(19)

(20)

獲得粒子所在等參單元的局部坐標后,其所在位置的電磁場可表示為

(21)

帶電粒子運動產(chǎn)生的電流作為電磁場方程的源項,假設宏粒子為具有一定大小的球形粒子,每個宏粒子的空間分布函數(shù)表示為

(22)

其中:α為空間分布多項式函數(shù)的指數(shù)項,R為設置的帶電粒子影響半徑,r為需計算電流分布所在位置的坐標 ;r0為粒子中心位置坐標。S(r-r0)在|r-r0|∈[0,R]范圍內(nèi)的積分為1,宏粒子運動產(chǎn)生的電流密度可表示為

I(r-r0)=evS(r-r0)

(23)

其中,e為宏粒子所帶電荷量。宏粒子電荷密度的空間分布可表示為

ρ(r-r0)=eS(r-r0)

(24)

3.3 真空二極管間斷有限元全電磁數(shù)值模擬研究

采用自行研制的2.5維時域間斷有限元全電磁粒子模擬軟件對同軸真空二極管進行了數(shù)值模擬,圖10為真空二極管計算模型的剖面圖。其中,二極管端注入電壓為1 kV,上升沿為0.5 ns,沿著軸向的引導磁場強度為1.0 T,二極管同軸注入端的內(nèi)外半徑分別為5 cm和5.9 cm,同軸二極管陰極環(huán)面的內(nèi)外半徑分別為2.4 cm和2.7 cm,陰陽極板之間的距離為0.5 cm。對二極管進行了2維非結(jié)構(gòu)網(wǎng)格的剖分,網(wǎng)格平均尺寸為1 mm。

(b) Triangular grid of vacuum diode圖10 真空二極管計算模型Fig.10 Simulation model of vacuum diode

對二極管陰極面的發(fā)射電流及陰極發(fā)射面r=2.65 cm位置處的法向電場強度進行了診斷,模擬結(jié)果如圖11所示。由圖11(a)可見,穩(wěn)定后的電流為0.892 A;由圖11(b)可見,穩(wěn)定后的電場強度為70.7 kV·m-1。 采用UNIPIC軟件對相同的物理模型進行了數(shù)值模擬得到,穩(wěn)定后的陰極電流為0.901 A,與間斷有限元全電磁粒子算法模擬結(jié)果的相對偏差為1%;計算得到穩(wěn)定后的電場強度為73.3 kV·m-1,與時域間斷有限元全電磁粒子算法模擬結(jié)果的相對偏差為3.6%。主要原因是兩套代碼在處理電流發(fā)射時均采用了唯象的Gauss電子爆炸發(fā)射模型,該模型需計算陰極面附近位置處的法向電場強度及陰極發(fā)射面位置處的電荷密度。在UNIPIC軟件中,法向電場值為距陰極面半網(wǎng)格長度位置處的電場值,同時電荷分配在正交網(wǎng)格的頂點處;在2.5維間斷有限元全電磁粒子模擬算法中,采用的是非結(jié)構(gòu)三角形網(wǎng)格。2種算法的計算結(jié)果仍能保持一致,表明時域間斷有限元全電磁粒子模擬算法能用于真空電子器件的數(shù)值研究。

(a) Beam current vs. time

(b) Electric field strength along z direction vs. time圖11 2.5維時域間斷有限元粒子模擬程序計算結(jié)果Fig.11 Simulation results by 2.5D particle-in-cell code basedon discontinuous Galerkin time-domain method

4 結(jié)論

20多年來,本文作者所在團隊在全電磁粒子模擬算法研究及軟件研制方面取得了系列進展:(1)研制了基于計算機輔助設計的2.5維全電磁粒子模擬軟件,開發(fā)了融合并行遺傳算法和全電磁粒子模擬算法真空電子器件的優(yōu)化軟件,具備對RBWO等真空電子器件的結(jié)構(gòu)參數(shù)及電學參數(shù)的優(yōu)化能力;(2)研制了3維并行共形粒子模擬軟件,具備復雜模型建模及數(shù)值模擬能力,能在國產(chǎn)超級計算機上對真空電子器件開展全尺寸3維結(jié)構(gòu)的數(shù)值模擬研究;(3)發(fā)展了時域間斷有限元全電磁粒子模擬技術(shù),采用非結(jié)構(gòu)三角網(wǎng)格實現(xiàn)了器件結(jié)構(gòu)的精確描述,大幅提高了全電磁粒子模擬模擬軟件計算真空電子器件的計算精度。

下一步,將會持續(xù)開展磁流體力學、高能射線碰撞電離、多物理場耦合及人工智能賦能等離子體數(shù)值模擬等先進數(shù)值模擬技術(shù)研究,進一步拓展粒子模擬軟件在系統(tǒng)電磁脈沖模擬、跨尺度等離子體模擬及高密度等離子體模擬等領(lǐng)域的應用研究。