基于二維PIC-DSMC的微間隙氣體放電等離子體演化過程研究*

徐　翱，尚紹環(huán)，金大志，談效華

（中國工程物理研究院電子工程研究所，四川綿陽621900）

基于二維PIC-DSMC的微間隙氣體放電等離子體演化過程研究*

徐翱*，尚紹環(huán)，金大志，談效華

（中國工程物理研究院電子工程研究所，四川綿陽621900）

準(zhǔn)確的理解微間隙氣體放電中非平衡等離子隨時(shí)間的演化過程對(duì)于設(shè)計(jì)氣體開關(guān)、微電子及其它等離子體器件有著非常大的幫助。通過二維PIC-DSMC耦合算法模擬了一個(gè)大氣壓氮?dú)猸h(huán)境下微間隙平板電極發(fā)生氣體放電時(shí)電子及離子的運(yùn)動(dòng)演化過程，得到了氣體放電過程中平板電極間電子和離子數(shù)密度分布隨時(shí)間變化的趨勢，討論了陽極附近電子云的形成與演變、陰極附近存在的鞘層以及電子和離子的速度及溫度分布，最后將模擬計(jì)算得到的擊穿電壓與帕邢曲線及相關(guān)實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比，為相關(guān)等離子體器件的進(jìn)一步發(fā)展提供了理論基礎(chǔ)。

放電；網(wǎng)格質(zhì)點(diǎn)法；直接模擬蒙特卡羅法；等離子體

隨著近些年來等離子體計(jì)算模型及診斷技術(shù)的不斷發(fā)展，直流或脈沖條件下的微間距氣體或真空擊穿研究受到了研究人員的廣泛關(guān)注［1-5］。但是由于微間隙的極間距一般在微米量級(jí)，通過傳統(tǒng)的等離子體診斷方法能夠得到的相關(guān)信息仍然非常有限，而電學(xué)測量方法只能得到擊穿過程的電壓與電流關(guān)系，因此通過粒子模擬的方法獲得微間隙氣體放電過程的行為特性就成為了相關(guān)研究的重點(diǎn)。Zhang、Lee和Radmilovic-Radjenovic等人［6-8］首先使用PIC方法實(shí)現(xiàn)了微間隙氣體放電的數(shù)值模擬，并將計(jì)算得到的擊穿電壓與實(shí)驗(yàn)結(jié)果進(jìn)行了比較。2012年美國Sandia國家實(shí)驗(yàn)室的Moore等人在PIC的基礎(chǔ)上引入了PIC-DSMC混合算法實(shí)現(xiàn)了對(duì)大氣條件下微間隙平板電極氣體放電的一維模擬［9］，得到了更多的氣體放電行為特性。2013年Moore等人又利用一維PIC-DSMC算法分析了微間隙平板電極中不同背景中性氣體密度對(duì)放電的影響［10］。但目前的這些研究基本都集中在如何實(shí)現(xiàn)數(shù)值模擬結(jié)果與帕邢曲線或?qū)嶒?yàn)測得的擊穿電壓相符，對(duì)微間隙氣體放電過程中陽極附近電子云的形成與演變、陰極附近的鞘層、間隙中正離子及電子的運(yùn)動(dòng)導(dǎo)致的等離子體通道形成等過程并沒有詳細(xì)的報(bào)道，而這些卻是能夠證實(shí)數(shù)值模擬能夠真正解釋微間隙氣體放電物理現(xiàn)象的重要部分。因此本文基于二維PIC-DSMC耦合算法模擬了一個(gè)大氣壓氮?dú)猸h(huán)境下微間隙平板電極發(fā)生氣體放電時(shí)電子及離子的運(yùn)動(dòng)行為和等離子體通道的形成過程，并對(duì)其中等離子體的相關(guān)特性進(jìn)行了分析，為相關(guān)等離子體器件的進(jìn)一步發(fā)展提供了堅(jiān)實(shí)的基礎(chǔ)。

1　計(jì)算模型

