粒子模擬中庫侖碰撞算法的實現(xiàn)

任澤平， 王 玥， 陳再高

(西北核技術研究所， 西安 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提出了適用于不同權重粒子之間的碰撞算法[9]。Sentoku在TA算法的基礎上考慮了相對論修正[6]。Nanbu提出了使用累積散射角的碰撞算法(簡稱Nanbu算法)，降低了碰撞算法對時間步長的要求[10]。Perez將Nanbu碰撞算法推廣到相對論情形，并考慮了低溫修正[11]。

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

1 帶電粒子庫侖碰撞算法

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

1.1 TA碰撞算法

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

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

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

3) 將網(wǎng)格單元內(nèi)的粒子進行配對。碰撞粒子對的選擇方法如圖 1所示。對于兩種粒子之間的碰撞，設網(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) 在質心系中看待兩個粒子之間的碰撞可以簡化對碰撞的處理[7]。設兩粒子在實驗室參考系中的質量和速度分別為m1，m2，v1，v2，則質心系的運動速度為

(1)

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

(2)

(3)

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

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

(4)

碰撞使兩個粒子的相對速度沿原來的方向發(fā)生偏轉和旋轉，變?yōu)閡′，相對速度的大小不變。為了便于計算碰撞之后的相對速度，首先需要將u變換到以它的方向為z軸的坐標系(簡稱z系)中?？梢则炞C，式(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

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

圖3 相對速度矢量的偏轉和旋轉 Fig.3 Scattering of relative velocity

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

(6)

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

(7)

(8)

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

(9)

(10)

5)碰撞角度的計算。旋轉角Φ可在[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不能取太小?？闪頤=(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)的位置無關。模擬庫侖碰撞，理論上應在等離子體德拜長度大小的網(wǎng)格內(nèi)選取碰撞粒子對，但只要粒子的溫度、密度梯度不太大，就可以降低對網(wǎng)格尺寸的要求[6]。

2 庫侖碰撞算法的應用

2.1 兩種碰撞算法的比較

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

(23)

(24)

(25)

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

(26)

(27)

(28)

(29)

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

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

(30)

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

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

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

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

(31)

(32)

(33)

(34)

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

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

(35)

(36)

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

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

同時，由圖5和圖6還可以看出，流體方程的解與PIC模擬結果的差異在于達到平衡的時間及溫度峰值的大小，這可能是因為流體方法是在假設電子滿足漂移麥克斯韋分布的基礎上進行的，而真實的情況是電子在減速過程中是偏離麥克斯韋分布的，平行于速度方向和垂直于速度方向上的溫度并不相等。圖 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 結論

詳細介紹了在PIC模擬程序中實現(xiàn)庫侖碰撞算法的過程，從質心系的角度介紹了碰撞對速度產(chǎn)生的改變，有利于理解庫侖碰撞算法的原理。給出了Nanbu碰撞算法中計算散射角的擬合公式，避免了求解超越方程的過程，可以提高Nanbu碰撞算法的計算效率。通過加入庫侖碰撞作用，可以將PIC方法推廣到對高密度等離子體的模擬。通過比較兩種碰撞算法的模擬結果，發(fā)現(xiàn)Nanbu碰撞算法在時間步長較大時依然能給出較為準確的模擬結果。最后，采用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)代應用物理，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.