大氣壓氬氣介質阻擋放電特性的一維粒子仿真

孫安邦，劉天旭，葉書榮，李昊霖

(西安交通大學電氣工程學院電力設備電氣絕緣國家重點實驗室，陜西 西安 710049)

近年來，低溫等離子體在材料表面改性、污染防治、生物醫學等方面有著越來越多的應用。其中大氣壓低溫等離子體的產生不需要昂貴復雜的真空設備，所以其結構簡單且低成本的特點正在受到人們的重視。介質阻擋放電(dielectric barrier discharge，DBD)是大氣壓低溫等離子體產生的一種常用方式。DBD是一種將絕緣介質插入放電空間的一種放電形式，可以在很寬的氣壓范圍(104Pa至106Pa)與頻率范圍(50 Hz至1 MHz)內啟動[1]，屬于高氣壓下的非熱平衡放電。在大氣壓下其放電外觀分為絲狀放電模式和均勻放電模式，前者存在大量絲狀的放電通道，會產生局部的熱效應，對材料表面產生一定的損傷；而后者放電現象分布在整個空間，類似湯森放電和輝光放電，不會損壞試品，非常適合于材料的處理[2]。研究表明，使用惰性氣體氛圍、施加氣流作用等措施能有效地使介質阻擋放電向著均勻模式轉變[2]。

近年來，Massine等人通過實驗測量以及CCD高速攝像對大氣壓下氬氣和氦氣的介質阻擋放電中的輝光放電現象進行了研究[3]；清華大學王新新從放電理論角度出發討論了DBD兩種模式的成因，并通過實驗實現了氦氣和氮氣的均勻模式放電[2]；空軍工程大學的宋飛龍等人采用光譜診斷法研究大氣壓氬氣介質阻擋放電特性[4]。然而，DBD尺寸小的特點給實驗診斷帶來了一定的困難，而且實驗無法以較高的時間與空間精度揭示放電過程中的關鍵物理過程。因此不少研究者對DBD進行了仿真研究。Hagelaar等人使用流體模型，對大氣壓下氬氣DBD射頻放電進行了仿真工作，揭示了放電從α模式向γ模式的轉變過程[5]；西安交通大學的姚聰偉等人搭建耦合了外電路與化學反應的一維流體動力學模型，對氬氣DBD過程中的物理參量變化與放電過程進行了討論[6]；李平等人建立二維軸對稱流體模型，通過有限元法研究氬氣在不同氣壓下DBD特性[7]。但由于DBD的非熱平衡性，電子能量呈現非麥克斯韋分布，流體模型的描述尚不夠精確。而粒子模擬直接跟蹤帶電粒子的運動，其物理過程是自洽的，可以提供等離子體任何詳細的運動信息[8]。

本文搭建了一維PIC/MCC(Particle-in-cell/Mont Carlo Collision)模型，仿真研究了氬氣在均勻模式下的DBD特性，通過耦合彈性碰撞、激發碰撞、電離碰撞以及復合等反應，實現了對放電過程的精確描述；通過使用粒子權重自適應(APM)以及基于OpenMPI的并行技術加速仿真過程，實現了500kHz正弦交流電源驅動下氬氣DBD過程的仿真，探究了其時空演化特性以及外施電壓、氣隙間距、介質板介電常數等參數對放電特性的影響。本文與參考文獻[6]的流體模型的放電形貌和物理量數值與變化趨勢具有較好的一致性。

1 模型的建立

1.1 PIC/MCC模型

PIC(Particle in cell)模型通過跟蹤宏粒子的運動并將宏粒子向格點的插值來獲得空間物理量的分布，并通過MCC(Monte Carlo Collision)過程，引入粒子的短程碰撞，使PIC算法可以處理粒子彈性碰撞、激發、電離等過程，二者結合就是PIC/MCC模型[8]。算法流程如圖1所示，圖中t為仿真時刻，tend為仿真結束時刻。

圖1 PIC/MCC算法流程

首先進行粒子初始化，之后每個時間步內依次執行累積電荷、求解電磁場、推動粒子、碰撞處理(如碰撞、激發、電離、復合等過程)四個步驟[9]，最后對粒子進行統計平均得到宏觀物理量。