PIC-DSMC（網(wǎng)格質(zhì)點(diǎn)法-直接模擬蒙特卡羅法）耦合算法是一種可以較好的模擬等離子體運(yùn)動(dòng)演化過程的數(shù)值模擬方法。本文模擬使用的PICDSMC耦合算法由PIC模塊與DSMC模塊組成，其中PIC模塊處理所有帶電粒子的運(yùn)動(dòng)、電場的計(jì)算、粒子與電極的作用、粒子間發(fā)生的電離、復(fù)合等化學(xué)反應(yīng)，而DSMC處理中性粒子的運(yùn)動(dòng)，這樣就能夠滿足氣體擊穿過程包含電子、離子和中性粒子等多種類復(fù)雜運(yùn)動(dòng)的模擬需要。

實(shí)際物理模型是間隙為d的兩個(gè)對(duì)稱的圓形平板電極，因此可以用二維圓對(duì)稱模型來進(jìn)行模擬，即只需模擬兩電極中心對(duì)稱軸與電極之間的二維區(qū)域來表征整個(gè)物理模型。假設(shè)模擬區(qū)域的電極半徑r=20 μm，電極間隙d=50 μm，模擬模型如圖1所示。

圖1　模擬模型

模擬初始時(shí)刻兩電極間充滿均勻分布的常溫（20℃）氮?dú)猓鋲簭?qiáng)為1個(gè)標(biāo)準(zhǔn)大氣壓，另外此時(shí)電極間還存在均勻分布的游離自由帶電粒子，假設(shè)其為電子和氮?dú)怆x子，數(shù)密度都為1×1015m-3。文獻(xiàn)［10-11］中一般要求設(shè)置初始電子離子數(shù)密度小于中性粒子密度的1/109，通過反復(fù)計(jì)算可以發(fā)現(xiàn)微間隙氣體擊穿模擬中初始電子離子數(shù)密度的設(shè)置值在1015m-3量級(jí)附近，變化±1個(gè)量級(jí)對(duì)計(jì)算得到的擊穿電壓值基本沒有影響，僅僅對(duì)等離子體通道形成時(shí)間有一定影響，因此設(shè)置初始電子離子密度為1×1015m-3）。同時(shí)使得陰極接地，陽極電壓設(shè)為570 V，考慮多種電子、離子及中性粒子間的碰撞及激發(fā)、電離、復(fù)合等化學(xué)反應(yīng)，如表1所示，所使用的相關(guān)碰撞及反應(yīng)截面數(shù)據(jù)見文獻(xiàn)［12-14］，而中性粒子間的彈性碰撞由DSMC模塊直接計(jì)算。

表1　碰撞及化學(xué)反應(yīng)

PIC的網(wǎng)格及時(shí)間步長的選取原則［15］一般為網(wǎng)格Δx≤λD，其中λD表示德拜長度；而時(shí)間步長Δt≤0.2/ωpe，其中ωpe表示等離子體角頻率。考慮到發(fā)生氣體擊穿時(shí)的數(shù)密度一般在1020m-3量級(jí)左右，而電子溫度在幾eV到幾十eV量級(jí)，因此模擬時(shí)選擇PIC模塊的網(wǎng)格為邊長0.5 μm的正方形網(wǎng)格，時(shí)間步長為1×10-13s。同時(shí)由于DSMC的網(wǎng)格及時(shí)間步長選取原則［16］一般為網(wǎng)格Δx≤λ/3，其中λ表示平均自由程；時(shí)間步長Δt≤Δx/vmax，其中vmax表示該方向粒子的最大速度，而發(fā)生氣體擊穿時(shí)平均自由程一般都大于德拜長度（假設(shè)電子溫度為1 eV，則模擬條件下的電子平均自由程約4.5 μm，而對(duì)應(yīng)的德拜長度為約0.74 μm），因此可以設(shè)置DSMC模塊與PIC模塊使用同一套網(wǎng)格，但DSMC的時(shí)間步長設(shè)置為1×10-12s，即每運(yùn)行10次PIC程序，運(yùn)行1次DSMC程序。而考慮發(fā)生擊穿時(shí)仿真分子的數(shù)量為每個(gè)網(wǎng)格有約80個(gè)代表電子、80個(gè)代表離子和40個(gè)代表中性粒子，則可以設(shè)置最大仿真分子數(shù)為80萬個(gè)。并且每個(gè)仿真分子代表的真實(shí)粒子權(quán)重在程序中可以自動(dòng)調(diào)整。

