楊雄 程謀森 王墨戈 李小康

(國(guó)防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長(zhǎng)沙 410073)

螺旋波等離子體放電三維直接數(shù)值模擬?

楊雄?程謀森 王墨戈 李小康

(國(guó)防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長(zhǎng)沙 410073)

(2016年5月10日收到;2016年11月1日收到修改稿)

在詳細(xì)考慮電化學(xué)反應(yīng)和粒子碰撞關(guān)系的基礎(chǔ)上,建立了螺旋波放電三維直接數(shù)值計(jì)算模型,舍棄以往模型小擾動(dòng)假設(shè),對(duì)Maxwell方程組直接求解以計(jì)算電磁場(chǎng)能量沉積份額,擴(kuò)展了螺旋波等離子體計(jì)算模型的精度和適用范圍.以Ar為工質(zhì)氣體的仿真結(jié)果顯示：密度躍升效應(yīng)和電子溫度與放電壓力的關(guān)系與Toki等和Chen的實(shí)驗(yàn)結(jié)果符合較好;與經(jīng)典鞘層理論對(duì)比,在粒子數(shù)密度、德拜長(zhǎng)度、電勢(shì)以及電子溫度的分布上取得高度一致,驗(yàn)證了模型的有效性和精度.利用本文模型研究了低場(chǎng)螺旋波放電過程中的磁場(chǎng)峰值現(xiàn)象,驗(yàn)證了放電室端面的波干涉機(jī)理,并發(fā)現(xiàn)波干涉的本質(zhì)是螺旋波分量與其端面回波疊加形成的駐波.

等離子體,螺旋波,直接數(shù)值模擬,波干涉

1 引 言

螺旋波等離子體具有電離效率高、無(wú)電極腐蝕的優(yōu)點(diǎn),在半導(dǎo)體制造以及空間推進(jìn)等領(lǐng)域有著廣泛的應(yīng)用前景.螺旋波電離過程中涉及的復(fù)雜電磁-粒子相互作用并沒有被完全地理解.早期的研究中,朗道阻尼[1]或者螺旋波被認(rèn)為是螺旋波電離的主要能量沉積方式,但隨著研究的深入,越來(lái)越多的研究認(rèn)為TG(Trivelpiece-Gould)波可能才是能量傳遞的主要因素[2-5].

近些年,對(duì)螺旋波等離子體源的放電模擬取得了長(zhǎng)足進(jìn)展,促進(jìn)了對(duì)螺旋波放電機(jī)理的認(rèn)識(shí).Arnush編寫了螺旋波放電一維HELIC代碼,通過電磁場(chǎng)理論解析地計(jì)算螺旋波和TG波能量沉積,但是并沒有考慮等離子體的電離和輸運(yùn)過程;Chen[6]運(yùn)用均勻放電模型,計(jì)算了均勻參數(shù)下螺旋波等離子體源在放電過程中的能量損失特征;Curreli和Chen[7,8]發(fā)展了徑向非均勻螺旋波放電模型,迭代得到等離子體的徑向分布參數(shù),并合理解釋了等離子體徑向分布的原因;Ahedo[9]在考慮等離子體質(zhì)量和動(dòng)量守恒的基礎(chǔ)上提出了徑向密度和速度分布的計(jì)算方法,隨后Ahedo和Jaume[10]將其擴(kuò)充為二維算法,但是其計(jì)算中并未計(jì)算能量沉積過程;國(guó)內(nèi)方面,成玉國(guó)等[11,12]采用與Chen相似的方法計(jì)算了一維電子溫度非均勻的螺旋波放電過程,對(duì)磁場(chǎng)在螺旋波放電過程的作用進(jìn)行了研究.

完善的螺旋波放電仿真應(yīng)包含對(duì)能量沉積和碰撞電離過程的模擬兩部分.對(duì)于能量沉積的模擬,前人的研究幾乎都是在線性小擾動(dòng)的假設(shè)下推導(dǎo)出磁化冷等離子體中螺旋波和TG波的解析解,進(jìn)而方便地計(jì)算電磁波能量沉積.這種方法僅適用于中低射頻能量密度和磁場(chǎng)強(qiáng)度的電離過程,其本質(zhì)上不是完全的直接數(shù)值模擬,計(jì)算精度受到限制.對(duì)于碰撞電離的輸運(yùn)過程模擬,目前的研究最多涉及二維或準(zhǔn)二維模型,且對(duì)于等離子體生成的復(fù)雜電化學(xué)反應(yīng)過程考慮較少,因此計(jì)算結(jié)果與實(shí)際放電狀態(tài)存在差異.本文建立了三維條件下直接數(shù)值模擬的螺旋波放電模型,在充分考慮電離過程中的粒子碰撞和電化學(xué)反應(yīng)因素的基礎(chǔ)上,通過漂移-擴(kuò)散輸運(yùn)方程計(jì)算等離子體各粒子密度和電子溫度,通過實(shí)時(shí)解算的麥克斯韋方程組計(jì)算能量沉積份額,利用COMSOL MultiphysicsTM所提供的高效有限元PDE求解器進(jìn)行全耦合求解.

在本文中,首先對(duì)三維螺旋波放電模型作了詳細(xì)的介紹,然后將模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果和經(jīng)典理論進(jìn)行了對(duì)比,最后利用該模型研究了有限空間螺旋波等離子體放電過程中的波干涉現(xiàn)象.

2 計(jì)算模型