粒子初始化中，本文假定起始時刻的電子密度與離子密度相等，均為8×1015m-3。電子與離子的位置在氣隙空間均勻分布，而速度呈麥克斯韋分布。累積電荷的過程中，本文采用一階插值函數，將帶電粒子的電荷映射到相鄰的兩個格點之上。一維DBD結構的電場求解使用泊松方程：

(1)

其中φ表示電勢，ρ表示電荷密度，ε表示介電常數，e表示元電荷量，ni與ne分別表示該點的離子數密度與電子數密度。DBD過程中，運動到介質板表面的帶電粒子會附著其上，形成表面電荷。介質板表面電荷及其形成的電場對于DBD的演化過程十分重要。介質表面電荷密度σ的變化可用下式表示：

(2)

其中，Γe與Γi分別為介質板處電子與離子的通量，γ為介質表面的二次電子發射系數。

在離散的條件下，泊松方程可以寫成如下格式：

(3)

其中Δx為空間格點的間距，ε為介電常數，i=1,2,...n，由左側極板接電源正極，右側極板接地的邊界條件可知：

(4)

假設左側與右側的介質板與等離子體的邊界的格點序號分別是d1與d2，左側、右側介質板的厚度分別是l1和l2，以左側介質板為例，有：

Dd1+1/2S-Dd1-1/2S=σ1S+ρd1SΔx/2

(5)

其中Dd1-1/2和Dd1+1/2分別是邊界左、右側Δx/2處的電位移矢量，S為包圍路徑法向的面積，σ1為邊界面介質板表面電荷密度，ρd1為緊鄰邊界面空間網格內的電荷密度，用格點電勢表示電位移矢量，有：

(6)

其中ε0和ε1分別為真空和左側介質板的介電常數，考慮兩側極板的邊界條件，有：

(7)

同理，右側介質板與等離子體邊界處的邊界條件為：

(8)

在仿真過程中等離子體區域的電場強度由外部電場與等離子體區域帶電粒子的自洽電場所疊加構成，外部電場由外施電場與介質板表面電荷電場所疊加，介質板表面電荷電場對放電的發生與熄滅有著重要的影響。

假設左右介質板表面電荷面密度分布為σ1、σ2，根據高斯定理可以求得表面電荷單獨作用在等離子體區域產生的電場為：

(9)

假設左、右介質板的厚度為d，其相對介電常數為ε1，氣隙間距為L，加在極板上的外施電壓為V(t)，左側介質板、氣隙、右側介質板中電場強度分別為EL、EV和ER，它們滿足如下關系：

(10)

可以得出外施電壓單獨作用在等離子體區域產生的電場為：

(11)

合成電場EO由外施電壓與介質板表面電荷單獨作用電場的疊加：

EO=Eσ+EV

(12)

帶電粒子的運動采用velocity verlet算法，利用粒子前一時刻的位置xn、速度vn、加速度an與該時刻的加速度an+1，通過迭代得到粒子下一時刻的位置xn+1與速度vn+1[10]：

(13)

蒙特卡洛碰撞過程中考慮的電子、離子碰撞類型如表1所示。采用偽碰撞(Null Collision)[11]和相應的碰撞截面數據σ(ε)[12]來確定具體碰撞的類型。

表1 電子、離子碰撞反應類型

1.2 復合過程

本文中考慮電子與離子的復合過程，定義了粒子源項Si(m-3s-1)，表示單位時間單位體積內有Si個電子-離子對結合變成中性原子[4]。

(14)

fre代表復合反應的總數量，ne代表參與復合反應的電子數密度，ni代表參與復合反應的離子數密度，nα代表參與復合反應的中性分子數密度，而kri代表了第i個復合反應的復合系數。在每個格點Vcell的體積中，δt時間內復合的電子-離子對數量Nr為：

Nr=SiVcellδt

(15)

表2給出本文考慮的復合反應類型及復合系數，反應系數來自文獻[13-17]。

表2 復合反應及系數