模擬考慮的兩個(gè)電極均為清潔表面的銅電極，假設(shè)中性粒子與電極表面的作用為反射率為1的漫反射，而使用的電子轟擊電極表面產(chǎn)生的二次電子率曲線見參考文獻(xiàn)17，離子轟擊電極表面產(chǎn)生的二次電子率由下式得到［18］：

其中Eiz，s表示入射離子的電離能，Eφ表示以eV為單位的功函數(shù)。最后在模擬過程中監(jiān)控電極間的傳導(dǎo)電流，如果發(fā)現(xiàn)傳導(dǎo)電流迅速增加則認(rèn)為等離子體通道已經(jīng)形成，氣體擊穿即將發(fā)生，模擬計(jì)算中止。

2　微間隙氣體放電等離子體演化過程分析

微間隙氣體放電的物理特征就是一個(gè)在微間隙中形成等離子體通道并導(dǎo)通的過程，因此通過PIC-DSMC耦合算法計(jì)算發(fā)生氣體擊穿時(shí)電子、N2+及N+的數(shù)密度分布隨時(shí)間變化就可以得到微間隙氣體放電的演化過程。假設(shè)電極加電壓的瞬間為初始時(shí)刻t=0，則計(jì)算得到的t=1 ns、3 ns、5 ns、8 ns、10 ns和12 ns時(shí)電子、N2+及N+的數(shù)密度在整個(gè)模擬計(jì)算區(qū)域中的分布云圖如圖2～圖4所示，這些圖的左側(cè)為陰極，右側(cè)為陽極，下端為對(duì)稱軸，與圖1所示一致。

圖2　電子數(shù)密度分布云圖

通過圖2可以發(fā)現(xiàn)，t=1 ns時(shí)陽極附近形成電子云，這應(yīng)該是由于間隙中的游離電子受陽極的吸引向陽極附近運(yùn)動(dòng)并且獲得了一定的能量，再由于電子轟擊陽極產(chǎn)生的二次電子，這會(huì)造成在陽極附近出現(xiàn)許多電子并且大量電離背景氣體產(chǎn)生更多的電子，從而形成陽極附近的電子云；而隨著時(shí)間增加到t=5 ns左右時(shí)陽極附近的電子數(shù)密度、電子分布區(qū)域都呈現(xiàn)減小的趨勢，這是由于電子不斷向陽極運(yùn)動(dòng)，而通過電子電離背景氣體產(chǎn)生的電子數(shù)量少于被陽極俘獲的電子數(shù)量，從而造成電子數(shù)密度及電子分布區(qū)域都減小；但隨著時(shí)間增加到t=8 ns左右時(shí)陽極附近的電子數(shù)密度繼續(xù)減小，但陰極附近出現(xiàn)電子并且電子分布區(qū)域出現(xiàn)增大的趨勢，這是由于正離子轟擊陰極產(chǎn)生了二次電子，這些二次電子又迅速向陽極運(yùn)動(dòng)且運(yùn)動(dòng)過程中又電離產(chǎn)生了新的電子從而造成電子分布區(qū)域增大，另外此時(shí)緊挨陽極的附近出現(xiàn)高電子數(shù)密度區(qū)域，這可能與電子轟擊陽極產(chǎn)生的二次電子又被陽極吸引從而又轟擊陽極產(chǎn)生新二次電子這樣的循環(huán)過程有關(guān)；最后時(shí)間增加到t=12 ns左右時(shí)電子已經(jīng)能夠形成通道，這是由于正離子轟擊陰極產(chǎn)生二次電子，二次電子向陽極運(yùn)動(dòng)時(shí)又不斷產(chǎn)生新的電子和正離子，而這些新的正離子會(huì)向陰極運(yùn)動(dòng)并轟擊陰極產(chǎn)生二次電子，如此循環(huán)從而導(dǎo)致電子通道的形成。

圖3　N2+數(shù)密度分布云圖

