黃 鰲,楊中海,黃 桃,胡 權(quán),金曉林

摘 要:利用粒子模擬方法(PIC),用Fortran語言進(jìn)行程序設(shè)計(jì),對(duì)無限大平行板間軸對(duì)稱電子束的運(yùn)動(dòng)進(jìn)行二維模擬。在該程序設(shè)計(jì)中,宏粒子采用環(huán)宏粒子模型,用體積加權(quán)模型將每一個(gè)宏粒子電量分配到格點(diǎn)上,相對(duì)于面積加權(quán)模型更加精確。與通常數(shù)值模擬方法相比,PIC方法得出各個(gè)宏粒子的瞬態(tài)運(yùn)動(dòng)特性。最后通過以Magic仿真軟件的計(jì)算結(jié)果為標(biāo)準(zhǔn)進(jìn)行驗(yàn)證,證明了該算法的準(zhǔn)確性。

關(guān)鍵詞:PIC;環(huán)宏粒子模型;體積加權(quán)模型;有限差分方法

中圖分類號(hào):TP311文獻(xiàn)標(biāo)識(shí)碼:B

文章編號(hào):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 引 言

近年來計(jì)算機(jī)技術(shù)的迅速發(fā)展,在傳統(tǒng)的解析理論分析和實(shí)驗(yàn)研究以外,發(fā)展了第3種強(qiáng)有力的重要研究方法:計(jì)算機(jī)模擬,即使用計(jì)算機(jī)模擬計(jì)算跟蹤大量微觀粒子的運(yùn)動(dòng),再對(duì)組成物體(包括氣態(tài)、液態(tài)、固態(tài)和等離子態(tài))的大量微觀粒子進(jìn)行統(tǒng)計(jì)平均,由此得到宏觀物體的物質(zhì)特征和運(yùn)動(dòng)規(guī)律。粒子模擬的特點(diǎn)在于直接從最基本的電磁運(yùn)動(dòng)規(guī)律和粒子運(yùn)動(dòng)的力學(xué)規(guī)律出發(fā),利用高速計(jì)算機(jī)直接求解完整的Maxwell方程組和Lorentz方程,追蹤每一個(gè)宏粒子的運(yùn)動(dòng),它能更真實(shí)的反映實(shí)際物理過程[1]。

這里將采用粒子模擬的方法對(duì)無限大平行板間軸對(duì)稱電子束進(jìn)行靜電PIC模擬程序設(shè)計(jì),算法的最終設(shè)計(jì)目標(biāo)為了實(shí)現(xiàn)任意形狀型號(hào)電子槍的通用軟件設(shè)計(jì)。

1 物理模型

本文所模擬的模型為無限大金屬平行板系統(tǒng)這一理想化模型。電子以某初速度從陰極板的某一圓形區(qū)域連續(xù)發(fā)射,到達(dá)陽極板后電子完全被吸收。

2 粒子模擬方法的二維靜電模型

無限大平行板間軸對(duì)稱電子束運(yùn)動(dòng)模擬采用二維靜電PIC模型。

根據(jù)以上的物理模型,將數(shù)值模擬分析分為以下4個(gè)部分:Poisson方程的求解得出空間電勢(shì)分布、空間電場(chǎng)的求解、粒子所在位置場(chǎng)的求解和粒子位置的更新。

2.1 求解Poisson方程

2φ=-ρ/ε0(1)

在旋轉(zhuǎn)對(duì)稱場(chǎng)的情況下,式(1)變?yōu)?

2φ祌2+1r?鄲摘祌+2φ祕(mì)2=-ρε0(2)

采用有限差分法(FDM),對(duì)式(2)進(jìn)行差分處理。采用如圖1所示的網(wǎng)格,軸線(z軸)處j=1,而0點(diǎn)落在第j行,網(wǎng)格為邊長(zhǎng)為h的正方形,通過對(duì)式(2)做等間距差分化后:

圖1 旋轉(zhuǎn)對(duì)稱場(chǎng)的等距網(wǎng)絡(luò)

2.1.1 對(duì)稱軸以外的點(diǎn)(j>1)

4φ0=φ2+φ4+[1+1/2(j-1)]φ1+

[1-1/2(j-1)]φ3+ρ/h2ε0(3)

2.1.2 對(duì)稱軸上(j=1)

4φ0=φ2+φ4+2φ1+ρ/h2ε0(4)

2.2 求解網(wǎng)格點(diǎn)上的電場(chǎng)

E0z=(φ2-φ4)/2h;E0r=(φ3-φ1)/2h(5)

2.3 求粒子所受到的力

采用如圖2所示面積加權(quán)法,粒子所受到的電場(chǎng)力為:

Fi(ri)=qiE=qi[E1A1+E2A2+E3A3+E4A4]/h2(6)

圖2 二維面積加權(quán)法

2.4 速度和位置的更新

粒子的洛侖茲運(yùn)動(dòng)方程:

dvi/dt=F(rj,t)/mi(7)

dri/dt=vi(t)(8)

用蛙跳算法來推動(dòng)粒子,將式(7),式(8)用時(shí)間中心差分得:

vin+1=vin+(Fin+1/2/mi)Δt(9)

rin+3/2=rin+1/2+vin+1Δt(10)

初始時(shí)刻t=0開始,由給定的粒子速度和位置分布v0i和r0i,求得r1/2i=r0i+v0i(Δt/2),進(jìn)而求得F1/2i。

2.5 模擬流程圖

模擬流程圖如圖3所示。

圖3 PIC模擬流程圖

