基于Matlab的二維PIC/MCC模型的實現

王俊杰 鄭錦華,2 魏新煦 吳 雙 許 璐

（1. 鄭州大學化工與能源學院熱能系統節能技術與裝備教育部工程研究中心，鄭州 450001；2. 鄭州瑞邦石油機械有限公司，鄭州 450001）

基于Matlab的二維PIC/MCC模型的實現

王俊杰1鄭錦華1,2魏新煦1吳 雙1許 璐1

（1. 鄭州大學化工與能源學院熱能系統節能技術與裝備教育部工程研究中心，鄭州 450001；2. 鄭州瑞邦石油機械有限公司，鄭州 450001）

為了實現限定條件下氣體放電的二維模擬，基于Matlab通過編程手段開發模擬程序,并用此模擬程序對直流低壓條件下氬氣的放電過程進行驗證。對電子速度進行麥克斯韋初始化，并進行結果驗證。提取50個時間步長內的電子數變化趨勢，驗證了氬氣放電過程中電子崩的發生；并且通過100個時間步長結束時的電子位置分布，可以觀察到陽極鞘層的產生；根據解泊松方程至200步和300步結果的對比研究，驗證在本文的模擬條件下200步解泊松方程是足夠的；對電子速度均方根處理，得到溫度分布示意圖，驗證了直流氬氣放電時陽極板附近溫度高的實驗現象。至此，基于Matlab開發的二維粒子/蒙特卡洛碰撞的耦合（PIC/MCC）程序，能夠進行氣體放電模擬。

Matlab；PIC/MCC；直流放電；數學模型；麥克斯韋分布

高電壓絕緣氣體放電模擬一直以來是電工學領域的熱點，PIC/MCC模型是低壓條件下常用的放電模型，該模型自20世紀50年代被南部等人提出，至今幾十年被不斷完善[1-2]。很多國內的研究者采用這種方法給電器技術領域帶來了發展，如：湖南大學的汪沨教授采用PIC模型對10%～90%SF6/N2混合氣體中絕緣子沿面放電特性進行研究[3]；上海交通大學的研究團隊采用改進的MCC方法研究SF6和CO2混合氣體電子崩的參數[4]；沈陽工業大學研究團隊采用PIC方法對斷路器短間隙氣體放電過程進行模擬研究[5]。

氣體放電數學模型有：流體模型即N-S模型也稱宏觀模型；粒子（PIC）模型也稱微觀模型；二者混合模型三種[6]。本文針對低壓氣體放電進行粒子（PIC）模擬，因為在低壓條件下（＜10mtorr）流體模型并不適用[7]。本文首次采用Matlab編程的方式，實現二維PIC/MCC模型。并對氬氣直流條件下的放電過程進行模擬驗證。氬氣是微電子器件表面處理的常用氣體，并在各種放電條件下被研究。因此選擇氬氣作為驗證氣體。

Matlab具有直接面向矩陣的特點，決定其在粒子模型實現方面有很大的優勢。Matlab程序簡潔且易于理解，并且可以與其他匯編語言進行混合編程，解決編程問題手段多變[8]，故采用Matlab來實現PIC/MCC模型的程序開發。

二維放電模型與一維模擬[9]相比有以下不同：①二維模型需要考慮平板邊界條件對放電過程的影響，而一維模型無需考慮；②二維模擬能夠直觀的反映出在放電過程中，電子與離子的運動情況和各個參數的變化；③對于二維模型而言，在相同條件下的粒子放電模擬，要達到同樣的參數精度，網格數目是一維時的二次方倍，計算時間成倍的增加，工作效率降低[10]。

本文從離散玻爾茲曼方程、建立完善的二維數學模型，到程序流程建立以及編程實現，再到氬氣放電驗證等角度，來討論直流氣體放電過程中的粒子模擬，并將模擬結果與實驗現象以及其他研究者所取得的相關模擬結果進行比較。