通過圖3可以發(fā)現(xiàn)t=1 ns時(shí)陽極附近同樣會(huì)形成N2+云，這說明陽極附近確實(shí)存在大量的電離；而隨著時(shí)間增加到t=5 ns左右時(shí)陽極附近的N2+云的高數(shù)密度中心出現(xiàn)向陰極方向的移動(dòng)，這是由于受陰極的吸引向陰極運(yùn)動(dòng)，但N2+的速度比電子速度慢，因此其數(shù)密度分布變化不是特別明顯；但隨著時(shí)間增加到t=8 ns左右時(shí)N2+云的高數(shù)密度中心出現(xiàn)明顯擴(kuò)散開的趨勢，靠近陰極附近的N2+數(shù)密度也明顯增加，而陽極附近卻出現(xiàn)顯著的N2+數(shù)密度減小，這些都與電子數(shù)密度分布的變化契合，進(jìn)一步表明了正離子轟擊陰極產(chǎn)生二次電子，這些二次電子迅速向陽極運(yùn)動(dòng)且運(yùn)動(dòng)過程中又電離產(chǎn)生了新的電子和正離子從而導(dǎo)致區(qū)域內(nèi)等離子體密度增大，同時(shí)陽極附近的電子轟擊陽極產(chǎn)生二次電子被吸引轟擊陽極又產(chǎn)生新的二次電子這樣的循環(huán)過程中電子難以獲得足夠的能量來電離背景氣體，所以正離子數(shù)量明顯減小；最后時(shí)間增加到t=12 ns左右時(shí)N2+已經(jīng)在微間隙內(nèi)形成了明顯的通道。

圖4　N+數(shù)密度分布云圖

對(duì)比圖3與圖4可以發(fā)現(xiàn)，N+的數(shù)密度分布隨時(shí)間變化趨勢與N2+比較類似，但是N+的變化稍稍快一些，這是由于其質(zhì)量較小，因此變化速度相對(duì)較快。另外N+的數(shù)密度明顯小于N2+，這是由于背景氣體主要是N2，而N需要由N2通過離解等反應(yīng)才能產(chǎn)生，因此相對(duì)較少，再將N電離才能產(chǎn)生N+的數(shù)量就更少。

通過分析1 ns、6 ns和12 ns時(shí)靠近對(duì)稱軸線區(qū)域（由于所使用的軟件無法提取對(duì)稱軸線上的相關(guān)數(shù)據(jù)，因此選擇非常靠近對(duì)稱軸線的區(qū)域，即r=0.3 μm來代表對(duì)稱軸區(qū)域）的電場及電勢變化，如圖5和圖6所示，可以發(fā)現(xiàn)陰極附近（靠近x=0）的電場和電勢都隨著時(shí)間t的增加而增大，并且其隨著x增大而變化的趨勢也更加劇烈；而陽極附近（靠近x=50 μm）的電場隨著時(shí)間t的增加而減小，電勢隨著時(shí)間t的變化不大，但其隨著x減小而變化的趨勢變得相對(duì)平緩。這說明了陰極附近應(yīng)該會(huì)形成鞘層，并且該鞘層隨著時(shí)間t增加而變得更加明顯。而陽極附近并沒有形成明顯的鞘層，這可能是由于陽極附近隨著時(shí)間t增加形成了靠近陽極的高數(shù)密度電子云，但該電子云沒有能夠電離出足夠多的正離子，從而使得該區(qū)域變成了“電子通道”而不是鞘層。

圖5　r=0.3 μm處的電場分布

圖6　r=0.3 μm處的電勢分布

通過分析12 ns時(shí)靠近對(duì)稱軸線區(qū)域（r=0.3 μm）的電子速度及N2+速度，如圖7所示，可以發(fā)現(xiàn)電子速度大于N2+速度2-3個(gè)數(shù)量級(jí)，這證明了離子運(yùn)動(dòng)確實(shí)比電子運(yùn)動(dòng)慢；另外電子速度在陰極附近有明顯突增，然后從陰極向陽極稍有減小，這可能是由于電子離開陰極后立即受到強(qiáng)靜電場作用使其速度突增，但隨后靜電場減小并且電子會(huì)不斷的與氮?dú)獍l(fā)生碰撞從而又使其速度有所降低；而N2+主要是由電子碰撞電離產(chǎn)生，產(chǎn)生后就會(huì)受到靜電場作用不斷使其加速，又由于其仍然會(huì)與氮?dú)獾扰鲎矎亩鴾p速，因此其在大部分區(qū)域的速度都比較相似，但隨著其靠近陰極，靜電場會(huì)明顯增加，因此靠近陰極的N2+速度明顯增大。而通過12 ns時(shí)靠近對(duì)稱軸線區(qū)域（r=0.3 μm）的電子溫度及N2+溫度，如圖7所示，可以發(fā)現(xiàn)電子溫度同樣大于N2+溫度1個(gè)數(shù)量級(jí)，并且它們的分布趨勢與速度分布趨勢類似。