3 計(jì)算模型

為節(jié)省計(jì)算機(jī)內(nèi)存和計(jì)算時(shí)間,在進(jìn)行計(jì)算機(jī)模擬時(shí)不必跟蹤每一個(gè)點(diǎn)粒子(微觀粒子),只需計(jì)算能代表這個(gè)荷電粒子群的一個(gè)粒子就可以,從而把這個(gè)荷電粒子群處理為一個(gè)特殊粒子,即所謂“宏粒子”。

兩板間的粒子束會(huì)產(chǎn)生空間電荷場(chǎng)(即自洽場(chǎng)),其大小由泊松方程求解,在由宏粒子的位置變化而得到電荷密度分布和求宏粒子所受到的電場(chǎng)力的過程中,需要借助面積加權(quán)模型[5] 來進(jìn)行電荷、電場(chǎng)分配。在本模型中,陰極板所發(fā)射出來的電子束為一軸對(duì)稱模型,可以采用二維的z-r坐標(biāo)系靜電模型來計(jì)算。在對(duì)電荷向各個(gè)網(wǎng)格點(diǎn)分配時(shí),面積加權(quán)模型在z-r坐標(biāo)系中的精度并不是很高,尤其在對(duì)稱軸附近的區(qū)域會(huì)產(chǎn)生較大誤差,在電子粒中心附近會(huì)產(chǎn)生振蕩,與實(shí)際不符[6]。程序設(shè)計(jì)采用一種體積加權(quán)模型,將宏粒子視為環(huán)形體,粒子在其內(nèi)均勻分布。隨著粒子的向前推進(jìn),環(huán)的半徑增大,但宏粒子環(huán)本身的z向?qū)挾群蛂向?qū)挾炔蛔?它的電量也不變。如圖2所示,環(huán)宏粒子的電量為qi,環(huán)橫截面的面積為A,則分配到1,2,3,4個(gè)格點(diǎn)所在環(huán)上的電量分別為qiA1/A,qiA2/A,qiA3/A,qiA4/A。進(jìn)而求得每一個(gè)格點(diǎn)上的電荷體密度。由泊松方程可得出每一格點(diǎn)的電勢(shì),從而得出格點(diǎn)上每一點(diǎn)的電場(chǎng)強(qiáng)度,采用面積加權(quán)的方法求出每個(gè)粒子上的電場(chǎng),進(jìn)而對(duì)粒子的速度和位置進(jìn)行更新,進(jìn)入下一個(gè)時(shí)間步長(zhǎng)。

4 模擬結(jié)果

這里所選用的物理模型是:相隔距離為s=0.01 m的兩無限大(只要兩板間的距離遠(yuǎn)小于板的線度就可視為無限大平行板)平行板,從陰極板的一半徑為r=0.001 m的圓形區(qū)域內(nèi)發(fā)射一束電子,電流密度為jcathode,粒子初速度為cell(i)%vz。程序設(shè)計(jì)的目標(biāo)是觀察粒子在不同時(shí)刻兩板間運(yùn)動(dòng)粒子位置狀態(tài)分布。

圖4和圖5是在如下初始條件分別通過Magic和該程序設(shè)計(jì)輸出的結(jié)果:

粒子初速度為:

cell(i)%vz= 8.381675d6 m/s

發(fā)射電流密度為:

jcathode=-2.5d3 A/m2

運(yùn)行總時(shí)間為:

t=4.5 d-9 s。

圖4 Magic輸出圖

圖5 本文程序設(shè)計(jì)輸出圖

通過輸出結(jié)果比較,其包絡(luò)面形狀基本一致。Magic結(jié)果輸出圖在陽極板的最大高度為:2.055 387 370 310 987 E-003 m,所設(shè)計(jì)程序結(jié)果在陽極板最大高度為:2.007 482 300 816 977 E-003 m。如果該程序?qū)⒕W(wǎng)格化分的更細(xì),兩者的差別將更小,證明本程序設(shè)計(jì)的正確性。

5 結(jié) 語

這里用粒子模擬的方法對(duì)無限大平行板間軸對(duì)稱粒子束運(yùn)動(dòng)進(jìn)行模擬,采用環(huán)宏粒子模型進(jìn)行程序設(shè)計(jì),在宏粒子電量向網(wǎng)格分配時(shí)采用比面積加權(quán)模型更加精確的體積加權(quán)模型(CIC)。相對(duì)于通常的計(jì)算機(jī)數(shù)值模擬方法,該程序設(shè)計(jì)通過追蹤每一個(gè)粒子,能觀察到任一時(shí)刻的粒子狀態(tài)(位置、速度)分布,能輸出每個(gè)宏粒子在任意時(shí)刻的位置和速度數(shù)據(jù),直接求解完整的Maxwell方程組和Lorentz方程,它能更真實(shí)地再現(xiàn)整個(gè)微觀過程。在相同的初始條件下,程序通過Magic通用軟件進(jìn)行驗(yàn)證是正確的。

參考文獻(xiàn)

[1]廖平,楊中海,雷文強(qiáng),等.微波管電子槍三維粒子模擬研究[J].強(qiáng)激光與粒子束,2004,16(3):353-357.

[2]王秉中.計(jì)算電磁學(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]李永東,何鋒,劉純亮.軸對(duì)稱平板二極管空間電荷限制流的二維效應(yīng)[J].強(qiáng)激光與粒子束,2005,17(6):913-916.

[6]盛新慶.計(jì)算電磁學(xué)要論[M].北京:科學(xué)出版社,2004.