在氬氣DBD過程中，存在Ar+和Ar2+兩種離子。Ar+由電離碰撞產生，文獻[6]的研究表明，由于電荷交換碰撞反應速率過高，放電之后每個時刻分子離子Ar2+的密度與總粒子密度的占比要明顯高于原子離子Ar+，說明分子離子Ar2+是正離子的主要成分。說明真實物理世界中，Ar電離產生的Ar+幾乎全部轉化為Ar2+，即Ar2+的數密度與原本產生的Ar+是高度一致的，且二者均是帶著一個正電荷，在本文中，原子離子Ar+的數密度時空分布可以等效被近似為Ar2+的數密度時空分布。

1.3 粒子權重自適應

PIC/MCC算法的程序運行時間與宏粒子數正相關。實際放電工況下帶電粒子數隨時間不斷變動，如果宏粒子數量過多，模型求解速度極慢且消耗的內存巨大；如果宏粒子數量過少則會引入模擬噪聲。本文中采用粒子自適應權重(Adapted particle management,APM)來動態調整粒子的權重以實現仿真區域中合適的粒子數量，常用的粒子合并方法有兩個粒子合并成為一個粒子，即2-1，或者三個粒子合成兩個粒子，即3-2等。本文采用2-1粒子合并算法，需要說明的是，雖然2-1粒子合并算法無法同時保證動量與動能守恒，但是可以通過合并速度近似相等的兩個粒子來盡可能實現結果精確。文獻[17]給出了具體實現過程：

首先設置每個網格內期望的粒子數Npcc，帶電粒子的期望權重為：

(16)

其中n為帶電粒子的數密度。對每個格點的期望權重與粒子實際權重進行比較，當粒子權重w<2/3wd時執行合并操作，首先使用k-d tree算法尋找距離同一網格內速度最接近的兩個同種粒子1和2，合成的新粒子A的權重為兩個舊粒子權重之和，速度隨機選擇合并前的一個粒子的速度，新粒子的位置由下式進行定義：

(17)

粒子權重w>1.5wd時執行分裂操作，將舊的粒子分裂成兩個同種類的新粒子，每個新粒子擁有與舊粒子相同的位置、速度、加速度以及一半的權重。

1.4 并行計算

并行計算是一種能夠在同一時間內執行多條指令的算法，適合快速求解大規模而復雜的計算問題。目前國內外的高性能計算機系統中，并行編程環境OpenMPI得到了廣泛的應用[18]。在PIC/MCC算法中，需要逐個處理每個宏粒子的運動、碰撞等過程，適合進行并行計算處理。本文將電子、離子均勻分配到多個核心同步處理，而電荷向格點的映射以及電磁場的求解等工作則由單一核心完成，以縮短程序運行時間。

1.5 DBD模型參數與仿真參數設置

本文中介質阻擋放電的仿真區域為平行平板電極結構，如圖2所示，左側極板接正弦電壓源U=Amsin(2πft)，Am為電壓幅值，頻率f=500 kHz，兩側極板上各覆蓋厚度d1=d2=1 mm的電介質，右側極板接地，盡管實際放電回路中的電阻、電容、電感會對極板電壓波形造成一定的畸變，但本文側重研究放電過程中的物理過程，故對模型作一定的簡化，不考慮外電路的影響，即認為左側極板電壓始終等于電源電壓。介質板邊界條件設置為電子完全吸收，離子的二次發射系數為0.4。兩側的介質板之中的氣隙充有純氬氣。空間網格與時間步長采用均勻劃分方法，空間網格間距Δx=10-6m，時間步長Δt=2×10-13s。

圖2 仿真區域結構示意

2 結果與分析

2.1 Ar放電形貌的時空演化過程

本小節設置氣體壓強為1 atm，背景氣體溫度設定為300 K，施加的正弦電壓峰值為3 kV，兩介質板間隙為1 mm，介質板介電常數為4.5。仿真共持續4.5 μs，即2.25個周期，結合仿真結果分析前兩周期的放電的形貌與參數特征。經實際仿真測試表明，在本文給定的工況下從第一周期負半周開始，各放電物理參量變化規律即可反應放電穩定后的狀態。