圖7　r=0.3 μm處的電子速度和N2+速度

圖8　r=0.3 μm處的電子溫度和N2+溫度

模擬計(jì)算得到的電極上收集到的電流與時(shí)間t的關(guān)系如圖9所示，可以發(fā)現(xiàn)t=12 ns左右時(shí)電流會(huì)迅速增加，因此可以認(rèn)為等離子體通道形成，微間隙內(nèi)將會(huì)發(fā)生放電。而基于相似的初始條件設(shè)置可以模擬計(jì)算得到不同間距微間隙的擊穿電壓，并與帕邢曲線計(jì)算得到的擊穿電壓及參考文獻(xiàn)4實(shí)驗(yàn)測試得到的擊穿電壓相比，如圖10所示，可以發(fā)現(xiàn)模擬計(jì)算得到的擊穿電壓比帕邢曲線及實(shí)驗(yàn)測試得到的擊穿電壓稍高一些，但曲線變化趨勢基本相似，造成這樣的原因可能是由于模擬中沒有考慮到激發(fā)態(tài)原子退激后產(chǎn)生的光輻射會(huì)造成光致電離，以及離子轟擊電極產(chǎn)生二次電子率公式比較簡單等。但是仍然可以認(rèn)為通過PIC-DSMC耦合算法能夠較為準(zhǔn)確的模擬計(jì)算得到微間隙氣體放電的物理過程。

圖9　電流與時(shí)間的關(guān)系

圖10　擊穿電壓與間距的關(guān)系

3　結(jié)論

通過與帕邢曲線及實(shí)驗(yàn)結(jié)果比較，可以認(rèn)為基于PIC-DSMC耦合算法能夠較為準(zhǔn)確的獲得氣體放電特性曲線，因此其模擬得到的放電形成過程中電子和離子的運(yùn)動(dòng)演化過程及相關(guān)等離子體參數(shù)時(shí)空分布規(guī)律能夠反映真實(shí)氣體放電過程的情況，所以通過對(duì)這些氣體放電等離子體演化過程研究就可以為設(shè)計(jì)氣體開關(guān)、微電子及其它等離子體器件打下堅(jiān)實(shí)的理論基礎(chǔ)。

［1］Hitchon W N.The Time History of Reakdown［J］.J Phys D：Appl Phys，2008，41（22）：222002.

［2］Sshnyder R，Howling A A，Bommottet D，et al.Direct Current Breakdown in Gases for Complex Geometries from High Vacuum to AtmosphericPtrssure［J］.JPhysD：ApplPhys，2013，46（28）：285205.

［3］Venkattraman A.Generalized Criterion for Thermo-Field Emission Driven Electrical Breakdown of Gases［J］.Applied Physics Letters，D2014，104（19）：194101.

［4］Radmilovic-Radjenovic M，Radjenovic B，et al.The Breakdown Phenomena in Micrometer Scale Direct-Current Gas Discharges［J］.Plasma Chem Plasma Process，2014，34：55-64.

［5］盧新培，嚴(yán)萍，任春生，等.大氣壓脈沖放電等離子體的研究現(xiàn)狀與展望［J］.中國科學(xué)：物理學(xué)力學(xué)天文學(xué)，2011，41（7）：801-815. Lu Xinpei，Yan Pin，Ren Chunsheng，et al.Review on Atmospheric Pressure Pulsed DC Discharge.Scientia Sinica Phys.Mech& Astron，2011，41（7）：801-815.

［6］酉小廣，成永紅，孟國棟，等.微機(jī)電系統(tǒng)微尺度間隙的脈沖擊穿特性［J］.強(qiáng)激光與粒子束，2013，25（7）：1867-1872. You Xiaoguang，Cheng Yonghong，Meng Guodong，et al.Pulse Breakdown Characteristics of Microscale Gap in Micro-Electro-Mechanical System.High Power Laser and Particle Beams，2013，25（7）：1867-1872.

