任澤平， 王 玥， 陳再高

(西北核技術(shù)研究所， 西安 710024)

粒子模擬(particle in cell，PIC)方法基于第一性原理，通過求解帶電粒子運動的Newton-Lorentz力方程和電磁場演化的Maxwell方程組[1]，可以實現(xiàn)對真空電子器件中束-場互作用[2-5]及激光與等離子體作用[6]等復雜非線性過程的模擬。為了減小數(shù)值噪聲，PIC方法采用有限大小粒子模型，各個宏粒子直接與網(wǎng)格上的電磁場分量進行相互作用，忽略了粒子之間近程碰撞作用。PIC方法在帶電粒子密度不高時能很好適用，但對于低溫、高密度等離子體，帶電粒子之間碰撞頻率很高，需要在PIC方法中加入直接碰撞對粒子速度進行修正[7]。

Takizuka和Abe最早在粒子模擬程序中實現(xiàn)庫侖碰撞算法[8](簡稱TA算法)，用來近似動理學方程中Fokker-Planck碰撞項的作用。Miller提出了適用于不同權(quán)重粒子之間的碰撞算法[9]。Sentoku在TA算法的基礎(chǔ)上考慮了相對論修正[6]。Nanbu提出了使用累積散射角的碰撞算法(簡稱Nanbu算法)，降低了碰撞算法對時間步長的要求[10]。Perez將Nanbu碰撞算法推廣到相對論情形，并考慮了低溫修正[11]。

本文首先介紹了粒子模擬程序中TA碰撞算法和Nanbu碰撞算法的實現(xiàn)方法，然后比較了這兩種碰撞算法的精度和計算效率，最后采用庫侖碰撞算法對等離子體雙流體方程中動摩擦因數(shù)計算公式的合理性進行了驗證。

1 帶電粒子庫侖碰撞算法

PIC程序中通?？梢栽诹Ｗ铀俣韧七M之前或之后調(diào)用庫侖碰撞算法，加入粒子之間近程碰撞作用，對粒子的速度進行修正。本文假設(shè)所有宏粒子采用相同的權(quán)重，介紹TA算法和Nanbu算法這兩種最常用的庫侖碰撞算法的實現(xiàn)流程。

1.1 TA碰撞算法

在TA算法中，首先根據(jù)粒子所處的網(wǎng)格單元對粒子進行分組，然后將同一網(wǎng)格單元內(nèi)的粒子進行隨機配對，最后進行粒子對的散射，對粒子的速度進行修正，碰撞的散射角滿足一定分布。TA算法的具體步驟如下：

1) 假設(shè)模擬的系統(tǒng)中有α和β兩種帶電粒子，根據(jù)粒子所處的網(wǎng)格單元對粒子進行分組，將粒子的指針存儲到網(wǎng)格單元內(nèi)該類粒子對應(yīng)的指針數(shù)組中。

2) 將指針數(shù)組元素的順序隨機打亂。設(shè)某個數(shù)組中有N個元素，從最后一個元素開始，從前面的元素中隨機選出一個與它進行交換，以此類推，對第N-1,N-2,…，直到第2個元素執(zhí)行類似的操作。

3) 將網(wǎng)格單元內(nèi)的粒子進行配對。碰撞粒子對的選擇方法如圖 1所示。對于兩種粒子之間的碰撞，設(shè)網(wǎng)格單元內(nèi)兩種粒子的數(shù)量分別為Nα和Nβ，則總碰撞次數(shù)N=max(Nα,Nβ)，將存儲粒子指針的數(shù)組看成是循環(huán)的，分別從兩個數(shù)組的第一個元素開始，依次選出N個進行碰撞的粒子對。同種粒子之間的碰撞只需要在數(shù)組內(nèi)部依次選擇出粒子對。

(a) α-β (b) α-α 圖1 碰撞粒子對的選擇Fig.1 Selection of collision pairs

4) 在質(zhì)心系中看待兩個粒子之間的碰撞可以簡化對碰撞的處理[7]。設(shè)兩粒子在實驗室參考系中的質(zhì)量和速度分別為m1，m2，v1，v2，則質(zhì)心系的運動速度為

(1)

質(zhì)心系中兩粒子動量之和始終為0，質(zhì)量較大的粒子在質(zhì)心系中運動速度較慢。對于彈性碰撞，質(zhì)心系的速度在碰撞前后不發(fā)生改變。在質(zhì)心系中兩個粒子的速度分別為

(2)

(3)

式中，u為粒子1相對于粒子2的運動速度，即

u=v1-v2=v1,c-v2,c

(4)

