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

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

基于Matlab的二維PIC/MCC模型的實現(xiàn)

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

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

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

Matlab；PIC/MCC；直流放電；數(shù)學(xué)模型；麥克斯韋分布

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

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

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

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

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

1 PIC/MCC模型的Matlab實現(xiàn)

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

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

圖1 Matlab實現(xiàn)模型的流程圖

1.1 滿足麥克斯韋初始速度分布的實現(xiàn)

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

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

圖2 統(tǒng)計結(jié)果與麥克斯韋分布的比較

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

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

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

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

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

圖3 面積權(quán)重分配示意圖

根據(jù)網(wǎng)格節(jié)點的電荷密度，利用二維泊松方程實現(xiàn)電場求解，具體方法是有限差分法。方程的離散形式為

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

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

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

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

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

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

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

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

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

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

2 模擬條件與結(jié)果分析

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

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

圖4 氬氣放電過程模擬

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

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

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

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

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

圖5 電子數(shù)隨時間步長的變化曲線

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

圖6 步長結(jié)束時電勢分布圖

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

圖7 速度的對數(shù)（溫度）分布

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

3 結(jié)論

本文基于Matlab編程實現(xiàn)了二維氣體放電的PIC/MCC模型，并對直流低壓條件下氬氣放電過程進行了模擬驗證：①編程實現(xiàn)麥克斯韋隨機速度分布，實現(xiàn)粒子速度的初始化。其次，通過觀察100個時間步長后的二維電子與離子位置分布圖，可以看到陽極鞘層的形成，此現(xiàn)象與實際氬氣放電過程相符。模擬過程考慮了邊界條件對模型的影響，相較與一維模擬，該方法更接近真是的實驗條件；②比較200個循環(huán)結(jié)束時的電勢分布與300個循環(huán)時的泊松方程的解，證明泊松方程200個循環(huán)的求解是足夠的。從此，根據(jù)模擬結(jié)果50個時間步長結(jié)束時超電子數(shù)達到110，并有直線增加的趨勢，可以說明電子崩的發(fā)生。最后，根據(jù)溫度分布示意圖可以得到溫度最高區(qū)域在陽極附近，此結(jié)果與氣體放電實驗相吻合。至此，可以證明通過Matlab編程基本實現(xiàn)了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]. 電工技術(shù)學(xué)報, 2011, 26(8): 220-226.

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

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

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

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

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

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

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

[11] 劉沛華, 魯華祥, 龔國良, 等. 基于FPGA的高速任意分布偽隨機數(shù)發(fā)生器[J]. 應(yīng)用科學(xué)學(xué)報, 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-），男，碩士研究生，研究方向為等離子體放電的模擬研究。