［7］Zhang W，F(xiàn)isher T S，Garimella S V.Simulation of Ion Generation and Breakdown in Atmospheric Air［J］.Journal of Applied Physics，2004，96（11）：6066-6072.

［8］Lee S M，Seo Y S，Lee J K.Paschen Breakdown Curve by One-Dimensional PIC-MCC Simulation［J］.Computer Physics Communications，2007，177：132.

［9］Radmilovic-Radjenovic M，Radjenovic B.A Particle-in-Cell Simulation of the Breakdown Mechanism in Microdischarges with an Improved Secondary Emission Model［J］.Contrib.Plasma Phys.，2007，47（3）：165-172.

［10］Moore C H，Hopkins M M，Crozier P S，et al.1D PIC-DSMC Simulations of Breakdown in Microscale Gaps［C］//28th International Symposium on Rarefied Gas Dynamics，2012：629-636.

［11］Moore C H，Hopkins M M，Boerner J J，et al.PIC Simulation of Microscale Breakdown in Gaps Having a Non-Uniform Background Neutral Gas Density［C］//31th ICPIG，2013：10-13.

［12］Yukikazu I.Cross Sectrions for Electron Collisions with Nitrogen Molecules［J］.J Phys Chem Ref Data，2006，35（1）：31-53.

［13］Phelps A，Pitchford L.Anisotropic Scattering of Electrons by N2and Its Effect on Electron Transport［J］.Phys Rev A，1985，31（5）：2932-2949.

［14］Straub H C，Renault P，Lindsay B G，et al.Absolute Partial Cross Sectrions for Electron-Impact Ionization of H2，N2and O2from Threshold to 1 000 eV［J］.Physical Reviwe A，1996，54（3）：2146-2153.

［15］Tskhakaya D，Matyash K，Schneider R，et al.The Particle-In-Cell Methode［J］.Contrib Plasma Phys，2007，47（8-9）：563-594.

［16］Bird G A.Molecular Gas Dynamics and the Direct Simulation of Gas Flows［M］.Oxford：Oxford University Press，1994：334-369.

［17］Walker C G H，El-Gomati M M，Assa'd A M D，et al.The Secondary Electron Emission Yield for 24 Solid Elements Excited by Primary Electrons in the Range 250 eV～5 000 eV：A Theory/Experiment Comparisom［J］.Scanning，2008，30（5）：365-380.

［18］Lieberman M，Lichtenberg A.Principles of Plasma Discharges and Materials Processing［M］.New Jersey：Wiley，2005：299-320.

徐翱（1982-）男，漢族，湖北武漢人，就職于中國工程物理研究院電子工程研究所，副研究員，博士。主要研究方向?yàn)榈入x子體技術(shù)及電真空器件研制，xuao@caep.cn。

Study on Plasma Evolution in Microscale Gas Discharge Based on Two Dimension PIC-DSMC Code*

XU Ao*，SHANG Shaohuan，JIN Dazhi，TAN Xiaohua
（Institute of Electronic Engineering，China Academy of Engineering Physics，Mianyang Sichun 621900，China）

It is helpful to exactly understanding non-equilibrium plasma evolution in microscale gas discharge for designing gas switch，microelectronics and other plasma device.The gas discharge in microscale gap size under one atmospheric pressure of N2environment is simulated by a PIC-DSMC code.The kinetic evolution process of electron and ions in this gas discharge is presented，and the distribution of electron density and ion density varied with discharge time is obtained.Then，the formation and change of the electron cloud near anode，the sheath near cathode，and the distribution of velociy and temperture are discussed.Finally，the simulation result is compared with Paschen curve and experiment data.It can provide theoretic help for furher development of plasma device.

discharge；PIC；DSMC；plasma

O461.1

A

1005-9490（2016）05-1025-06

項(xiàng)目來源：中國工程物理研究院科學(xué)技術(shù)發(fā)展基金項(xiàng)目（2015B0102013）；中國工程物理研究院電子工程研究所科技創(chuàng)新基金項(xiàng)目（S20150807）

2015-10-22修改日期：2015-11-10

EEACC:231510.3969/j.issn.1005-9490.2016.05.001