孫安邦 李晗蔚 許鵬 張冠軍

(西安交通大學(xué),電力設(shè)備電氣絕緣國家重點(diǎn)實(shí)驗(yàn)室,西安 710049)

流注放電低溫等離子體中電子輸運(yùn)系數(shù)計算的蒙特卡羅模型?

孫安邦?李晗蔚 許鵬 張冠軍

(西安交通大學(xué),電力設(shè)備電氣絕緣國家重點(diǎn)實(shí)驗(yàn)室,西安 710049)

流體或者粒子-流體混合數(shù)值仿真是研究流注放電基本物理機(jī)制的常用手段,而精確的電子輸運(yùn)系數(shù)是保證其仿真正確性的必要前提.鑒于現(xiàn)有電子輸運(yùn)系數(shù)求解工具存在一定缺陷,本文開發(fā)了采用蒙特卡羅方法求解低溫等離子體中電子輸運(yùn)系數(shù)的仿真工具,測試表明其準(zhǔn)確性和精確度均較高.研究了氮氧氣體混合比及大氣壓下三體碰撞吸附對電子輸運(yùn)系數(shù)的影響.氮?dú)庵辛髯⒎烹姺抡姹砻?流體仿真中采用本模型改進(jìn)后的電子輸運(yùn)系數(shù)可顯著改善流注通道內(nèi)部的等離子體參數(shù)分布.

流注,低溫等離子體,電子輸運(yùn)系數(shù),蒙特卡羅模型

1 引 言

流注放電(streamer discharge)通常呈現(xiàn)為明亮的絲狀結(jié)構(gòu),是大氣壓下氣體放電初始階段的一種常見現(xiàn)象.由于流注頭部強(qiáng)烈的電荷聚集效應(yīng),使其可以向局部電場低于擊穿電場的區(qū)域傳播［1］.在自然界中,流注在閃電的形成階段起著重要作用,為其產(chǎn)生和傳播提供了預(yù)電離通道［2］.另外,在高空低氣壓大氣層中,流注通常以瞬態(tài)發(fā)光現(xiàn)象(terrestrial luminescence events,TLEs)的形式存在,比如精靈閃電(sprite)等［3］.流注放電在廢氣/廢水凈化、生物醫(yī)學(xué)以及輔助燃燒等方面都有著廣闊的應(yīng)用前景［4］.

自20世紀(jì)40年代流注放電理論被提出以來,科研工作者對流注放電的物理機(jī)制進(jìn)行了大量的數(shù)值仿真工作.基于微觀描述的粒子模型和基于宏觀描述的流體模型是流注放電的兩種主要模擬方法.粒子模型可以較為精確地表述流注放電中的隨機(jī)行為,但是由于需要跟蹤百萬甚至千萬量級的粒子,對計算機(jī)的性能要求較高,耗時較長.因此,以往的流注放電數(shù)值模以計算效率較高的流體模型為主.近年來,隨著計算機(jī)性能的提高和數(shù)值計算方法的發(fā)展,基于空間位置的三維粒子-流體混合模型也逐步建立起來［5,6］,其中,流注內(nèi)部電荷密度高且電場梯度較小的區(qū)域采用流體模型,而其外圍以及頭部電場和電荷密度梯度較高的區(qū)域采用粒子模型進(jìn)行描述,既保證了流注頭部敏感區(qū)域仿真的可靠性,又很大程度上降低了采用純粒子方法的計算消耗,是研究流注放電現(xiàn)象的一種高效的仿真手段.