本文以國(guó)防科技大學(xué)束能與電磁推進(jìn)實(shí)驗(yàn)室小型螺旋波等離子體源(SHPS@BEEMP)為研究對(duì)象,建立了三維條件下直接數(shù)值模擬的螺旋波放電計(jì)算模型,如圖1所示.本文計(jì)算了幾種不同射頻天線的電離過程,在第三節(jié)的驗(yàn)證算例中分別重構(gòu)了名古屋III型天線(Nagoya III,m=±1)和三匝環(huán)形天線(m=0),在第四節(jié)波干涉研究中采用Shoji型(雙臂半波長(zhǎng)螺旋形,m=+1)天線,背景磁場(chǎng)和氣體壓力為均勻分布.除非特別說(shuō)明,本文所采用的放電參數(shù)分別為：射頻頻率f=13.56 MHz,放電壓力P0=10 mTorr(1 Torr=133.322 Pa),氣體溫度Tg=300 K,天線長(zhǎng)度La=124 mm,放電室長(zhǎng)度Ld=350 mm,放電室半徑rd=16 mm.

圖1 仿真模型示意圖Fig.1.Schematic diagram of the helicon discharge simulation.

在本節(jié)中,分別從電磁場(chǎng)方程、電子輸運(yùn)方程、重粒子守恒方程、電化學(xué)反應(yīng)、邊界條件以及求解處理等方面對(duì)模型進(jìn)行具體介紹.

2.1 電磁場(chǎng)方程

螺旋波放電過程中的電磁場(chǎng)可以分為兩部分,一部分是以天線為源的交變電磁場(chǎng)(包括誘導(dǎo)電流感生的電磁場(chǎng)),另一部分是由等離子體雙極擴(kuò)散形成的靜電場(chǎng).在本文中,對(duì)于交變電磁場(chǎng)部分在頻域下以磁矢勢(shì)為變量計(jì)算安培定律,對(duì)靜電場(chǎng)部分以電勢(shì)為變量計(jì)算泊松方程.將安培定律在頻域下改寫可得

其中j為虛數(shù)單位,ω為射頻角頻率,A為磁矢勢(shì),σ為電導(dǎo)率,ε0為真空介電常數(shù),μ0為真空磁導(dǎo)率.電導(dǎo)率在等離子體以外的計(jì)算域?yàn)槌Ｖ?在磁化等離子體中可表示為張量形式：

其中電子碰撞頻率νe是電子溫度、電子數(shù)密度以及中性粒子數(shù)密度的函數(shù),具體表達(dá)式參見本文附錄,Bx,By以及Bz為背景磁場(chǎng)的三個(gè)分量,ne為電子數(shù)密度,e為電荷元量.

等離子體電勢(shì)按照穩(wěn)態(tài)條件下的泊松方程來(lái)描述：

其中V為等離子體電勢(shì),ρv為等離子體電荷密度.

方程(1)和(3)描述了等離子體中的電磁場(chǎng)關(guān)系,電磁波在等離子體中傳播存在阻性損耗,表征了放電過程中從電磁波沉積到等離子體中的能量為

其中*為共軛相乘,J為等離子體電流,包括感生電流和位移電流兩部分：

Emf為電場(chǎng)的交變分量,

2.2 電子輸運(yùn)方程

嚴(yán)格的電子輸運(yùn)過程應(yīng)該使用相空間(r,v)的動(dòng)力學(xué)玻爾茲曼方程來(lái)描述：

其中fe為電子相空間分布函數(shù),ue為電子漂移速度,等式右邊為碰撞項(xiàng).方程(7)是含有非線性碰撞項(xiàng)的積分--微分方程,通常難以直接求解.盡管螺旋波等離子體具有較高的電離效率,但其仍然屬于非熱力學(xué)平衡的冷等離子體范疇,大量實(shí)驗(yàn)和理論研究都表明,除非極端的放電條件[13](放電壓力約1 mTorr且同時(shí)射頻功率不低于數(shù)kW/cm2),在絕大多數(shù)情況下螺旋波等離子體中的中性離子數(shù)密度要比電子高約1-2個(gè)數(shù)量級(jí)以上[14-16],因此可以按照弱電離等離子體處理,在此基礎(chǔ)上可對(duì)碰撞項(xiàng)線性化,并將動(dòng)力學(xué)方程在速度空間進(jìn)行積分為流體方程,從而在時(shí)域下的三維空間進(jìn)行高效率求解.忽略放電室中流場(chǎng)的整體運(yùn)動(dòng),電子的粒子數(shù)守恒可以寫成如下形式：

其中Re為電子生成源項(xiàng),Γe為電子密度通量矢,

其中E為靜態(tài)電場(chǎng),是等離子體電勢(shì)的空間梯度,μe為電子遷移率,De為電子擴(kuò)散率.假設(shè)電子偏離麥克斯韋分布不遠(yuǎn),將(7)式碰撞項(xiàng)線性化后可以導(dǎo)出：

其中〈?〉是指對(duì)?在電子分布函數(shù)的速度空間內(nèi)做積分,下同.

其中Te為電子溫度,在本文中以V為計(jì)量單位.

假設(shè)有M個(gè)反應(yīng)凈生成電子,電子生成源項(xiàng)Re與化學(xué)反應(yīng)速率rk的關(guān)系為：

