黃 鰲,楊中海,黃 桃,胡 權(quán),金曉林
摘 要:利用粒子模擬方法(PIC),用Fortran語言進行程序設(shè)計,對無限大平行板間軸對稱電子束的運動進行二維模擬。在該程序設(shè)計中,宏粒子采用環(huán)宏粒子模型,用體積加權(quán)模型將每一個宏粒子電量分配到格點上,相對于面積加權(quán)模型更加精確。與通常數(shù)值模擬方法相比,PIC方法得出各個宏粒子的瞬態(tài)運動特性。最后通過以Magic仿真軟件的計算結(jié)果為標(biāo)準(zhǔn)進行驗證,證明了該算法的準(zhǔn)確性。
關(guān)鍵詞:PIC;環(huán)宏粒子模型;體積加權(quán)模型;有限差分方法
中圖分類號:TP311文獻標(biāo)識碼:B
文章編號:1004-373X(2009)12-108-03
PIC Simulated Program Design of Two-dimensional Axial Symmetry
Electron Beam in Infinite Parallel
HUANG Ao,YANG Zhonghai,HUANG Tao,HU Quan,JIN Xiaolin
(School of Physical Electronics,University of Electronic Science and Technology of China,Chengdu,610054,China)
Abstract:The movement of axial symmetry electron beam located between infinite parallel by using particle simulation method in Fortran platform.Ring structure macro particle is distributed on grid by the volume-weighted model method and can be more accurate than face-weighted method.Compared with traditional numerical simulation method,the position and velocity of macro electron is calculated at any time step in this PIC method.According to the comparison with the calculation result from MAGIC,the accuracy of such algorithm mentioned in this paper has been verified.
Keywords:PIC;ring marco particle model;volume-weighted mode;finite difference method
0 引 言
近年來計算機技術(shù)的迅速發(fā)展,在傳統(tǒng)的解析理論分析和實驗研究以外,發(fā)展了第3種強有力的重要研究方法:計算機模擬,即使用計算機模擬計算跟蹤大量微觀粒子的運動,再對組成物體(包括氣態(tài)、液態(tài)、固態(tài)和等離子態(tài))的大量微觀粒子進行統(tǒng)計平均,由此得到宏觀物體的物質(zhì)特征和運動規(guī)律。粒子模擬的特點在于直接從最基本的電磁運動規(guī)律和粒子運動的力學(xué)規(guī)律出發(fā),利用高速計算機直接求解完整的Maxwell方程組和Lorentz方程,追蹤每一個宏粒子的運動,它能更真實的反映實際物理過程[1]。
這里將采用粒子模擬的方法對無限大平行板間軸對稱電子束進行靜電PIC模擬程序設(shè)計,算法的最終設(shè)計目標(biāo)為了實現(xiàn)任意形狀型號電子槍的通用軟件設(shè)計。
1 物理模型
本文所模擬的模型為無限大金屬平行板系統(tǒng)這一理想化模型。電子以某初速度從陰極板的某一圓形區(qū)域連續(xù)發(fā)射,到達陽極板后電子完全被吸收。
2 粒子模擬方法的二維靜電模型
無限大平行板間軸對稱電子束運動模擬采用二維靜電PIC模型。
根據(jù)以上的物理模型,將數(shù)值模擬分析分為以下4個部分:Poisson方程的求解得出空間電勢分布、空間電場的求解、粒子所在位置場的求解和粒子位置的更新。
2.1 求解Poisson方程
2φ=-ρ/ε0(1)
在旋轉(zhuǎn)對稱場的情況下,式(1)變?yōu)?
2φ祌2+1r?鄲摘祌+2φ祕2=-ρε0(2)
采用有限差分法(FDM),對式(2)進行差分處理。采用如圖1所示的網(wǎng)格,軸線(z軸)處j=1,而0點落在第j行,網(wǎng)格為邊長為h的正方形,通過對式(2)做等間距差分化后:
圖1 旋轉(zhuǎn)對稱場的等距網(wǎng)絡(luò)
2.1.1 對稱軸以外的點(j>1)
4φ0=φ2+φ4+[1+1/2(j-1)]φ1+
[1-1/2(j-1)]φ3+ρ/h2ε0(3)
2.1.2 對稱軸上(j=1)
4φ0=φ2+φ4+2φ1+ρ/h2ε0(4)
2.2 求解網(wǎng)格點上的電場
E0z=(φ2-φ4)/2h;E0r=(φ3-φ1)/2h(5)
2.3 求粒子所受到的力
采用如圖2所示面積加權(quán)法,粒子所受到的電場力為:
Fi(ri)=qiE=qi[E1A1+E2A2+E3A3+E4A4]/h2(6)
圖2 二維面積加權(quán)法
2.4 速度和位置的更新
粒子的洛侖茲運動方程:
dvi/dt=F(rj,t)/mi(7)
dri/dt=vi(t)(8)
用蛙跳算法來推動粒子,將式(7),式(8)用時間中心差分得:
vin+1=vin+(Fin+1/2/mi)Δt(9)
rin+3/2=rin+1/2+vin+1Δt(10)
初始時刻t=0開始,由給定的粒子速度和位置分布v0i和r0i,求得r1/2i=r0i+v0i(Δt/2),進而求得F1/2i。
2.5 模擬流程圖
模擬流程圖如圖3所示。
圖3 PIC模擬流程圖
3 計算模型
為節(jié)省計算機內(nèi)存和計算時間,在進行計算機模擬時不必跟蹤每一個點粒子(微觀粒子),只需計算能代表這個荷電粒子群的一個粒子就可以,從而把這個荷電粒子群處理為一個特殊粒子,即所謂“宏粒子”。
兩板間的粒子束會產(chǎn)生空間電荷場(即自洽場),其大小由泊松方程求解,在由宏粒子的位置變化而得到電荷密度分布和求宏粒子所受到的電場力的過程中,需要借助面積加權(quán)模型[5] 來進行電荷、電場分配。在本模型中,陰極板所發(fā)射出來的電子束為一軸對稱模型,可以采用二維的z-r坐標(biāo)系靜電模型來計算。在對電荷向各個網(wǎng)格點分配時,面積加權(quán)模型在z-r坐標(biāo)系中的精度并不是很高,尤其在對稱軸附近的區(qū)域會產(chǎn)生較大誤差,在電子粒中心附近會產(chǎn)生振蕩,與實際不符[6]。程序設(shè)計采用一種體積加權(quán)模型,將宏粒子視為環(huán)形體,粒子在其內(nèi)均勻分布。隨著粒子的向前推進,環(huán)的半徑增大,但宏粒子環(huán)本身的z向?qū)挾群蛂向?qū)挾炔蛔?它的電量也不變。如圖2所示,環(huán)宏粒子的電量為qi,環(huán)橫截面的面積為A,則分配到1,2,3,4個格點所在環(huán)上的電量分別為qiA1/A,qiA2/A,qiA3/A,qiA4/A。進而求得每一個格點上的電荷體密度。由泊松方程可得出每一格點的電勢,從而得出格點上每一點的電場強度,采用面積加權(quán)的方法求出每個粒子上的電場,進而對粒子的速度和位置進行更新,進入下一個時間步長。
4 模擬結(jié)果
這里所選用的物理模型是:相隔距離為s=0.01 m的兩無限大(只要兩板間的距離遠小于板的線度就可視為無限大平行板)平行板,從陰極板的一半徑為r=0.001 m的圓形區(qū)域內(nèi)發(fā)射一束電子,電流密度為jcathode,粒子初速度為cell(i)%vz。程序設(shè)計的目標(biāo)是觀察粒子在不同時刻兩板間運動粒子位置狀態(tài)分布。
圖4和圖5是在如下初始條件分別通過Magic和該程序設(shè)計輸出的結(jié)果:
粒子初速度為:
cell(i)%vz= 8.381675d6 m/s
發(fā)射電流密度為:
jcathode=-2.5d3 A/m2
運行總時間為:
t=4.5 d-9 s。
圖4 Magic輸出圖
圖5 本文程序設(shè)計輸出圖
通過輸出結(jié)果比較,其包絡(luò)面形狀基本一致。Magic結(jié)果輸出圖在陽極板的最大高度為:2.055 387 370 310 987 E-003 m,所設(shè)計程序結(jié)果在陽極板最大高度為:2.007 482 300 816 977 E-003 m。如果該程序?qū)⒕W(wǎng)格化分的更細,兩者的差別將更小,證明本程序設(shè)計的正確性。
5 結(jié) 語
這里用粒子模擬的方法對無限大平行板間軸對稱粒子束運動進行模擬,采用環(huán)宏粒子模型進行程序設(shè)計,在宏粒子電量向網(wǎng)格分配時采用比面積加權(quán)模型更加精確的體積加權(quán)模型(CIC)。相對于通常的計算機數(shù)值模擬方法,該程序設(shè)計通過追蹤每一個粒子,能觀察到任一時刻的粒子狀態(tài)(位置、速度)分布,能輸出每個宏粒子在任意時刻的位置和速度數(shù)據(jù),直接求解完整的Maxwell方程組和Lorentz方程,它能更真實地再現(xiàn)整個微觀過程。在相同的初始條件下,程序通過Magic通用軟件進行驗證是正確的。
參考文獻
[1]廖平,楊中海,雷文強,等.微波管電子槍三維粒子模擬研究[J].強激光與粒子束,2004,16(3):353-357.
[2]王秉中.計算電磁學(xué)[M].北京:科學(xué)出版社,2002.
[3]Freeman J C.Preliminary Study of Electron Emission for Use in the PIC Portion of MAFIA[R].NASA/TM220012210890,Glenn Research Center,2001.
[4]Birdsall C K,Fuss D.Clouds-in-Clouds,Clouds-in-Cell′s Physics for Many Body Plasma Simulation[J].Comput.Phys.,1969(3):494-511.
[5]李永東,何鋒,劉純亮.軸對稱平板二極管空間電荷限制流的二維效應(yīng)[J].強激光與粒子束,2005,17(6):913-916.
[6]盛新慶.計算電磁學(xué)要論[M].北京:科學(xué)出版社,2004.