為了保證流體/混合模型的正確性,電子輸運(yùn)系數(shù)的精確獲取是其中的一個重要因素,包括電子平均能量,漂移、擴(kuò)散系數(shù)以及電離系數(shù)等.電子輸運(yùn)系數(shù)的獲得主要有以下三種手段:1)實(shí)驗(yàn)測量;2)求解Boltzmann方程;3)跟蹤粒子運(yùn)動的蒙特卡羅模擬方法.在已有的研究中,BOLSIG+和MAGBOLTZ是兩種比較有代表性的電子輸運(yùn)系數(shù)求解工具［7,8］,長期以來已被廣泛采用.然而,兩種工具的缺點(diǎn)也較為明顯:BOLSIG+基于求解兩項近似(“two-term”approximation)的Boltzmann方程,其主要缺點(diǎn)為計算精度較低,特別是強(qiáng)電場下電子的能量分布函數(shù)呈現(xiàn)出強(qiáng)烈各向異性的情況下,誤差較大［9］;MAGBOLTZ采用蒙特卡羅方法求解電子輸運(yùn)系數(shù),其采用內(nèi)嵌式的粒子碰撞截面,使得用戶很難驗(yàn)證其碰撞截面的可靠程度,并無法求解其未內(nèi)嵌的氣體組分.此外,上述兩種工具都只能求解電子輸運(yùn)的通量輸運(yùn)系數(shù)( fl ux transport coefficient),而無法求解其體輸運(yùn)系數(shù)(bulk transport coefficient).近期,Rabie和Franck［10］開發(fā)了一套基于MATLAB平臺的蒙特卡羅模型METHES,用于計算低溫等離子體中電子的輸運(yùn)系數(shù).與MAGBOLTZ比較,METHES可以讀取基于‘LXCAT’數(shù)據(jù)庫的電子碰撞截面［11］,且可以同時求解電子的通量和體輸運(yùn)系數(shù).然而,METHES在處理電負(fù)性氣體三體碰撞時存在一定的問題,其一次只能計算某一特定電場的特點(diǎn)使得數(shù)據(jù)的獲取較為繁瑣,且采用商業(yè)軟件Matlab開發(fā),經(jīng)作者測試發(fā)現(xiàn),其僅適用于Matlab 2014及更新的版本,使用具有一定的局限性.

鑒于上述原因,本文基于Fortran 95平臺獨(dú)立開發(fā)了一套用于計算放電低溫等離子體中電子輸運(yùn)系數(shù)的蒙特卡羅模型,為流注放電的精確數(shù)值仿真提供便利.本文開發(fā)的工具可以方便地獲取任意氣體及混合比、任意電場范圍內(nèi)的電子輸運(yùn)系數(shù)和能量分布;此外,模型可以方便地采用不同來源的電子碰撞截面,有利于驗(yàn)證數(shù)據(jù)的正確性,這也是精確模擬流注放電過程的先決條件.

2 蒙特卡羅粒子模擬方法

2.1 電子運(yùn)動與碰撞處理

電子的運(yùn)動遵循牛頓-洛倫茲定律,在不考慮磁場作用下,其運(yùn)動方程為:

式中,me為電子質(zhì)量,q為電子電荷,v和x分別為電子的速度和位置矢量,E為電場.

模型采用 ‘Velocity Verlet’格式［12］來描述電子的運(yùn)動,其計算公式為:

其中,vk和xk分別表示電子運(yùn)動k個時間步長后的速度和位置,?t為時間步長,ak=qEk/me表示電子在k?t時刻電子的加速度大小.

由 于 ‘Velocity Verlet’格 式 比 傳 統(tǒng) 的 “蛙 跳(leapfrog)”格式對時間步長的要求較低［12］,且可較為方便地處理一個時間步長內(nèi)粒子的多次碰撞問題,因此在流注放電［13,14］以及分子動力學(xué)［15］粒子模擬領(lǐng)域被廣泛采用.在本文的模型中,空間電荷對電場的影響較小可以忽略,所以可假設(shè)電子的加速度為恒定值a0,(3)式可簡化為

電子與中性氣體之間的碰撞主要包括彈性碰撞、激發(fā)以及電離碰撞,若存在電負(fù)性氣體,需包含電子與中性氣體之間的吸附作用.在?t時間內(nèi),電子與中性氣體碰撞概率可表示為

式中,n0表示中性氣體數(shù)密度,vinc為電子速度,σT(εinc)代表電子與中性氣體之間碰撞總截面積,代表電子能量.假設(shè)目標(biāo)粒子總共有N種碰撞,則

由于在每個時間步長內(nèi)都要計算所有電子的動能以及碰撞概率,計算量很大.本文采用Birdsall和Vahedi提出的偽碰撞(null collision)方法［16］,主要計算步驟如下:

1)計算粒子的最大碰撞頻率vc,max,即

2)計算在?t時間發(fā)生碰撞的最大概率,

3)計算發(fā)生碰撞的最大次數(shù)Ncoll,

其中Ntotal為模擬粒子總數(shù);

4)選取介于［0,1］之間的一個隨機(jī)數(shù)RF,若

發(fā)生第一種碰撞;

發(fā)生第二種碰撞;不做任何處理.

依次循環(huán)步驟4),直至發(fā)生碰撞的次數(shù)達(dá)到Ncoll,碰撞循環(huán)結(jié)束.對于電子發(fā)生碰撞后的速度分配,主要采用下述方式處理.