其中NA為阿伏伽德羅常數(shù),cj為反應(yīng)式中j組分摩爾濃度,kf,k為化學(xué)反應(yīng)速率系數(shù),由反應(yīng)的碰撞截面定義：

其中me為電子質(zhì)量,σk(ε)為化學(xué)反應(yīng)碰撞截面,它是電子能量的函數(shù),其離散的數(shù)據(jù)可以通過查詢得到,f(ε)為電子能量分布函數(shù)(EEDF),本文以麥克斯韋分布進(jìn)行計(jì)算.

電子平均能量守恒式可以寫成：

其中nε為平均電子能量密度,Γε為電子能量通量矢,Sen為碰撞造成的能量損失(包括能彈性碰撞和非彈性碰撞),Qdep為電磁波沉積到等離子體中的能量：Qdep=Qrh.電子能量通量矢采用與電子密度通量矢相似的定義：

其中με為電子能量遷移率,Dε為電子能量擴(kuò)散率.當(dāng)EEDF符合麥克斯韋分布時(shí),通過玻爾茲曼方程的積分推導(dǎo)可以證明輸運(yùn)系數(shù)之間滿足相互關(guān)系[17]：με=5μe/3,Dε=5De/3.

由碰撞導(dǎo)致的能量交換Sen考慮全部化學(xué)反應(yīng)過程,包含彈性和非彈性碰撞：

其中F為法拉第常數(shù),dek為能量變化值.電子在彈性碰撞過程中交換的能量為

上式中mAr為Ar的原子質(zhì)量,彈性碰撞過程中交換能量dek的單位以伏特計(jì).對(duì)于非彈性碰撞,能量的損失等于具體的化學(xué)反應(yīng)中能量的變化,詳見2.5節(jié).

2.3 重粒子守恒方程

在等離子體中除電子外還存在基態(tài)原子、激發(fā)態(tài)原子、離子等重粒子,假設(shè)由N個(gè)反應(yīng)式生成Q種重粒子,那么存在Q-1個(gè)相互獨(dú)立的粒子守恒方程：

其中wk為k粒子的質(zhì)量分?jǐn)?shù),jk為重粒子擴(kuò)散通量矢,Rk為k粒子生成速率源項(xiàng).螺旋波放電是一種非熱平衡電離過程,重粒子的溫度通常在放電過程中的變化較小,忽略重粒子的熱擴(kuò)散效應(yīng),則k組分粒子擴(kuò)散通量矢為

其中括號(hào)里第一項(xiàng)為濃度擴(kuò)散項(xiàng),第二項(xiàng)為電場(chǎng)遷移項(xiàng).Dk,f為k組分粒子擴(kuò)散系數(shù),zk為k組分電荷量,μk,f為k組分遷移率：

重粒子生成源項(xiàng)由化學(xué)反應(yīng)決定：

υk,j為化學(xué)反應(yīng)計(jì)量系數(shù),rj為j反應(yīng)的反應(yīng)速率,Mk為重粒子摩爾質(zhì)量.

方程(18)由Q-1個(gè)相互獨(dú)立的粒子守恒方程構(gòu)成,加上質(zhì)量守恒約束關(guān)系式,形成對(duì)所有重粒子的完整描述：

2.4 邊界條件

為了對(duì)全部計(jì)算域等離子體形成統(tǒng)一的解析,本文推導(dǎo)了各態(tài)粒子的固壁通量邊界條件.電子在等離子體的邊界上通過熱擴(kuò)散造成巨大損失,當(dāng)不考慮電子撞擊壁面形成的反射時(shí),壁面處垂直向外的電子通量矢為

其中n為法向元量,vn為電子運(yùn)動(dòng)速度垂直壁面的分量,ve,th為電子熱速度.壁面處損失電子能量通量矢為

對(duì)于重粒子,其粒子數(shù)的變化不是由于熱運(yùn)動(dòng)速度導(dǎo)致,而是受到表面化學(xué)反應(yīng)或壁面處的電場(chǎng)控制,因此其壁面處的擴(kuò)散通量矢可以定義為

其中Mk,rsurf,k,ck,μm,k分別為k組分質(zhì)量、表面化學(xué)反應(yīng)速率、摩爾濃度以及擴(kuò)散系數(shù).等式右邊第一項(xiàng)為化學(xué)反應(yīng)作用,第二項(xiàng)為電場(chǎng)遷移作用,當(dāng)重粒子為中性粒子時(shí),只受到第一項(xiàng)的作用;當(dāng)重粒子為離子時(shí),還要受到第二項(xiàng)的控制.括號(hào)里為方向控制,保證當(dāng)電場(chǎng)與壁面垂直時(shí),由于遷移的作用粒子流向外流出.

除通量邊界條件外,在本算例中計(jì)算了導(dǎo)電邊界條件下的電離,放電室壁面電勢(shì)V=0,在端面上則應(yīng)用了懸浮電勢(shì)的邊界條件：

其中D為位移電場(chǎng),以上條件形成了位移電場(chǎng)不能穿透的絕緣邊界.

2.5 化學(xué)反應(yīng)

螺旋波電離是一種非熱力學(xué)平衡的感應(yīng)電離過程,其中涉及的二級(jí)電離可以忽略不計(jì),在本文算例中考慮了Ar一級(jí)電離過程中所發(fā)生的7組反應(yīng),涉及基態(tài)Ar原子、激發(fā)態(tài)Ar原子(Ars)以及Ar離子(Ar+)三種重粒子,化學(xué)反應(yīng)表達(dá)式如表1所示,其中反應(yīng)5對(duì)于低氣壓的放電維持的等離子體具有最重要的作用.此外,在等離子體的邊界上還考慮了激發(fā)態(tài)Ar原子和Ar離子的復(fù)合反應(yīng)(反應(yīng)I和II).