圖3展示了放電過程中的極板間電流密度、氣隙電壓與外施電壓的關系，按照放電電流密度曲線特征劃分為A-G七個時刻，其具體時間與對應的放電的物理過程如表2所示，A-G階段，左側介質板與等離子體邊界處的電勢始終大于右側介質板與等離子體邊界處電勢，即左側介質板邊界為陽極，右側介質板邊界為陰極。

t/s

表3 不同標識對應的時間及物理過程

可以看出，每半個周期電流密度曲線有一個明顯的尖峰，伴隨著氣隙電壓的明顯下降。這是因為當氣隙電壓升至一定數值時，產生了強烈的電子崩并引起氣隙擊穿，氣隙中充斥著大量的帶電粒子使得電導率上升，電流密度也隨之上升，氣隙電壓也有明顯的下降。但是氣隙擊穿后電流密度沒有一直維持在一個很高的數值，而是迅速下降，這是因為放電過程中大量電荷會積聚在介質板上形成反向電場，抑制外施電壓形成的正向電場，使得放電熄滅。

在后半個放電周期中，之前介質板表面積累的電荷產生的電場與外施電壓產生的電場同向相疊加，會促進放電的提前發生。由此導致后續放電過程電流密度尖峰的相位提前于電壓峰值約0.63π。放電電流相位變化原因是介質板表面電荷產生的電場對于外施電場在等離子體區域的影響。

從圖4中可以看出介質板表面電荷電場對于放電的影響。在t=0起始時刻，兩側介質板均無電荷，外部電場與外施電場保持一致，隨外施電壓的升高而升高，直至發生擊穿，首次主放電的電流尖峰相位落后于電壓零點，此時大量帶電粒子從氣隙中產生并在電場的作用下向介質板運動并吸附在其表面形成介質板表面電荷，迅速形成反向電場，使得外部總電場迅速減小至零并反向增加，使正半周的放電迅速熄滅，需要說明的是，該工況下介質板表面電荷產生的最大電場的強度約為外施電場強度最大值的3倍，占外部電場的主要部分。從第一個周期的負半周開始，主放電發生于電壓前半周期的下降段，超前于電壓零點相位，這是正向的電場減小的過程中，對介質板表面電荷的反向電場的削弱作用減弱，從而積累出一個更強的反向外部電場，進而發生反向的擊穿過程。

t/s

從圖5的電子密度與圖6的電子能量的時空分布可以看出，介質阻擋放電前期存在一個湯森放電的過程，從A時刻開始，電子崩從陰極向陽極發展；而之后從B時刻到D時刻，電子崩又快速從陽極發展到陰極附近，這個過程中發生強烈的電離作用，電子密度和離子密度的峰值快速增大，如圖5，圖7所示。之后由于介質板表面電荷的積累產生的反向電場的抑制作用，電子崩運動速度逐漸降低。而大質量的陽離子緩慢向陰極運動，在這個階段表現為積聚在陰極附近形成鞘層結構，在鞘層區域內，電子與離子在較強的電場作用下加速，獲得更高的能量。

圖5 電子數密度的時空演化圖

圖6 電子能量的時空演化圖

圖7 離子數密度的時空演化圖

與此同時，由于復合作用，由公式14可知，在電子與離子數密度較高的區域，復合速率較大，大量的電子和離子結合成為中性氣體原子，電子離子的數密度的峰值隨后快速下降。

D時刻到F時刻是一個從湯森放電到輝光放電的轉化過程。此前B-D產生的主電子崩的頭部位置由于強烈的電離作用產生大量的陽離子，向陰極運動一段距離之后撞擊陰極介質板產生二次電子。這些二次電子在D-E時間段內在由鞘層電場加速后，這一部分高能電子向陽極運動，發生一個二次湯森放電的過程。E-F時刻，電子崩向陰極運動，和B-D過程相似，只是強度有所減弱。此后鞘層空間正電荷減少，介質板表面電荷增加，鞘層電場減弱，直到無法再建立起足夠強的電子崩，放電便進入了穩定的輝光放電模式。