1)若為彈性碰撞,由于電子和中性粒子的質(zhì)量和速度差異較大,忽略碰撞后電子的能量損失,且碰撞后電子速度服從各向同性分布［17］;

2)若為激發(fā)碰撞,碰撞后的電子能量為

其中,εex為氣體激發(fā)閥能,εinc,e和εscat,e分別為碰撞前后電子的動能,且碰撞后的電子速度按照彈性碰撞后相同的方法處理.

3)若為電離碰撞,由能量守恒得到

其中,εscat,e為初始電子碰撞后的能量,εse為電離產(chǎn)生的二次電子能量,εionz為氣體電離閥能.對于碰撞后電子的能量分配,一種是采用平均分配的方法,即初始電子和二次電子均分電離后的剩余能量;第二種方法是初始電子碰撞后的能量為

其中,RF為一個介于0到1之間的隨機(jī)數(shù),B為已知函數(shù)(單位:eV).碰撞后的二次電子的能量為

模型對兩種不同的能量分配方式進(jìn)行了測試,其對電子輸運(yùn)系數(shù)基本無影響,在本文的計算中均采用第一種即能量平均分配的方法.兩個電子的速度均采用與彈性碰撞相同的方法處理.

4)若電子與電負(fù)性氣體發(fā)生吸附作用,則在模型中直接刪除該電子的所有信息.

2.2 固定模擬電子數(shù)目方法

由于非守恒碰撞(non-conservative collisions)的存在,比如,在電場較高時,電離碰撞引發(fā)電子數(shù)目呈指數(shù)增長,大大減緩了計算收斂的時間,且很容易使得粒子數(shù)目超過計算機(jī)的能力范圍,使得程序中斷;相反,在電場較低時,電負(fù)性氣體的存在使得電子大量吸附在氣體分子上,電子數(shù)的迅速下降會引起計算誤差變大或者出錯.圖1給出了大氣壓下純氧中電子數(shù)目在不同約化電場(reduced electric fi eld)下的變化曲線,可以看出,由于電負(fù)性氧分子的吸附作用,在電場為20 Td(1 Td=10?21V/m2)時,電子數(shù)目從106減少到零的時間約為8 ns,而在1 Td時僅需2 ns.

為了防止計算中電子個數(shù)快速增加或者衰減,本文采用以下策略來維持電子數(shù)目穩(wěn)定:首先設(shè)定模擬中所期望的電子數(shù)目(Ne),若電子數(shù)目增長超過(Ne+?Ne),隨機(jī)刪除?Ne個電子;若電子數(shù)減少至(Ne??Ne),則新增加?Ne個模擬電子,新增電子隨機(jī)拷貝剩余電子中?Ne個電子的屬性.本文模型可根據(jù)需要靈活的設(shè)定?Ne的數(shù)值,經(jīng)測試發(fā)現(xiàn),?Ne設(shè)定為(0.5–1)Ne與不采用固定電子數(shù)目模擬相對比,對計算結(jié)果的影響基本可以忽略.此外,在本文所采用的電場范圍內(nèi),電離所引起的電子數(shù)目增長基本不會超出計算機(jī)的能力范圍,因此,此電子數(shù)控制策略主要在于解決電負(fù)性氣體的吸附效應(yīng)帶來的電子數(shù)劇減問題,且對輸運(yùn)系數(shù)的計算精度影響較小.需要說明的是,以往研究工作者還采用人工電離的方法,即一個電子被吸附則人為添加一個新的電子,以到達(dá)類似的控制目的.若模型中取?Ne=1,則與文獻(xiàn)［18］中的處理基本一致.

圖1 (網(wǎng)刊彩色)大氣壓下純氧中電子數(shù)目在不同電場下隨時間的變化Fig.1.(color online)Decay of the number of electrons for di ff erent reduced electric fi eld in atmospheric oxygen.

2.3 電子輸運(yùn)系數(shù)計算

2.3.1 電子的體輸運(yùn)和通量輸運(yùn)系數(shù)

體輸運(yùn)速度表征在外部電場作用下所有電子的質(zhì)心位置變化,而通量輸運(yùn)速度則表示電子質(zhì)心位置的變化率［6,18,19］.在實(shí)際統(tǒng)計中,質(zhì)心位置不僅僅由電子的漂移決定,非守恒碰撞(電離碰撞或者吸附)所引起的電子能量變化對其影響也較大.相比較而言,體輸運(yùn)系數(shù)更易在實(shí)驗(yàn)中測得,而通量輸運(yùn)系數(shù)在常壓下特別是高電場的情況下較難測量.兩者的關(guān)系可表示為