表1 化學(xué)反應(yīng)關(guān)系式Table 1.Collisions and chemical reactions.

2.6 求解處理

2.6.1 阻抗邊界條件

當(dāng)射頻功率源頻率為13.56 MHz時(shí),電流在銅質(zhì)天線中的趨膚深度為17.94μm.對(duì)于3D模型,過于細(xì)化的網(wǎng)格會(huì)呈級(jí)數(shù)增加PDEs的自由度,由于實(shí)際計(jì)算資源的限制,要?jiǎng)澐殖瞿芙馕龈哳l電流的網(wǎng)格通常很困難.在本文模型的計(jì)算中,采用阻抗邊界條件定義射頻天線的激勵(lì)方式,只計(jì)算天線表面的電流和電場(chǎng)而不對(duì)天線劃分三維網(wǎng)格,這種處理方式并不影響實(shí)際電磁場(chǎng)分布和天線能量損耗,但能極大簡(jiǎn)化計(jì)算過程.阻抗邊界條件下由電流源激發(fā)的電磁場(chǎng)可以寫成如下形式：

其中εr,μr以及σ分別為天線材料的相對(duì)介電常數(shù)、相對(duì)磁導(dǎo)率以及電導(dǎo)率,Ja為天線表面電流密度(注：此處單位為A/m).

2.6.2 對(duì)數(shù)穩(wěn)定求解

由于漂移-擴(kuò)散輸運(yùn)方程本身具有高度的非線性,尤其是在等離子體鞘層中,等離子體的電子密度、電子平均能量密度等參數(shù)具有極大的空間梯度.為了解決這一點(diǎn)給數(shù)值計(jì)算上帶來(lái)的困難,在實(shí)際的計(jì)算中以電子密度、電子平均能量密度和重粒子質(zhì)量分?jǐn)?shù)的對(duì)數(shù)為自變量來(lái)求解輸運(yùn)控制方程,令Ne=lnne,En=lnnε,Wk=lnwk重新改寫輸運(yùn)方程：

(1)式、(3)式,(22)式,(28)式-(30)式包含7組方程以及7個(gè)因變量：磁矢勢(shì)A、電勢(shì)V、電子數(shù)密度對(duì)數(shù)Ne、電子平均能量對(duì)數(shù)En,以及重粒子質(zhì)量分?jǐn)?shù)對(duì)數(shù)WAr,WArs,WAr+,這7組方程以及所述邊界條件構(gòu)成封閉的PDEs,但要對(duì)其進(jìn)行求解還應(yīng)包括合適的初始條件：ne0=1×1014m-3,Te=0 V,A=[0 0 0]Wb/m,V=0 V,wAr=1,wArs=wAr+=0以及初始電中性限制條件ne=ni.

在空間網(wǎng)格方面,求解域的主體采用四面體網(wǎng)格,由于等離子體參數(shù)在鞘層內(nèi)變化劇烈,在放電室的壁面處采用了六面體邊界層網(wǎng)格,邊界層網(wǎng)格的厚度呈等比數(shù)列,比例系數(shù)為1.4,共6層,其中第一層的厚度為0.08 mm,小于當(dāng)?shù)氐入x子體Debye長(zhǎng)度(參照?qǐng)D5),對(duì)空間網(wǎng)格點(diǎn)采用線性函數(shù)的有限元格式進(jìn)行離散.時(shí)間離散為向后差分格式,具體采用了COMSOL MultiphysicsTM求解器的自適應(yīng)時(shí)間步長(zhǎng)技術(shù),初始時(shí)間步長(zhǎng)為1×10-14s.

3 計(jì)算結(jié)果與分析

3.1 實(shí)驗(yàn)算例驗(yàn)證

Toki等[14]開展了小型螺旋波Ar等離子體源實(shí)驗(yàn)探索,放電腔直徑2.5 cm,放電天線為Nagoya III型,利用射頻補(bǔ)償?shù)睦士姞栯p探針研究了不同射頻頻率、射頻功率、氣壓等參數(shù)下等離子體特性的變化趨勢(shì).Chen[15,19]設(shè)計(jì)并實(shí)驗(yàn)測(cè)量了三匝環(huán)形天線螺旋波Ar等離子體源的放電特性,放電室直徑和長(zhǎng)度均為5.0 cm,天線靠近放電室一側(cè)放置.為驗(yàn)證模型的可靠性,本文針對(duì)上述兩種等離子體源分別重構(gòu)了驗(yàn)證算例進(jìn)行計(jì)算,圖2分別顯示了模型的仿真計(jì)算值與兩種等離子體源實(shí)驗(yàn)值隨射頻功率的變化曲線(注：圖2(b)中實(shí)驗(yàn)記錄的功率值為射頻電源輸出值,因此在計(jì)算的過程中考慮了外電路的損耗,外電路電阻值采用Chen建議的1 Ω[15]).