1 PIC/MCC模型的Matlab實現

對直流低壓條件下氬氣PIC/MCC放電模型的研究基于以下幾點假設：①忽略電子運動形成的磁場；②忽略帶電粒子間的庫倫作用，即忽略庫倫碰撞；③邊界為良導體，電子到達表面會被吸收；④為了達到驗證目的并減少工作量，采用矩形放電區域，網格為正方形網格。

圖1為編程實現模型的計算流程圖，先要給各個粒子賦初始值。（其中包括粒子所帶電荷，粒子的位置，初始速度等）進入第一個循環，按照面積權重將粒子所帶電荷分配到各自鄰近的網格節點上。接下來將網格節點的帶電荷與外電場進行耦合，解泊松方程，得到各個節點的電場強度。帶電粒子在得到的新的電場的作用下運動，然后利用MCC方法處理離子之間的碰撞過程。最后判斷是否結束，沒有結束進入下一個時間步長，直到循環結束，并進行結果分析。

圖1 Matlab實現模型的流程圖

1.1 滿足麥克斯韋初始速度分布的實現

按照前面提到的模擬流程，第一步要實現參數的初始化。在Ar氣放電的開始階段，粒子滿足麥克斯韋速度分布，需要對粒子速度進行初始化即速度分布滿足方程為

式中，T為氣體溫度；m為電子質量；v為粒子速度；k為玻爾茲曼常數。

圖2 統計結果與麥克斯韋分布的比較

Matlab具有產生滿足一定分布隨機數的內部函數，例如：滿足正態分布隨機數的產生，平均分布隨機數的產生等，但是沒有產生滿足麥克斯韋分布的隨機數的函數。基于式（1）麥克斯韋分布函數，采用舍選法[11]，編程實現滿足麥克斯韋分布函數的隨機速度的生成，并編程對隨機數進行驗證，統計在模擬壓力和溫度的條件下，隨機數的分布情況。對200個隨機數的統計結果進行八次多項式擬合，圖2中的擬合結果與麥克斯韋分布曲線一致，驗證了編程的正確性，說明利用此方法，可以完成粒子速度大小的初始化。

1.2 玻爾茲曼方程離散之PIC部分

玻爾茲曼方程是氣體動力學基本方程，PIC/MCC方法是基于玻爾茲曼離散方程的一種模擬方法。根據文獻[1-2]玻爾茲曼方程可以離散為粒子部分和碰撞部分。其中粒子部分包括電荷分配模型、電場模型以及運動模型。碰撞部分用MCC模型進行研究。

在粒子賦予初值之后，利用權重函數（面積權重）進行二維電荷分配，實現電荷向網格節點的分配。圖3為權重向節點分配的示意圖，利用式（2）[1]計算節點A的面積權重。圖3中（Xp,Yp）為粒子點的位置坐標，正方形單元網格的邊長AB=Δh。與式（2）類似可以得到B、C、D點的權重式。

郜教授提出，乘的本質是縮放。讀者可反復讀“乘的本質是縮放”這一段，對乘法會有新的認識。張景中院士的文章也提出，現代數學中，并不一定從加法出發引入乘法，“幾個幾”的說法，不是數學語言，而是直觀描述，從集合概念出發，建立乘法的概念。然而對于小學生來說，直觀可描述建立乘法概念也是很好的，作為教師當然不應該只滿足于這些，所以郜教授給我們提供的解讀也是非常必要的，對于解讀教材與教學設計應是有幫助的。郜教授的解讀說明教師掌握本體性知識的重要性，提醒教師看教材要有一種居高臨下的能力。

圖3 面積權重分配示意圖

根據網格節點的電荷密度，利用二維泊松方程實現電場求解，具體方法是有限差分法。方程的離散形式為