其中,vb和vf分別表示電子的體以及通量漂移速度,S表示非守恒碰撞效應(yīng).很顯然,在非守恒碰撞不存在的條件下,vb和vf應(yīng)具有相同的數(shù)值.

在蒙特卡羅模型中,電子的體遷移(μb)和擴(kuò)散(Db)系數(shù)分別表示為:

其中,r和〈r〉分別表示電子實(shí)時位置和平均位置.

類似地,電子的通量遷移及擴(kuò)散系數(shù)為:

在實(shí)際應(yīng)用中,體輸運(yùn)系數(shù)和通量輸運(yùn)系數(shù)的適用范圍存在差異.如,體輸運(yùn)系數(shù)一般可由實(shí)驗(yàn)測量直接獲得,適用于從實(shí)驗(yàn)數(shù)據(jù)反推粒子的碰撞截面;而通量傳輸系數(shù)則已證明在流注放電以及阻性板室中放電行為的流體以及混合模擬中較為適用［6,18］.兩種輸運(yùn)系數(shù)的異同及其適用范圍將在后續(xù)的工作中詳細(xì)探討,本文中不做贅述.

2.3.2 電離系數(shù)

若模型中不采用固定數(shù)量的模擬電子個數(shù),電離系數(shù)(α)可直接由電子個數(shù)的增長率獲得,即

其中,Ne(t2)和Ne(t1)分別表示t2和t1兩個時刻的電子總數(shù).

若采用如2.2節(jié)介紹的固定模擬電子個數(shù)的方法,則通過統(tǒng)計平均電離和吸附頻率的方法獲得電離系數(shù),可表示為

2.3.3 電子平均能量

電子的平均能量可表示為

其中,Ne表示電子數(shù)目.

2.4 計算流程

本文計算電子輸運(yùn)參數(shù)的流程如圖2所示.初始時刻設(shè)定所需計算的電場范圍、氣體種類、溫度以及混合比等;當(dāng)前模型可以直接讀取來源于‘LXCAT’數(shù)據(jù)庫的電子與任意氣體的碰撞截面;電子初始位置設(shè)置為0,初始能量假設(shè)為0.1 eV且速度服從Maxwell分布.

圖2 計算流程圖Fig.2.Flow chart of the calculation.

計算過程中,跟蹤電子在電場中的運(yùn)動以更新其速度和位置信息,進(jìn)而通過蒙特卡羅法處理電子與中性氣體之間的彈性以及非彈性碰撞.由于非守恒碰撞的存在,模型設(shè)定固定數(shù)量的計算粒子以避免電子數(shù)的快速增加或者減少.最后,通過模擬足夠長的時間(電子的能量分布達(dá)到穩(wěn)態(tài))并比較當(dāng)前時刻與上一時刻計算參數(shù)的差異判斷是否達(dá)到穩(wěn)態(tài),若電子輸運(yùn)系數(shù)達(dá)到穩(wěn)態(tài),則輸出具體信息,否則循環(huán)上述步驟.

3 模擬結(jié)果與討論

模擬中所采用的參數(shù)如下:氣體壓強(qiáng)p=105Pa,氣體溫度T=293 K.模擬開始時采用的電子個數(shù)Ne=10000.電場范圍取0.1–500 Td.電子的碰撞截面從LXCAT數(shù)據(jù)庫中獲得［11］.

3.1 模型準(zhǔn)確性驗(yàn)證

為了驗(yàn)證模型的準(zhǔn)確性,本文與BOLSIG+以及METHES比較了計算得到的氮?dú)饧把鯕庵械耐總鬏斚禂?shù).選取氮?dú)夂脱鯕馐蔷哂写硇缘?分別表征普通和電負(fù)性兩種不同氣體.電子與氮?dú)夥肿雍脱鯕夥肿拥呐鲎步孛婢鶃碜訠iagiv8.9［11］.圖3給出了電子漂移系數(shù),平均能量,橫向、縱向擴(kuò)散系數(shù)以及電離系數(shù)隨約化電場的變化情況.對比可知,本模型計算得到的電子通量漂移系數(shù)(圖3(a))、電子平均能量(圖3(b))與BOLSIG+及METHES求解所得具有高度的一致性;模型得到的橫向以及縱向擴(kuò)散系數(shù)(圖3(c)和圖3(d))與METHES也較為一致,但是與BOLSIG+有較大差別,特別是在電場大于100 Td后,差別較為突出.這主要是由于BOLSIG+通過求解兩項近似的Boltzmann方程獲得電子輸運(yùn)系數(shù),其精確度對于分子氣體來說較低,且在高電場下電子速度偏離Maxwell分布下其計算誤差更大［18］.此外,由BOLSOG+得到的電離系數(shù)與本模型以及METHES也有細(xì)微差別,如圖3(e)所示.值得說明的是,在此三類模型中,都忽略或未恰當(dāng)處理空氣中的三體碰撞,其影響將在3.3節(jié)中詳細(xì)討論.