圖2 等離子體源實(shí)驗(yàn)平均電子密度與計(jì)算值對(duì)比(a)Toki[14]實(shí)驗(yàn)值 (B0=45 G,p0=50 mTorr,f=67.8 MHz);(b)Chen[15]實(shí)驗(yàn)值 (B0=80 G,p0=20 mTorr,f=13.56 MHz)Fig.2.The average electron density of simulation vs.experimental results：(a)Experiment by Toki[14](B0=45 G,p0=50 mTorr,f=67.8 MHz);(b)by Chen[15](B0=80 G,p0=20 mTorr,f=13.56 MHz).

在螺旋波電離中存在所謂的“密度躍升效應(yīng)”,即等離子體密度隨著射頻功率的增加存在突然增加的階段,Chen[15]利用螺旋波等離子體放電電阻的非單調(diào)變化對(duì)這一現(xiàn)象進(jìn)行了解釋.圖2顯示了仿真計(jì)算與實(shí)驗(yàn)結(jié)果不僅在數(shù)值上相接近,并且密度躍升區(qū)間也保持了較好的符合,驗(yàn)證了仿真計(jì)算的準(zhǔn)確性.

在螺旋波放電中,電子溫度隨著氣體壓力的升高而下降,這是由于隨著氣體壓力的升高,放電室中的中性粒子數(shù)增加,電子與中性粒子間因碰撞導(dǎo)致的能量交換更加頻繁,最終導(dǎo)致電子能量下降.圖3顯示了電子溫度與氣體放電壓力的對(duì)應(yīng)關(guān)系,仿真計(jì)算的結(jié)果與實(shí)驗(yàn)值保持相同的變化趨勢(shì),但電子溫度的絕對(duì)值略高.

圖3 Toki[14]等離子體源實(shí)驗(yàn)電子溫度與計(jì)算值對(duì)比(B0=45 G,Prf=3-12 W,f=67.8 MHz)Fig.3.The electron temperature of simulation vs.experimental results by Toki[14](B0=45 G,Prf=3-12 W,f=67.8 MHz).

3.2 等離子體鞘層理論驗(yàn)證

有界等離子體邊緣存在等離子體鞘層,具有滯留帶電粒子使等離子體保持穩(wěn)定的作用.鞘層內(nèi)并不滿足準(zhǔn)中性條件,電勢(shì)和電子密度發(fā)生急劇變化,并且這種變化的空間尺度極小,一般僅為數(shù)個(gè)德拜長(zhǎng)度.因此要精確解析鞘層內(nèi)部結(jié)構(gòu)往往存在數(shù)值處理上的困難,在Chen[6]以及Ahedo和Jaume[10]以往螺旋波模型中,采取了一種近似化解決方案：以離子速度達(dá)到Bohm聲速為邊界條件進(jìn)行計(jì)算.這種方法下只解析了等離子體從主體到預(yù)鞘層區(qū)域的特性,對(duì)于鞘層的特性是未知的.在本文的計(jì)算模型中,通過控制實(shí)際固壁邊界的通量條件可以解析真實(shí)鞘層的詳細(xì)結(jié)構(gòu).

根據(jù)等離子體靜電鞘層理論,離子從等離子體主體向壁面進(jìn)行漂移擴(kuò)散,漂移速度隨著電勢(shì)的下降而逐步增加,當(dāng)離子速度滿足Bohm判據(jù)時(shí)達(dá)到預(yù)鞘層,由于質(zhì)量慣性的影響,等離子體開始出現(xiàn)電荷分離,電子密度與離子密度出現(xiàn)顯著差異.從圖4可以看出電子密度與離子密度都隨著與壁面距離減小而下降,其中電子密度在最后的薄層內(nèi)下降尤為明顯.

圖4 鞘層內(nèi)電子與離子密度徑向分布Fig.4.Distribution of electron and ion density in the sheath.

圖5顯示了離子漂移速度、當(dāng)?shù)谺ohm速度和Debye長(zhǎng)度隨半徑的變化曲線.從圖中可以看出,離子漂移速度在離壁面約0.4 mm處達(dá)到Bohm速度,約為數(shù)個(gè)當(dāng)?shù)氐掳蓍L(zhǎng)度,與經(jīng)典鞘層理論相符.

圖6分別顯示了天線電流為11 A時(shí)電子溫度和電勢(shì)在放電室中達(dá)到穩(wěn)態(tài)時(shí)的分布.電子溫度分布相對(duì)較為均勻,但在靠近天線壁面的區(qū)域內(nèi)較高,對(duì)應(yīng)此處的能量沉積較多,電子被加熱到更高溫度.等離子體電勢(shì)在放電室呈現(xiàn)中間高邊緣低的分布,符合由鞘層穩(wěn)定的等離子體電勢(shì)特性.

根據(jù)靜電鞘層理論,Ar等離子體在預(yù)鞘層損失的電勢(shì)為0.5Te,在鞘層勢(shì)降為4.7Te,這兩部分構(gòu)成等離子體的完整勢(shì)降.圖7顯示了等離子體平均電子溫度與電勢(shì)隨半徑的變化曲線,從圖中可以看出,等離子體電勢(shì)隨著半徑的增大而減小,在臨近壁面的鞘層區(qū)域內(nèi)迅速降低;電子溫度隨著半徑的增大出現(xiàn)非單調(diào)的變化,在靠近壁面處取值較大.電子溫度在鞘層內(nèi)約為2.7 V,在預(yù)鞘層內(nèi)約為2.0 V,而等離子體總勢(shì)降為13.65 V,與理論達(dá)到高度的匹配,證明了本文模型計(jì)算的精確度.