碰撞使兩個粒子的相對速度沿原來的方向發(fā)生偏轉(zhuǎn)和旋轉(zhuǎn)，變?yōu)閡′，相對速度的大小不變。為了便于計算碰撞之后的相對速度，首先需要將u變換到以它的方向為z軸的坐標系(簡稱z系)中。可以驗證，式(5)中的矩陣滿足這樣的變換。這個變換不用顯式進行，只需要計算出cosθ，sinθ，cosφ和sinφ的值即可。

(5)

其中，角度θ和φ的定義如圖 2所示，cosθ=uz/u；sinθ=u⊥/u； cosφ=ux/u⊥； sinφ=uy/u⊥；u⊥=(ux2+uy2)1/2。

圖2 相對速度矢量方位角 Fig.2 Angle of relative velocity vector

相對速度矢量的偏轉(zhuǎn)和旋轉(zhuǎn)如圖 3所示。在z系中碰撞后的散射角和旋轉(zhuǎn)角分別為Θ和Φ，兩個粒子碰撞之后的相對速度變?yōu)閡′=(usinΘcosΦ,usinΘsinΦ,ucosΘ)。

圖3 相對速度矢量的偏轉(zhuǎn)和旋轉(zhuǎn) Fig.3 Scattering of relative velocity

通過式(6)中的反變換，將z系中碰撞后的相對速度再變換回質(zhì)心坐標系。

(6)

計算出碰撞之后的相對速度u′，就能得到碰撞后兩個粒子在質(zhì)心系中的速度：

(7)

(8)

進一步可以得到碰撞后兩粒子在實驗室參考系中的速度：

(9)

(10)

5)碰撞角度的計算。旋轉(zhuǎn)角Φ可在[0, 2π]區(qū)間隨機選取。散射角Θ可根據(jù)庫侖碰撞散射角所滿足的分布函數(shù)抽樣得到[8]，令δ=tan(Θ/2)，則δ滿足高斯分布，且該分布的方差為

δ2

(11)

其中，m12=(m1m2)/(m1+m2)；q1和q2分別為兩個物理粒子的電荷量；ε0為真空介電常數(shù)；lnΛ為庫侖對數(shù)，大小一般在10左右；nL=min(n1,n2)為兩種粒子數(shù)密度較小的一個，選擇nL是為了在碰撞過程中能以相同的時間步長推進所有粒子[10]。

由Box-Muller方法可以抽樣得到δ值[12]：

(12)

其中，R1和R2為 [0, 1]區(qū)間內(nèi)均勻分布的隨機數(shù)。進而可得散射角的正弦值和余弦值：

(13)

(14)

TA碰撞算法實現(xiàn)起來較為簡單，但要求模擬的時間步長遠小于平均碰撞時間間隔。當?shù)入x子體碰撞頻率很高時，采用TA碰撞算法的時間步長會非常短。

1.2 Nanbu碰撞算法

Nanbu碰撞算法使用帶電粒子發(fā)生多次小角度碰撞后累積的散射角計算速度的改變，計算效率較TA算法更高，允許使用較大的時間步長[10]。Nanbu碰撞算法中粒子分組、配對和碰撞過程的處理方法與TA算法中基本一致，但對散射角的計算方法不同。散射角的計算步驟如下：

1) 計算碰撞參數(shù)s：

(15)

2) 解式 (16)中的方程，計算碰撞參數(shù)A：

(16)

3) 累積散射角χ滿足分布函數(shù)：

(17)

根據(jù)累積散射角分布函數(shù)，抽樣得到散射角的值：

(18)

其中，U為[ 0, 1 ]區(qū)間均勻分布的隨機數(shù)。

可以根據(jù)文獻[11]給出的擬合公式，簡化步驟2)中對方程的求解：

1) 如果s<0.1，則

cosχ=1+slnU

(19)

為了避免計算得到的cosχ<-1，要限制隨機數(shù)U不能取太小。可令U=(1-1×10-12)×R+0.5×10-12，其中，R為[0，1]區(qū)間均勻分布的隨機數(shù)。

2) 如果0.1≤s<3，則

1/A=0.005 695 8+0.956 020 2s-0.508 139s2+

0.479 139 06s3-0.127 889 75s4+0.023 895 67s5

(20)

3) 如果3≤s<6，則

A=3exp(-s)

(21)

4) 如果s≥6，則

cosχ=2U-1

(22)

Nanbu碰撞算法采用累積的散射角計算碰撞，時間步長最大可取為與平均碰撞時間間隔一致，因而比TA碰撞算法效率更高[13]。