圖3 (網(wǎng)刊彩色)當(dāng)前模型計算得到的電子通量輸運(yùn)系數(shù)與METHES以及BOLSIG+計算結(jié)果對比 (a)遷移系數(shù);(b)平均能量;(c),(d)橫向及軸向擴(kuò)散系數(shù);(e)電離系數(shù)Fig.3.(color online)Comparison of the electron fl ux transport coefficients between current model,METHES and BOLSIG+:(a)Drift coefficient;(b)mean energy;(c)transverse and longitudinal di ff usion coefficient;(e)ionization rate coefficient.

3.2 氮氧氣體混合比對輸運(yùn)系數(shù)的影響

不同氮氧混合比氣體中流注的放電形態(tài)近年來被廣泛研究,其中一個重要特征是在含氧量較高的氮氧混合氣體中流注呈現(xiàn)出較大的流注直徑,而在含氧量較小時則呈現(xiàn)出較多類似于樹枝狀的放電細(xì)絲［20］.

本文計算分析了氮氧混合比對電子的輸運(yùn)系數(shù)的影響,如圖4所示.由圖4(a)可以看出,在低電場區(qū),電子遷移系數(shù)隨含氧量的增加而增大,而在高電場區(qū)含氧量對電子遷移率的影響較小;電子能量在電場小于100 Td時呈現(xiàn)出交錯變化(見圖4(b)),這主要是由于電子的主要碰撞類型由彈性碰撞向非彈性激發(fā)碰撞轉(zhuǎn)變造成的,在電場較高時(＞200 Td),氣體混合比對電子能量的影響基本可忽略;由圖4(c)和圖4(d)可看出,在電場低于10 Td時,氣體混合比例對電子擴(kuò)散系數(shù)的影響較小,且縱向和橫向擴(kuò)散系數(shù)基本一致,而在電場大于10 Td后,電子擴(kuò)散系數(shù)含氧量的增加而增大,并且橫向擴(kuò)散系數(shù)要顯著大于縱向,在流注放電的流體模型中應(yīng)注意區(qū)分;由于氧氣的電離能比氮?dú)庖?因此電離系數(shù)隨著含氧量的增加而增大,見圖4(e).

圖4 (網(wǎng)刊彩色)氮氧氣體混合比對電子通量輸運(yùn)系數(shù)的影響Fig.4.(color online)Variation of the electron fl ux transport coefficients with reduced electric fi elds for various N2-O2mixtures.

3.3 三體碰撞吸附對電子輸運(yùn)系數(shù)的影響

需要指出的是,當(dāng)電負(fù)性氣體存在時,低能電子與其之間的吸附作用是必須要考慮的一個重要因素.以純氧為例,在BOLSIG+以及METHES中,基本忽略了三體碰撞吸附效應(yīng):e?+O2+O2→O?2+O2,其電子的損失主要通過附著解離(dissociative attachment)進(jìn)行,即e?+O2→O?+O.通常情況下,三體碰撞吸附與氣體密度與電場大小有關(guān).對純氧而言,氣體密度增加,三體碰撞吸附增強(qiáng)進(jìn)而電子損失加劇;電場增加,電子能量增大,而氧氣的三體碰撞吸附截面隨電子能量增加而減少,因此三體碰撞減弱,電子吸附損失降低.而附著解離與壓強(qiáng)無關(guān),氧氣中附著解離效應(yīng)僅隨電場的增大而增強(qiáng)［19］.圖5給出了本模型計算得到的大氣壓下純氧中三體碰撞吸附效應(yīng)對電子通量輸運(yùn)系數(shù)的影響.可以看出,在低電場時(E/n0＜10 Td),三體碰撞吸附對電子輸運(yùn)系數(shù)的影響較大、不可忽略.因此,由METHES和BOLSIG+所計算得到的電子通量傳輸系數(shù)僅僅適用于壓強(qiáng)較低或者電場較大,從而附著解離占主導(dǎo)地位的情形,這也是使用METHES和BOLSIG+計算電負(fù)性氣體中電子輸運(yùn)系數(shù)的又一個缺陷.