圖5 離子漂移速度和當(dāng)?shù)谺ohm速度徑向分布以及Debye長(zhǎng)度徑向分布Fig.5.The radial distribution of ion diffusion velocity and the local Bohm velocity,and the distribution of the local Debye length.

圖6 (網(wǎng)刊彩色)等離子體參數(shù)分布云圖 (a)電子溫度;(b)等離子體電勢(shì)Fig.6.(color online)Contour of plasma parameters：(a)Electron temperature;(b)electric potential.

圖7 電子溫度和等離子體電勢(shì)徑向分布Fig.7.The radial distribution of electron temperature and electric potential.

4 螺旋波等離子體放電的波干涉

Chen等[20]在1989-1992年間的多次實(shí)驗(yàn)中發(fā)現(xiàn)了低場(chǎng)螺旋波放電的磁場(chǎng)峰值現(xiàn)象：當(dāng)射頻放電的能量不太高時(shí),電子密度隨著磁場(chǎng)的增加存在峰值,當(dāng)射頻能量較高磁場(chǎng)密度峰值現(xiàn)象消失,如圖8所示.這一結(jié)果后來(lái)被實(shí)際應(yīng)用于等離子體的發(fā)生：由于不需要維持像電子回旋共振等離子體那樣高的磁場(chǎng),低場(chǎng)螺旋波等離子體能以較低的成本獲得良好的等離子體發(fā)生收益.但是,在隨后相當(dāng)長(zhǎng)的時(shí)間內(nèi)對(duì)于低場(chǎng)放電磁場(chǎng)峰值形成的機(jī)理一直缺乏合理的解釋.2003年,Chen[21]用端面干涉機(jī)理對(duì)這種現(xiàn)象做了解釋,認(rèn)為是電磁波在有限長(zhǎng)的等離子體放電室端面上形成干涉,最終導(dǎo)致磁場(chǎng)峰值現(xiàn)象的形成.

在本文的計(jì)算模型下,同樣觀察到了低場(chǎng)放電的磁場(chǎng)峰值現(xiàn)象,圖9顯示了等離子體密度隨磁場(chǎng)變化關(guān)系.當(dāng)輸入電流不高于5.5 A(射頻功率約100 W)時(shí),發(fā)現(xiàn)峰值磁場(chǎng)隨著射頻輸入能量的增加而升高.當(dāng)射頻輸入電流達(dá)到7.7 A(射頻功率約300 W)時(shí),磁場(chǎng)峰值現(xiàn)象消失,等離子密度隨磁場(chǎng)增加而單調(diào)升高,但在磁場(chǎng)強(qiáng)度約為100 G時(shí)出現(xiàn)了微弱的波動(dòng),這種分布與前面所述因端面干涉所造成的低場(chǎng)放電磁場(chǎng)峰值現(xiàn)象表現(xiàn)出不同的特征,其峰值磁場(chǎng)和峰值等離子體密度都要明顯低于前者.Arnush和Chen[4,22]研究了磁等離子體中的螺旋波和TG波,發(fā)現(xiàn)在磁場(chǎng)很低時(shí)兩者表現(xiàn)出了非常類似的波系結(jié)構(gòu)并容易相互作用,因此此處的弱波動(dòng)很有可能是螺旋波和TG波的共振造成的.

圖8 電子密度隨磁場(chǎng)變化曲線(Rd=4 cm,均勻磁場(chǎng),Chen實(shí)驗(yàn)數(shù)據(jù))Fig.8.Electron density vs.magnetic field strength(Rd=4 cm,uniform magnetic field,by Chen’s experimental results).

圖9 電子密度隨磁場(chǎng)變化曲線(Rd=1.6 cm,均勻磁場(chǎng),仿真結(jié)果)Fig.9.Electron density vs.magnetic field strength(Rd=1.6 cm,uniform magnetic field,by simulation results).

為了對(duì)低場(chǎng)放電磁場(chǎng)峰值現(xiàn)象的端面干涉機(jī)理進(jìn)行驗(yàn)證,本文設(shè)計(jì)了驗(yàn)證算例：計(jì)算天線電流為3.3 A,放電室分別為200,350以及500 mm時(shí)電子密度與磁場(chǎng)強(qiáng)度對(duì)應(yīng)關(guān)系.假設(shè)端面干涉是低場(chǎng)放電磁場(chǎng)峰值現(xiàn)象形成的主要機(jī)理,那么當(dāng)放電室的長(zhǎng)度越大時(shí),干涉的現(xiàn)象越不明顯,所導(dǎo)致的磁場(chǎng)峰值應(yīng)該越小.計(jì)算結(jié)果如圖10所示,可以看出放電室長(zhǎng)度對(duì)磁場(chǎng)峰值的影響與假設(shè)高度符合.

圖10 不同放電室長(zhǎng)度下電子密度隨磁場(chǎng)變化曲線Fig.10.Electron density vs.magnetic field strength at different-length discharge tube.

圖11 軸向射頻磁場(chǎng)等高線圖 (a)B0=300 G左端面;(b)B0=300 G右端面;(c)B0=185 G左端面;(d)B0=185 G右端面Fig.11.Contour line of axial RF magnetic field：(a)Left endplate atB0=300 G;(b)right endplate atB0=300 G;(c)left endplate atB0=185 G;(d)right endplate atB0=185 G.