式中，ρi,j是網格節點（Xi,Yj）（任意節點）處的電荷密度。在泊松方程的求解過程中，每個步長內用200步就可以實現網格節點電勢的求解（具體內容在本文第三部分論述）。由泊松方程解得網格節點的電勢分布，電勢向網格節點線性差值，可以實現節點電場強度的求解。在編程的過程中，需要注意邊界線性差值公式與中間節點的區別。每一個時間步長求解泊松方程，會降低程序的運行效率，但根據相關文獻，采取每一步長內對方程求解，能提高方程的收斂性。

解泊松方程得到網格節點電場的分布，帶電粒子在極板電場和帶電粒子的耦合電場作用下運動。運動方程的離散為

式（4）、式（5）分別是位移和速度在x方向上的離散，而y方向以及vy,vz的求解，與x方向類似，其中p代表第p個粒子。這里是耦合電場，在計算vx時需要通過解得新一個步長的電場分布，新一個步長的電場的求解需通過式（4）計算出新的粒子位置，進而算出t+Δt時刻電場的大小，再利用式（5），得到新時間步長下的速度分布。這樣便實現了電場與運動方程的耦合。

1.3 玻爾茲曼方程離散之MCC模型

氣體在放電的過程中，離子會發生劇烈的碰撞。碰撞形式有：離子或電子與分子間的彈性碰撞、電荷交換碰撞、電離碰撞、激發碰撞等。在電場中，發生何種形式的碰撞可以用MCC方法進行判斷。MCC模型假設碰撞發生在一個時間步長的結束時刻，這樣對于編程是可以實現的。對粒子先后用PIC和MCC模型處理，便實現了粒子模型與碰撞模型的耦合[1-2,12-13]。在利用Matlab編程的過程中需注意碰撞過程中電子的消失與出現的問題，可以利用Matlab提供的空函數和賦值函數來實現。

MCC模型在一個時間步長內，一個粒子最多只能處理一次碰撞，這就需要考慮時間步長的選取不能太大。參考birdsall提出的方法[1,12]，選取時間步長Dt=10-8s。

在放電過程中發生碰撞的概率，采用以下式確定：

式中，n為粒子密度；v為速度；ε為電子能量；σi為總碰撞截面，是各種類型的碰撞截面的總和。碰撞截面是與能量有關的量，隨著粒子能量的變化而變化。意味著，粒子在不同時刻，所具有的總能量不同，決定發生碰撞的概率不同。

根據動量和能量守恒定律，中子和電子碰撞前的速度分別為V, v，碰撞后的速度分別為V', v'。當處理電子與中性粒子發生碰撞時，假設M+m≈M，V'≈V，即二者質量不在同一數量級，忽略中性粒子碰撞前后的速度差。電子碰撞后的速度可由下式計算得到[12-13]：

式中，h為笛卡爾分量；χ為散射角。散射角的計算用各項同性散射假設，參考式cosχ=1-2R2，R2為[0, 1]之間的隨機數。關于散射角的處理可以參考其他論文的方法[1-2]。

2 模擬條件與結果分析

模擬條件：上表面加100V的直流電壓，極間距為1mm，極板寬度為1mm；初始粒子權重為1010；氣體壓力：1Pa；初始時刻溫度：300K。其余3個壁面接地，進行模擬。

模擬氬氣放電過程，只考慮電子與中性粒子之間的彈性碰撞、電離碰撞和激發碰撞。經過100個時間步長之后，得到如圖4所示的電子圖4（a）和離子圖4（b）位置分布，與實拍實驗照片圖4（c）對比。由圖可以看到，在陽極附近區域有很多電子積累，電子的密度也呈現有一定規律的不均勻分布。

圖4 氬氣放電過程模擬