圖5 大氣壓下純氧中三體碰撞吸附對電子通量輸運(yùn)系數(shù)的影響Fig.5.Variation of the electron fl ux transport coefficients due to the three-body attachment in atmospheric pure O2.

圖6 (網(wǎng)刊彩色)大氣壓下三體碰撞吸附對電子體及通量漂移速度的影響Fig.6.(color online)In fl uence of the three-body attachment on the electron transport coefficients in atmospheric pressure.

另外,若考慮三體碰撞吸附的影響,氧氣中電子的體漂移速度和通量漂移速度差別也較大,見圖6.通常來說,吸附可引起電子能量以及體漂移速度的增大或者減少,稱為吸附加熱(attachment heating)或者吸附冷卻(attachment cooling).如圖6所示,對氧氣而言,在E/n0＜0.7 Td時,電子的體漂移速度要小于通量漂移速度,因?yàn)樵诖朔秶鷥?nèi)三體碰撞吸附主要為頭部能量相對較大的電子,從而降低了整體的電子能量和速度,即吸附冷卻;在0.7 Td＜E/n0＜20 Td范圍內(nèi),電子體漂移速度要大于通量漂移速度,因?yàn)榇藭r發(fā)生吸附的主要為尾部能量較低的電子,電子能量和體漂移速度增加,即吸附加熱效應(yīng)發(fā)生;在20 Td＜E/n0＜200 Td范圍內(nèi),電子體和通量漂移速度基本一致,因?yàn)榇朔秶鷥?nèi)由電離碰撞消耗的頭部高能量電子和三體碰撞吸附消耗的尾部低能量電子基本抵消;而在E/n0＞200 Td時,電子體漂移速度又大于通量漂移速度,這主要是由于高電場下高能量電子大量增加,電離碰撞消耗的頭部高能量電子比較于尾部被吸附的低能量電子較少造成的.需要說明的是,在非電負(fù)性氣體中,由于無吸附效應(yīng),因此在電場小于擊穿電壓時體和通量漂移速度是一致的;然而,當(dāng)電場大于擊穿電場并逐漸升高時,其體漂移速度也逐漸高于通量漂移速度.

3.4 流注放電仿真示例

為了進(jìn)一步說明電子輸運(yùn)系數(shù)精確度對模擬結(jié)果的影響,采用流體模型對大氣壓下氮?dú)庵辛髯⒎烹娺M(jìn)行了一維仿真,并與一維PIC/MCC模擬結(jié)果進(jìn)行比較.所采用的流體和粒子仿真模型見文獻(xiàn)［21］.如圖7所示,流體模型和粒子模型采用相同的初始條件:初始電子數(shù)密度呈高斯分布,最大值為1015m?3;背景電場取6×106V/m,氣體壓強(qiáng)為105Pa,溫度為293 K.

圖7 初始電子數(shù)數(shù)密度分布Fig.7.Initial electron number distribution.

圖8給出了電子數(shù)密度以及電場在不同時刻隨空間位置的變化.流體模型中電子的輸運(yùn)系數(shù)分別由BOLSIG+以及本文模型計算獲得.比較粒子模擬和流體模擬結(jié)果可知:隨著時間的增長,最大電子數(shù)密度由初始時刻的1015m?3增大到約4 ns時的3×1018m?3,放電逐漸由電子崩發(fā)展為負(fù)流注放電,且流注通道內(nèi)部電場逐漸降低;流體模擬和粒子模擬所得到的電子數(shù)密度和電場分布曲線變化趨勢基本相同,在流注通道內(nèi)部,使用本文模型計算得到的電子輸運(yùn)系數(shù)可以獲得和粒子模擬基本一致的電子數(shù)密度及電場分布,而采用BOLSIG+中的電子輸運(yùn)系數(shù)則結(jié)果相差較大,很好地體現(xiàn)了準(zhǔn)確電子輸運(yùn)系數(shù)的重要性;另外可以看出,粒子模擬中流注傳播速度比流體模擬中稍大,且在流注頭部區(qū)域,流體模擬得到的電子數(shù)密度與粒子模擬仍有不小的差異,其原因主要是由于流體模型中電子輸運(yùn)系數(shù)為局部區(qū)域電場的函數(shù),而在流注頭部敏感區(qū)域電場變化較小,電子數(shù)密度以及能量的變化卻較為顯著,已有的研究表明,若在流體方程中耦合電子數(shù)密度梯度的變化,可顯著改善流體仿真的精確度［22,23］.