一維計(jì)算程序[21]無(wú)法對(duì)波干涉的具體結(jié)構(gòu)進(jìn)行解析,利用本文的計(jì)算模型,可以更加直觀地研究螺旋波等離子體放電過程中波干涉的詳細(xì)結(jié)構(gòu).圖11顯示了不同穩(wěn)態(tài)磁場(chǎng)條件下放電室端面附近由天線激發(fā)出來(lái)的軸向射頻磁場(chǎng)等值線圖.對(duì)于開放端面的放電過程,天線激發(fā)的磁場(chǎng)應(yīng)該隨著遠(yuǎn)離天線中心而快速衰減,但從圖中可以看出,在靠近左端面處出現(xiàn)了明顯的磁場(chǎng)峰值,即射頻天線發(fā)出的電磁波在端面處發(fā)射后與回波形成了干涉疊加.由于雙臂半波長(zhǎng)螺旋形天線激發(fā)m=+1模電磁波,在計(jì)算域中向右傳導(dǎo)的能量遠(yuǎn)小于向左傳導(dǎo)的能量,因此在右端面的干涉效應(yīng)也遠(yuǎn)小于左端面.

為了進(jìn)一步研究電磁波在軸向上的干涉結(jié)構(gòu),圖12和圖13分別顯示了Ld=200 mm和Ld=350 mm條件下時(shí)間平均的磁場(chǎng)能量密度隨軸向位置的變化關(guān)系,可以看出當(dāng)B0=185 G時(shí)干涉現(xiàn)象最為嚴(yán)重,此時(shí)干涉的區(qū)域包含了整個(gè)放電室軸向范圍,且在干涉區(qū)域中出現(xiàn)了數(shù)個(gè)磁場(chǎng)能量密度的峰值.

圖12 時(shí)間平均的磁場(chǎng)能量密度沿軸向分布,Ld=200 mmFig.12.The axial distribution of the magnetic energy density averaged by time atLd=200 mm.

圖13 時(shí)間平均的磁場(chǎng)能量密度沿軸向分布,Ld=350 mmFig.13.The axial distribution of the magnetic energy density averaged by time atLd=350 mm.

根據(jù)成玉國(guó)[16]提供的一維無(wú)阻尼條件下TG波和螺旋波傳導(dǎo)公式,計(jì)算了在B0=185 G,電子密度為2×18 m-3時(shí)自由傳播的TG波和螺旋波Bz幅值(如圖14所示),可以看出此時(shí)兩種波表現(xiàn)出不同的傳播特性,螺旋波的波長(zhǎng)遠(yuǎn)大于TG波.特別地,螺旋波的波長(zhǎng)約等于圖12和圖13中磁場(chǎng)能量的峰峰值距離的2倍,證明端面干涉的本質(zhì)實(shí)際上是螺旋波分量與其端面回波疊加形成了駐波.

圖14 一維無(wú)阻尼等離子體中螺旋波和TG波磁場(chǎng)無(wú)量綱幅值比較Fig.14.Dimensionless amplitude of helicon wave and TG wave with one-dimensional undamped condition.

5 結(jié) 論

本文基于碰撞動(dòng)力學(xué)模型簡(jiǎn)化建立了直接數(shù)值模擬的螺旋波三維瞬態(tài)放電模型,利用該模型對(duì)經(jīng)典螺旋波放電實(shí)驗(yàn)進(jìn)行了重構(gòu)計(jì)算,二者的電子溫度和電子密度達(dá)到較為準(zhǔn)確的符合;將仿真結(jié)果與經(jīng)典鞘層理論進(jìn)行了對(duì)比,發(fā)現(xiàn)電子溫度、等離子體電勢(shì)等參數(shù)符合鞘層穩(wěn)定的等離子體特性,證明該模型能夠正確解析等離子體鞘層結(jié)構(gòu).進(jìn)一步地,在本文中利用該模型研究了Shoji型天線低場(chǎng)放電的磁場(chǎng)峰值現(xiàn)象,通過設(shè)計(jì)算例驗(yàn)證了導(dǎo)致該現(xiàn)象的電磁波端面干涉機(jī)理,并發(fā)現(xiàn)波干涉的本質(zhì)是螺旋波分量與其端面回波疊加形成了駐波,這一結(jié)論對(duì)于發(fā)展局部等離子體增強(qiáng)技術(shù)具有指導(dǎo)意義.同時(shí)由于本文計(jì)算模型具有較高的精度,可以用于設(shè)計(jì)改進(jìn)螺旋波等離子體源.

感謝COMSOL公司Dr.Kevin以及王剛博士提供的求解器算法技術(shù)支持.

附錄A 電子碰撞關(guān)系

電子碰撞頻率是復(fù)雜的電子密度、電子溫度函數(shù),對(duì)于氣體電離過程起著決定性的作用.電子的碰撞包含了電子與離子碰撞和電子與中性粒子碰撞兩部分:

電子與離子的庫(kù)侖碰撞頻率可以表示為νei=Keine,庫(kù)侖碰撞速率系數(shù)為

其中l(wèi)nΛ為等離子體庫(kù)侖對(duì)數(shù),

電子與中性粒子碰撞考慮電離碰撞與彈性碰撞兩種情形:

其中電離碰撞和彈性碰撞的速率常數(shù)分別寫成如下形式:

在本文的計(jì)算中采用文獻(xiàn)[10]的截面數(shù)據(jù):電離截面σion=2.8×10-20m2,電子-中性粒子碰撞截面σen=15×10-20,Ar的一級(jí)電離能Eion=15.76 eV.