庫侖碰撞算法一般只對同一網(wǎng)格單元內(nèi)的粒子進行碰撞，碰撞不改變粒子的位置，只改變粒子的速度。對網(wǎng)格單元內(nèi)的粒子進行隨機排序并組成碰撞粒子對，主要是為了模擬各種可能的碰撞狀態(tài)，碰撞粒子對的選擇與粒子在網(wǎng)格單元內(nèi)的位置無關(guān)。模擬庫侖碰撞，理論上應(yīng)在等離子體德拜長度大小的網(wǎng)格內(nèi)選取碰撞粒子對，但只要粒子的溫度、密度梯度不太大，就可以降低對網(wǎng)格尺寸的要求[6]。

2 庫侖碰撞算法的應(yīng)用

2.1 兩種碰撞算法的比較

等離子體中的電子有時會出現(xiàn)沿x軸方向的溫度和垂直于x軸方向的溫度不相等，此時，對各個溫度分別定義為

(23)

(24)

(25)

電子通過碰撞可使其溫度分布趨于各向同性。由于能量守恒，T的值保持為常量。當|T//-T⊥|

(26)

(27)

(28)

(29)

其中，ΔT=T⊥-T∥。由于T為常量，可得兩個方向溫差隨時間的變化關(guān)系為

ΔT(t)=ΔT0exp(-t/τ)

(30)

分別采用兩種碰撞算法模擬等離子體中電子溫度平衡的過程，模擬只考慮電子之間的碰撞。初始參數(shù)設(shè)置為：電子密度ne=1026/m3，垂直x軸方向的溫度T⊥=80 eV，沿x軸方向的溫度T∥=100 eV。選取v0Δt=0.07和v0Δt=0.37兩種時間步長，其中，v0=1/τ0為特征碰撞頻率。模擬得到電子兩個方向上的溫差隨時間的變化關(guān)系如圖4所示，其中還給出了相應(yīng)的理論解以茲比較。從圖4可以看出，在相同的時間步長下，采用Nanbu碰撞算法得到的結(jié)果與理論解更加吻合。

圖4 電子兩個方向上的溫差隨時間的變化 Fig.4 Difference of electron temperatures along two directions

2.2 雙流體方程中動摩擦因數(shù)計算公式的驗證

流體模擬方法是等離子體數(shù)值模擬的一個重要分支，適用于研究等離子體宏觀輸運現(xiàn)象[7]。當?shù)入x子體中只包含兩種帶電粒子成分(電子e和離子i)時，可以采用雙流體模型進行描述，兩種流體通過庫侖碰撞發(fā)生耦合。假設(shè)兩種流體的各物理量在空間均勻分布，且只在一個方向有相對運動，當只考慮它們之間碰撞引起的動量和能量交換時，可以用式(31)--式(34)描述其速度和溫度的平衡過程[14]：

(31)

(32)

(33)

(34)

其中，ne，ue，Te分別為電子的密度、速度和溫度；mei=memi/(me+mi)，me和mi分別為電子和離子的質(zhì)量；νei和μei分別為動摩擦因數(shù)和熱平衡因數(shù)；式(31)中-νei(ue-ui)表示摩擦阻力項。式(32)中等號右側(cè)第一項表示能量耗散項(歐姆加熱)，第二項表示兩種粒子之間的熱交換作用，當兩種流體速度不一致時會存在耗散，動能轉(zhuǎn)化為內(nèi)能；當兩種流體之間存在溫差時又會有熱交換作用。式(33)和式(34)中帶上標“0”的量為初始時刻參數(shù)。

輸運系數(shù)是雙流體方程求解的關(guān)鍵，一般由動理學理論推導而來。Decoster給出的υei和μei的計算公式分別為[15]

(35)

(36)

為了檢驗式(35)和式(36)的準確性，設(shè)離子電荷數(shù)Zi=1，離子質(zhì)量為電子質(zhì)量的20倍，即mi=20me，初始時刻Te=Ti=5 eV，流體沿x軸正方向運動，電子初速度ue=107m·s-1，離子初速度ui=0，電子與離子密度相等，即ne=ni=1026m-3，求解式(31)—式(34)得到電子和離子運動速度和溫度隨時間的變化；同時，用PIC庫侖碰撞算法在二維x-y坐標系對該過程進行模擬，模擬結(jié)果用來檢驗流體方程解的正確性。兩種途徑所得結(jié)果的對比如圖5和圖6所示。

由圖5和圖6可以看出，采用Decoster所給輸運系數(shù)的流體方程得到的結(jié)果與采用PIC模擬得到的結(jié)果吻合較好；電子和離子的動量很快達到平衡，溫度則緩慢達到平衡；歐姆加熱主要作用于質(zhì)量較輕的電子，電子的溫度起初由于歐姆加熱而迅速升高，隨著電子和離子相對運動速度逐漸減小，歐姆加熱作用逐漸減弱，電子與離子之間通過熱交換緩慢達到溫度平衡。