圖8 (網(wǎng)刊彩色)電子數(shù)密度以及電場強(qiáng)度在不同時刻隨空間位置的變化Fig.8.(color online)Temporal evolution of the electron number density and the electric fi eld as a function of the spatial position.

4 結(jié) 論

本文采用自主開發(fā)的蒙特卡羅模型求解低溫等離子體中的電子輸運(yùn)系數(shù),研究了氮氧氣體混合比以及三體碰撞吸附對電子輸運(yùn)系數(shù)的影響,并在氮?dú)饬髯⒎烹娨痪S仿真中測試電子輸運(yùn)系數(shù)的重要性.主要結(jié)論如下:

1)比較本模型以及現(xiàn)有工具BOLSIG+和METHES所得到的電子輸運(yùn)系數(shù)和其他參數(shù),驗(yàn)證了本模型的正確性和精確性;而基于兩項近似的BOLSIG+則在高電場時誤差較大;

2)本模型可以方便獲取任意氣體混合比中的電子輸運(yùn)數(shù)據(jù),且研究表明,氮氧混合比對電子輸運(yùn)系數(shù)的影響較大,軸向和橫向擴(kuò)散系數(shù)差異顯著,在實(shí)際流體仿真時應(yīng)予以區(qū)分;

3)在氣體壓強(qiáng)較高時,電子與電負(fù)性氣體分子的三體碰撞吸附對電子的體和通量輸運(yùn)系數(shù)影響均較大,而BOLSIG+和METHES中未充分考慮此因素,因此其適用范圍受限;

4)氮?dú)庵辛髯⒎烹姷囊痪S仿真表明,采用本模型改進(jìn)精度后的電子輸運(yùn)系數(shù)可以顯著改善流注通道內(nèi)電子數(shù)密度和電場仿真結(jié)果與粒子模擬結(jié)果的符合度,而在流注頭部與粒子仿真還存在差異,分析其原因是由于本模型中電子輸運(yùn)系數(shù)僅表示為局部電場的函數(shù),未考慮流注頭部的電子數(shù)密度梯度造成的,還有待于后期進(jìn)一步研究.

［1］RaizerYP1991GasDischargePhysics(Berlin:Springer)pp324–343

［2］Ebert U,Sentman D 2008J.Phys.D:Appl.Phys.41 230301

［3］Ebert U,Nijdam S,Li C,Luque A,Briels T,van Veldhuizen E 2010J.Geophys.Res.115 A7

［4］Shao T,Yan P 2015Atmospheric Gas Discharge and Its Plasma Applications(Beijing:Science Press)pp385–611(in Chinese)［邵濤,嚴(yán)萍2015大氣壓氣體放電及其等離子體應(yīng)用 (北京:科學(xué)出版社)第385—611頁］

［5］Li C,Teunissen J,Nool M,Hundsdorfer W,Ebert U 2012Plasma Sources Sci.Technol.21 055019

［6］Li C,Ebert U,Hundsdorfer W 2012J.Comput.Phys.231 1020

［7］www.bolsig.laplace.univ-tlse.fr version:07/2015［2017-6-13］

［8］Biagi S F 1999Nucl.Instrum.Methods A421 234

［9］Bankovic A,Dujko S,White R D,Buckman S J,Petrovic A L 2012Eur.Phys.J.D66 174

［10］Rabie M,Franck C M 2016Comput.Phys.Commun.203 268

［11］http://lxcatnet［2017-6-13］

［12］Sun A B,Becker M M,Lo ff hagen D 2016Comput.Phys.Commun.206 35

［13］Sun A B,Teunissen J,Ebert U 2013Geophys.Res.Lett.40 2417

［14］Sun A B,Teunissen J,Ebert U 2014J.Phys.D:Appl.Phys.47 445205

［15］Lippert R A,Bowers K J,Dror R O,Eastwood M P,Gregersen B A,Klepeis J L,Kolossvary I,Shaw D E 2007J.Chem.Phys.126 046101

［16］Birdsall C K 1991IEEE Trans.Plasma Sci.19 65

［17］Yous fiM,Hennad A,Alkaa A 1994Phys.Rev.E9 4

［18］Bo?njakovi? D,Petrovic Z L,White R D,Dujko S 2014J.Phys.D:Appl.Phys.47 435203

［19］Dujko S,Ebert U,White R D,Petrovic Z L 2011Jpn.J.Appl.Phys.50 08JC01