圖4（a）為電子的位置分布，可以看到明顯的陽極鞘層的形成。在放電初始階段，由于電子的質量比離子質量小的多，在電場的作用下，高速向陽極運動被陽極大量吸收。離子的滯后性，會在電場停留積累。隨著時間的推移，電子被吸收的越多，留下的正離子也越多，這時離子與陰極板形成的電場會中和掉一部分的極板電場，在放電區域與陰極板間會形成電壓梯度很大的鞘層區域，稱為陰極鞘層。由于離子運動速度比電子慢的多，甚至會在陽極附近出現離子正電勢大于基板正電勢的情況，這樣便在陽極附近形成與陰極類似的區域，稱為陽極鞘層，與圖4（c）比較，可以看出電子分布與實驗圖片吻合。

圖4（b）為離子的位置分布。經過一段時間的運動，離子被陰極吸附。由于陰極鞘層的作用，陽極會發生二次電子發射，來維持等離子體區域的電子平衡，這是等離子體形成的重要過程。本文模擬沒有考慮二次電子發射的影響，這是模型需要改進的地方。實驗照片圖4（c）中，可以明顯看出陰極板上方有很明亮的放電區域，即為陰極鞘層。此現象與圖4（b）中離子分布的結果一致。

在二維模擬電子和離子分布的情況下，受邊界0電荷條件的限制，電場分布并不是理想的直線分布。分布在不同位置的電子所受到的電場力是不同的。速度和位置的不同，決定粒子能量不同，發生碰撞的概率也不相同。因此二維模擬相比與一維更能反映真實電場中的情況。

取模擬步長的前50步，得到模擬電子數在前50個步長中的變化趨勢，如圖5所示。模擬區域中電子數目是反映放電激烈程度的主要參數。

在模擬開始階段電子數目有一些波動，這是電子運動到陽極板和粒子電離出電子的雙重作用的結果。在模擬進行到一段時間后，由于陰極附近離子的增加，被陰極吸收的電子數將顯著減少，這時原子電離占主要作用，電子數目迅速增加。此結果與相關論文[9]中的模擬結果是一致的。

圖5 電子數隨時間步長的變化曲線

根據步長結束時電勢分布圖（如圖6所示），可以看出，電勢分布在二維區域中的分布是光滑的曲面。與300個循環的結果對比發現，200個循環解本模型中的泊松方程是足夠的。并且能夠明顯地看到邊界電勢與中間部位的電勢分布趨勢是不一樣的。中間區域變化較平滑，即電場強度比較接近。而在放電區域邊緣，電勢變化較劇烈，即場強變化較大。再一次驗證了二維模擬與一維模擬是有顯著區別的，故不能忽略二維邊界條件對放電過程的影響。

圖6 步長結束時電勢分布圖

對速度均方根處理，得到50個時間步長結束時，速度的對數分布示意圖如圖7所示。由于速度與溫度的正相關性，速度分布可以反映出溫度分布。由圖可以看出在距離陽極板0.4mm左右，橫向中間的位置，半徑為0.13mm的區域內溫度最高，電子

圖7 速度的對數（溫度）分布

運動速度最快，放電最激烈。所以陽極板容易受到周圍放電反應的影響。此結果與氬氣放電反應的實驗相吻合，也驗證了模型的正確性。

3 結論

本文基于Matlab編程實現了二維氣體放電的PIC/MCC模型，并對直流低壓條件下氬氣放電過程進行了模擬驗證：①編程實現麥克斯韋隨機速度分布，實現粒子速度的初始化。其次，通過觀察100個時間步長后的二維電子與離子位置分布圖，可以看到陽極鞘層的形成，此現象與實際氬氣放電過程相符。模擬過程考慮了邊界條件對模型的影響，相較與一維模擬，該方法更接近真是的實驗條件；②比較200個循環結束時的電勢分布與300個循環時的泊松方程的解，證明泊松方程200個循環的求解是足夠的。從此，根據模擬結果50個時間步長結束時超電子數達到110，并有直線增加的趨勢，可以說明電子崩的發生。最后，根據溫度分布示意圖可以得到溫度最高區域在陽極附近，此結果與氣體放電實驗相吻合。至此，可以證明通過Matlab編程基本實現了PIC/MCC模型。