同時，由圖5和圖6還可以看出，流體方程的解與PIC模擬結(jié)果的差異在于達到平衡的時間及溫度峰值的大小，這可能是因為流體方法是在假設(shè)電子滿足漂移麥克斯韋分布的基礎(chǔ)上進行的，而真實的情況是電子在減速過程中是偏離麥克斯韋分布的，平行于速度方向和垂直于速度方向上的溫度并不相等。圖 7是模擬得到的電子在這兩個方向上的溫度隨時間變化曲線，它們一開始并不相等。

圖5 電子和離子速度隨時間的變化 Fig.5 Velocity of electron and ion species

圖6 電子和粒子溫度隨時間的變化 Fig.6 Temperature of electron and ion species

圖7 電子兩個方向上的溫度隨時間的變化 Fig.7 Time histories of parallel and perpendicular temperatures of electron

3 結(jié)論

詳細介紹了在PIC模擬程序中實現(xiàn)庫侖碰撞算法的過程，從質(zhì)心系的角度介紹了碰撞對速度產(chǎn)生的改變，有利于理解庫侖碰撞算法的原理。給出了Nanbu碰撞算法中計算散射角的擬合公式，避免了求解超越方程的過程，可以提高Nanbu碰撞算法的計算效率。通過加入庫侖碰撞作用，可以將PIC方法推廣到對高密度等離子體的模擬。通過比較兩種碰撞算法的模擬結(jié)果，發(fā)現(xiàn)Nanbu碰撞算法在時間步長較大時依然能給出較為準確的模擬結(jié)果。最后，采用Nanbu碰撞算法驗證了Decoster給出的等離子體雙流體方程動摩擦因數(shù)計算公式的合理性。

[1]BIRDSALL C K, LANGDON A B. Plasma Physics via Computer Simulation[M]. New York: McGraw-Hill, 1985.

[2]BENFORD J, SWEGLE J, SCHAMILOGLU E. High Power Microwaves[M]. New York: Taylor and Francis, 2007.

[3]WANG J G, ZHANG D H, LIU C L, et al. UNIPIC code for simulations of high power microwave devices [J]. Phys Plasmas, 2009, 16(3): 033108.

[4]王建國. 真空電子器件的粒子模擬方法[J]. 現(xiàn)代應(yīng)用物理，2013, 4(3): 251-262. (WANG Jian-guo. Particle simulation method of vacuum electronic devices[J]. Modern Applied Physics, 2013, 4(3): 251-262. )

[5]WANG J G, CHEN Z G, WANG Y, et al. Three-dimensional parallel UNIPIC-3D code for simulations of high-power microwave devices [J]. Phys Plasmas, 2010, 17(7): 073107.

[6]SENTOKU Y, KEMP A J. Numerical method for particle simulations at extreme densities and temperatures[J]. J Comput Phys, 2008, 227(14): 6 846-6 861.

[7]鄭春開. 等離子體物理[M]. 北京: 北京大學出版社, 2009.(ZHENG Chun-kai. Physics of Plasmas[M]. Beijing: Peking University Press, 2009.)

[8]TAKIZUKA T, ABE H. A binary collision model for plasma simulation with a particle code[J]. J Comput Phys, 1997, 25(3): 205-219.

[9]MILLER R H, COMBI M R. A coulomb collision algorithm for weighted particles simulations [J]. Geophysical Reasarch Letters, 1994, 21(16): 1 735-1 738.

[10]NANBU K, YONEMURA S. Weighted particles in coulomb collision simulations based on the theory of a cumulative scattering angle[J]. J Comput Phys, 1998, 145(2): 639-654.

[11]PEREZ F,GREMMILLET L, DECOSTER A, et al. Improved modeling of relativistic collision and colliosional ionization in particle in cell code[J]. Phys Plasma, 2012, 19(8): 083104.

[12]WILLIAM T. Fortran Numerical Recipes [M]. 2nd ed. Cambridge University Press, 1992.

[13]WANG C, LIN T, CAFLISCH R, et al. Particle simulation of coulomb collision: Comparing the method of Takizuka & Abe and Nanbu[J]. J Comput Phys, 2008, 227(9): 4 308-4 329.

[14]RAMBO P, PROCASSINI R J. A comparation of kinetic and multifluid simulations of laser-produced colliding plasmas[J]. Phys Plasma, 1995, 2(8): 3 130-3 145.

[15]JONES M, LEMONS R. A grid-based coulomb collision model for PIC code[J]. J Comput Phys, 1996, 123: 169-181.