［20］Wormeester G,Pancheshnyi S,Luque A,Nijdam S,Ebert U 2010J.Phys.D:Appl.Phys.43 505201

［21］Teunissen J,Sun A B,Ebert U 2014J.Phys.D:Appl.Phys.47 365203

［22］Li C,Brok W J M,Ebert U,van der Mullen J J A M 2007J.Appl.Phys.101 123305

［23］Li C,Ebert U,Hundsdorfer W 2010J.Comput.Phys.229 200

Monte Carlo simulations of electron transport coefficients in low temperature streamer discharge plasmas?

Sun An-Bang?Li Han-WeiXu Peng Zhang Guan-Jun
(State Key Laboratory of Electrical Insulation and Power Equipment,Xi’an Jiaotong University,Xi’an 710049,China)

Streamer is usually present at the initial stage of atmospheric pressure air discharge,which occurs in nature as a precursor to lightning,transient luminous events in upper atmosphere and has much potential applications in industry,such as the treatment of polluted gases/liquids,assisted combustion,plasma enhanced deposition etc.Streamer is a multi-scale problem both in time and in space,which brings much difficulty to the conventional diagnostic approaches.In past decades, fl uid or particle- fl uid hybrid models have been frequently used for understanding the mechanisms of streamer discharges because of their high efficiencies of calculations.Accuracies of the electron transport coefficients(including drift/di ff usion coefficient,ionization/attachment coefficient,electron mean energy and extra)play a key role in ensuring the correctness of the fl uid or hybrid simulations.As far as we know,BOLSIG+and MAGBOLTZ are two typical tools for obtaining the electron transport coefficients and have been widely utilized in previous models.BOLSIG+uses“two-term”approximation which is not sufficient for some molecular gases,MAGBOLTZ cannot calculate the bulk transport coefficients:these data are required for some models.METHES is an additional tool for computing electron transport coefficients,however,speci fi c platform is required which is not very user-friendly.As sorts of drawbacks exist in currently available calculating tools,in the paper,a Monte Carlo model is developed for computing the electron transport coefficients in gases,the model is fl exible to choose any type of gas mixture and its accuracy has been validated by comparing with BOLSIG+and METHES.Furthermore,the in fl uences of N2-O2mixture and three-body attachment process in high gas pressures on the transport coefficient are investigated.It is worth mentioning that three-body attachment process can signi fi cantly change the electron transport properties at a relatively low reduced electric fi eld.Thus,speci fi c attention must be paid to the transport coefficients if simulation is performed at a high pressure.In addition,di ff erences between the bulk and fl ux coefficients are analyzed which are not distinguished in some previous models.Finally,we further validate the present Monte Carlo model by performing simulation of streamer discharge in atmospheric N2,which shows that the improved electron transport coefficient from our Monte Carlo model can improve the simulated plasma properties,in particular at the interior of the streamer channel.The existence of divergence at the tip of the streamer channel might be due to our local fi eld approximation;if a density gradient term is included in the impact ionization term and local electron energy approximation of the electron transport coefficients is used,the accuracy of the fl uid can be improved further.

streamer,low temperature plasma,electron transport coefficient,Monte Carlo model

13 June 2017;revised manuscript

17 July 2017)

(2017年6月13日收到;2017年7月17日收到修改稿)

10.7498/aps.66.195101

?國家自然科學(xué)基金(批準(zhǔn)號:51777164)、西安交通大學(xué)“青年拔尖人才支持計劃”(批準(zhǔn)號:DQ1J008)、電力設(shè)備電氣絕緣國家重點(diǎn)實(shí)驗(yàn)室(批準(zhǔn)號:EIPE17311)和中央高?；究蒲袠I(yè)務(wù)費(fèi)專項資金(批準(zhǔn)號:1191329723)資助的課題.

?通信作者.E-mail:anbang.sun@xjtu.edu.cn

?2017中國物理學(xué)會Chinese Physical Society

PACS:51.10.+y,51.50.+v,31.15.–p

10.7498/aps.66.195101

*Project supported by the National Natural Science Foundation of China(Grant No.51777164),the “Young Talent Plan”of Xi’an Jiaotong University,China(Grant No.DQ1J008),State Key Laboratory of Electrical Insulation and Power Equipment,China(Grant No.EIPE17311),and the Fundamental Research Funds for the Central Universities,China(Grant No.1191329723).

?Corresponding author.E-mail:anbang.sun@xjtu.edu.cn