從F至G點屬于輝光放電階段，如圖9所示，左側等離子體正柱區的電子與離子密度幾乎一致，而右側的鞘層區域內離子密度大于電子密度。此外，左側正柱區的電場強度相對較低，進入鞘層區域，電場強度快速升高。整個輝光放電階段，鞘層區域的寬度較穩定，維持在0.25～0.27 mm。隨著放電過程的發展，鞘層內的最大電場強度呈逐漸降低的趨勢，這是因為隨著電荷在介質板上的積累造成反向電場的增強，以及復合作用，鞘層內的帶電粒子逐漸減少，空間場強也隨之降低。

2.2 外施電壓對放電特性的影響

在該小節中，保持電源的頻率為500 kHz不變，改變外施電壓依次為2，3，4 kV。定義歸一化電壓為電源電壓與電源電壓峰值的比值，用以判別放電電流脈沖相位。從圖7看，隨著外施電壓的升高，電流密度峰值隨之升高，電流尖峰的相位隨之提前。這是因為外施電壓幅值越高，相同時間內電場變化越快，就能夠在越短的時間內達到擊穿場強引發放電；同時，更高的電壓對應更大的電場強度，能夠給予電場中的帶電粒子更大的能量并將更多的氣體電離。當外施電壓為4 kV時，和外施電壓為2 kV和3 kV相比，每半個周期主放電后的后續放電脈沖的峰值更大，且有明顯的波動性，這是因為電壓較高時帶電粒子的能量較大，能夠在主放電熄滅后產生更加明顯的二次電子湯森放電現象，并帶來一定的放電不穩定性。

t/s

2.3 氣隙間距對放電特性的影響

該小節保持外施電壓頻率與幅值不變，改變氣隙的間距依次為1、2、3 mm。從圖11中可以看出，放電電流密度峰值對應的相位隨著氣隙間距的增大而依次滯后。

t/s

圖12給出了不同氣隙間距下的介質板表面電荷和外施電場情況，介質板表面電荷電場強度和外施電場強度均隨間隙距離的增加而減小。在該小節的工況下，氣隙間距為1和2 mm時電源電壓為負值時候即發生正向擊穿，而間距為3 mm時電源電壓為正值時發生正向擊穿。為達到擊穿場強，越長的間隙越需要減弱外施電場的反向抑制作用乃至增強正向的促進作用，使放電電流的相位更加滯后。

t/s

2.4 介質板介電常數對放電特性的影響

該小節中保持電壓與氣隙間距不變，依次設置介質板的相對介電常數分別為4.5，8，12。從圖9可以看出，放電電流的峰值隨著介質板介電常數的變大而變大，而電流尖峰對應的相位隨介電常數的增大而略有提前。這是因為介質板的介電常數越大，可以積累更多的表面電荷以產生更大的氣隙電場，使放電能夠更早發生，同時，在半周期內釋放更多的介質板表面電荷也使放電電流宏觀上的增大。在放電形貌上，介質板介電常數較大時，后續的放電脈沖更加明顯，這是介質板間場強增大使電子能量增大所引起。

t/s

3 結論

本文建立了大氣壓下氬氣介質阻擋放電的一維PIC/MCC全粒子仿真模型，模型中考慮了介質板表面電荷積累效應與電子-離子復合過程。此外，還使用了粒子權重自適應和并行計算策略加快模型的運行速度，實現了微秒時間尺度與毫米空間尺度的仿真計算。主要結論如下：

(1)對放電的時空演化過程分析表明，DBD的過程前后分為湯森放電和輝光放電兩個部分，其中介質板表面電荷的積累會促進本次放電熄滅和下一次放電的發生，此外，主放電之后還會有幅值較低的二次電子崩的過程。

(2)電源電壓的幅值升高使放電電流密度增大并使放電的相位有所提前，而高能電子的增多會使放電呈現出一定的不穩定性。

(3)氣隙距離對于放電電流密度幅值影響較小，但對相位有較明顯的影響，氣隙距離的增大會使放電相位有所滯后。

(4)介質板介電常數增大，放電電流密度幅值有所增大且放電電流相位有所提前，同時主放電后續的二次電子湯森放電也更加明顯。