[1] Birdsall C K. Particle-in-Cell Charged-Particle simulations, plus Monte Carlo collision with neutral atoms, PIC-MCC[J]. IEEE Trans.Plasma Sci, 1991, 19(2): 65-85.

[2] Nanbu K. Probability theory of electron-molecule, ion-molecule, molecule-molecule, and coulomb collisions for particle modeling of materials processing plasma and gases[J]. IEEE Trans. on Plasma Science, 2000, 28(3): 971-989.

[3] 汪沨, 肖曉林, 張憲標, 等. 基于PIC法SF6/N2混合氣體中絕緣子沿面放電特性研究[J]. 電工技術學報, 2011, 26(8): 220-226.

[4] 吳變桃, 肖登明. 用改進的蒙特卡羅法模擬SF6和CO2混合氣體電子崩參數[J]. 電工技術學報, 2007, 22(1): 13-16.

[5] 李靜, 曹云東, 王爾智, 等. 斷路器短間隙氣體擊穿過程的粒子模擬[J]. 中國電機工程學報, 2010(16): 125-130.

[6] 汪沨, 李錳, 潘雄峰, 等. 基于FEM-FCT算法的SF6/N2混合氣體中棒-板間隙電暈放電特性的仿真研究[J]. 電工技術學報, 2013, 28(9): 261-267.

[7] 馮明, 張鑒, 黃慶安. ICP刻蝕中等離子體分布的模擬[J]. 電子器件, 2006, 28(3): 529-531.

[8] 張志涌. MATLAB教程[M]. 5版. 北京: 北京航空航天大學出版社, 2003.

[9] 李靜, 曹云東, 鄒積巖, 等. 直流氣體放電過程中的離子分子碰撞模型[J]. 高電壓技術, 2009(7): 1677-1682.

[10] 邵福球. 等離子體離子模擬[M]. 北京: 科學出版社, 2002.

[11] 劉沛華, 魯華祥, 龔國良, 等. 基于FPGA的高速任意分布偽隨機數發生器[J]. 應用科學學報, 2012, 30(3): 306-310.

[12] Nanbu K. Particle modeling of nonequilibrium plasmas and gases for materials processing[J]. Vacuum Science Technology, 2004(5): 1-57.

[13] Birdsall C K. Plasma physics via computer simulation[M]. New York: McGraw-Hill, 1985.

Realization of Two-Dimensional PIC/MCC Model based on Matlab

Wang Junjie1Zheng Jinhua1,2Wei Xinxu1Wu Shuang1Xu Lu1
(1. Engineering Research Center of Energy-saving Technology & Equipment of Thermal Energy System, Ministry of Education, School of Chemical Engineering and Energy, Zhengzhou University, Zhengzhou 450001;
2. Zhengzhou Reborn Petroleum Machinery Co., Ltd, Zhengzhou 450001)

In order to achieve the two-dimensional simulation of gas discharge under limited conditions, a simulation program based on Matlab has been developed in the paper, and utilized to verify the discharge process of argon under direct current and low pressure. Electronic speed was initialized through Maxwell equations and the results were verified. The occurrence of electron avalanche in the process of argon discharge was verified by extracting the change trend of the number of electrons in the 50 time steps, and the production of anode sheath was observed through the distribution of electron position at the end of the 100 time step. According to the solution of Poisson’s equation to the 200 step and 300 step, the solution of Poisson's equation in the 200 step was sufficient under the simulated conditions in this paper. The diagram of temperature distribution given by the root mean square treatments of the electron velocity successfully verified the experimental phenomena that the temperature near the anode plate was higher in the DC argon discharge. Thus, the two-dimensional PIC/MCC program based on Matlab can satisfactorily simulate the gas discharge.

Matlab; PIC/MCC; DC discharge; mathematical model; maxwell distribution

王俊杰（1992-），男，碩士研究生，研究方向為等離子體放電的模擬研究。