[1]Chen F F 1991Plasma Phys.Controlled Fusion33 339

[2]Suwon C 1996Phys.Plasmas3 4268

[3]Suwon C 1997Phys.Plasmas4 4167

[4]Arnush D,Chen F F 1998Phys.Plasmas5 1239

[5]Arnush D 2000Phys.Plasmas7 3042

[6]Chen F F 2008IEEE Trans.Plasma Sci.36 2095

[7]Curreli D,Chen F F 2011Phys.Plasmas18 113501

[8]Chen F F,Curreli D 2013Phys.Plasmas20 057102

[9]Ahedo E 2009Phys.Plasmas16 113503

[10]Ahedo E,Jaume N C 2013Phys.Plasmas20 043512

[11]Cheng Y G,Cheng M S,Wang M G,Li X K 2014Acta Phys.Sin.63 035203(in Chinese)[成玉國(guó),程謀森,王墨戈,李小康2014物理學(xué)報(bào)63 035203]

[12]Cheng Y G,Cheng M S,Wang M G,Li X K,Yang X 2014Plasma Sources Sci.Technol.16 1111

[13]Boswell R W 1984Plasma Phys.Controlled Fusion26 1147

[14]Toki K,Shinohara S,Tanikawa T,Shamrai K P 2006Thin Solid Films506-507 597

[15]Chen F F 2007Plasma Sources Sci.Technol.16 593

[16]Cheng Y G 2015Ph.D.Dissertation(Changsha：National University of Defense Technology)(in Chinese)[成玉國(guó)2015博士學(xué)位論文(長(zhǎng)沙：國(guó)防科技大學(xué))]

[17]Hagelaar G J M,Pitchford L C 2005Plasma Sources Sci.Technol.14 722

[18]Dimitris P L,Demetre J E 1995J.Res.Nat.Inst.Stand.Technol.100 473

[19]Chen F F,Torreblanca H 2007Plasma Phys.Controlled Fusion49 81

[20]Chen F F 1992J.Vac.Sci.Technol.A10 1389

[21]Chen F F 2003Phys.Plasmas10 2586

[22]Chen F F,Arnush D 1997Phys.Plasmas4 3411

PACS:52.50.Qt,41.20.Jb,52.70.Ds,02.60.-x DOI:10.7498/aps.66.025201

Three-dimensional direct numerical simulation of helicon discharge?

Yang Xiong?Cheng Mou-Sen Wang Mo-GeLi Xiao-Kang
(College of Aerospace Science and Engineering,National University of Defense Technology,Changsha 410073,China)

10 May 2016;revised manuscript

1 November 2016)

With the detailed consideration of electrochemical reactions and collision relations,a direct numerical simulation model of helicon plasma discharge with three-dimensional fluid-dynamic equations is proposed in the present work.It can improve the precision of results and widen the model applicability by discarding the small perturbation theory in previous helicon models which are partially analytical in essence.Under the assumption of weak ionization,the Maxwell equations coupled with the plasma parameters are directly solved in the whole computational domain.Thus the energy deposited from electromagnetic wave to plasma can be then easily calculated.The values of plasma parameters which include electron density,mean electron energy and heavy species density are obtained by solving a set of drift-diffusion equations.Meanwhile,seven kinds of chemical reactions in the plasma and two kinds of surface reactions on the wall are taken into account.All of the partial differential equations are solved by the finite element solver of COMSOL MultiphysicsTMwith the full coupled method.

The results of numerical cases employing argon as the working medium show that there exists a sharp density jump from a low to high value as the radiofrequency power is raised.The density jump phenomenon is in accordance with the experimental results of Toki(Toki K,Shinohara S,Tanikawa T,Shamrai K P 2006Thin Solid Films506-507 597)and Chen(Chen F F 2007Plasma Sources Sci.Technol.16 593).The electron temperature decreases with an increase of the gas pressure,which is similar to Toki’s(Toki K,Shinohara S,Tanikawa T,Shamrai K P 2006Thin Solid Films506-507 597)measurement by a RF compensation probe.In comparison with the classical sheath theory,the simulation result demonstrates that the distribution of parameters such as particle number density,the Deby length,electric potential and electron temperature can be solved exactly.In addition,the phenomenon of low-field density peak in helicon discharge was studied in the work.Previous research by Chen(Chen F F 2003Phys.Plasmas10 2586)suggests that this peak is caused by constructive interference from the reflected wave.The effect of length of the discharge chamber on the relation of electron density and background magnetic field is investigated numerically.The results validate the mechanism of wave interference reflected by endplates of the discharge chamber.Furthermore,the time-averaged magnetic energy density has more than one peak on the axial direction.Comparing the distribution of the magnetic energy density to that of the dimensionless amplitude of the helicon wave and the TG wave in the one-dimensional undamped condition,it found that the length of peak to peak of the helicon wave is just as twice as that of the magnetic energy density,which indicates that the substance of wave interference is involved in the standing wave generated by the helicon wave and its reflected wave from endplates.

plasma,helicon wave,direct numerical simulation,wave interference

:52.50.Qt,41.20.Jb,52.70.Ds,02.60.-x

10.7498/aps.66.025201

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào)：11305265)資助的課題.

?通信作者.E-mail：yx12321@126.com

*Project supported by the National Natural Science Foundation of China(Grant No.11305265).

?Corresponding author.E-mail：yx12321@126.com