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

楊雄 程謀森 王墨戈 李小康

(國防科學技術大學航天科學與工程學院,長沙 410073)

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

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

(國防科學技術大學航天科學與工程學院,長沙 410073)

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

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

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

1 引 言

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

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

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

在本文中,首先對三維螺旋波放電模型作了詳細的介紹,然后將模型計算結果與實驗結果和經典理論進行了對比,最后利用該模型研究了有限空間螺旋波等離子體放電過程中的波干涉現象.

2 計算模型

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

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

在本節中,分別從電磁場方程、電子輸運方程、重粒子守恒方程、電化學反應、邊界條件以及求解處理等方面對模型進行具體介紹.

2.1 電磁場方程

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

其中j為虛數單位,ω為射頻角頻率,A為磁矢勢,σ為電導率,ε0為真空介電常數,μ0為真空磁導率.電導率在等離子體以外的計算域為常值,在磁化等離子體中可表示為張量形式：

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

等離子體電勢按照穩態條件下的泊松方程來描述：

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

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

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

Emf為電場的交變分量,

2.2 電子輸運方程

嚴格的電子輸運過程應該使用相空間(r,v)的動力學玻爾茲曼方程來描述：

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

其中Re為電子生成源項,Γe為電子密度通量矢,

其中E為靜態電場,是等離子體電勢的空間梯度,μe為電子遷移率,De為電子擴散率.假設電子偏離麥克斯韋分布不遠,將(7)式碰撞項線性化后可以導出：

其中〈?〉是指對?在電子分布函數的速度空間內做積分,下同.

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

假設有M個反應凈生成電子,電子生成源項Re與化學反應速率rk的關系為：

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

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

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

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

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

由碰撞導致的能量交換Sen考慮全部化學反應過程,包含彈性和非彈性碰撞：

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

上式中mAr為Ar的原子質量,彈性碰撞過程中交換能量dek的單位以伏特計.對于非彈性碰撞,能量的損失等于具體的化學反應中能量的變化,詳見2.5節.

2.3 重粒子守恒方程

在等離子體中除電子外還存在基態原子、激發態原子、離子等重粒子,假設由N個反應式生成Q種重粒子,那么存在Q-1個相互獨立的粒子守恒方程：

其中wk為k粒子的質量分數,jk為重粒子擴散通量矢,Rk為k粒子生成速率源項.螺旋波放電是一種非熱平衡電離過程,重粒子的溫度通常在放電過程中的變化較小,忽略重粒子的熱擴散效應,則k組分粒子擴散通量矢為

其中括號里第一項為濃度擴散項,第二項為電場遷移項.Dk,f為k組分粒子擴散系數,zk為k組分電荷量,μk,f為k組分遷移率：

重粒子生成源項由化學反應決定：

υk,j為化學反應計量系數,rj為j反應的反應速率,Mk為重粒子摩爾質量.

方程(18)由Q-1個相互獨立的粒子守恒方程構成,加上質量守恒約束關系式,形成對所有重粒子的完整描述：

2.4 邊界條件

為了對全部計算域等離子體形成統一的解析,本文推導了各態粒子的固壁通量邊界條件.電子在等離子體的邊界上通過熱擴散造成巨大損失,當不考慮電子撞擊壁面形成的反射時,壁面處垂直向外的電子通量矢為

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

對于重粒子,其粒子數的變化不是由于熱運動速度導致,而是受到表面化學反應或壁面處的電場控制,因此其壁面處的擴散通量矢可以定義為

其中Mk,rsurf,k,ck,μm,k分別為k組分質量、表面化學反應速率、摩爾濃度以及擴散系數.等式右邊第一項為化學反應作用,第二項為電場遷移作用,當重粒子為中性粒子時,只受到第一項的作用;當重粒子為離子時,還要受到第二項的控制.括號里為方向控制,保證當電場與壁面垂直時,由于遷移的作用粒子流向外流出.

除通量邊界條件外,在本算例中計算了導電邊界條件下的電離,放電室壁面電勢V=0,在端面上則應用了懸浮電勢的邊界條件：

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

2.5 化學反應

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

表1 化學反應關系式Table 1.Collisions and chemical reactions.

2.6 求解處理

2.6.1 阻抗邊界條件

當射頻功率源頻率為13.56 MHz時,電流在銅質天線中的趨膚深度為17.94μm.對于3D模型,過于細化的網格會呈級數增加PDEs的自由度,由于實際計算資源的限制,要劃分出能解析高頻電流的網格通常很困難.在本文模型的計算中,采用阻抗邊界條件定義射頻天線的激勵方式,只計算天線表面的電流和電場而不對天線劃分三維網格,這種處理方式并不影響實際電磁場分布和天線能量損耗,但能極大簡化計算過程.阻抗邊界條件下由電流源激發的電磁場可以寫成如下形式：

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

2.6.2 對數穩定求解

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

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

在空間網格方面,求解域的主體采用四面體網格,由于等離子體參數在鞘層內變化劇烈,在放電室的壁面處采用了六面體邊界層網格,邊界層網格的厚度呈等比數列,比例系數為1.4,共6層,其中第一層的厚度為0.08 mm,小于當地等離子體Debye長度(參照圖5),對空間網格點采用線性函數的有限元格式進行離散.時間離散為向后差分格式,具體采用了COMSOL MultiphysicsTM求解器的自適應時間步長技術,初始時間步長為1×10-14s.

3 計算結果與分析

3.1 實驗算例驗證

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

圖2 等離子體源實驗平均電子密度與計算值對比(a)Toki[14]實驗值 (B0=45 G,p0=50 mTorr,f=67.8 MHz);(b)Chen[15]實驗值 (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).

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

在螺旋波放電中,電子溫度隨著氣體壓力的升高而下降,這是由于隨著氣體壓力的升高,放電室中的中性粒子數增加,電子與中性粒子間因碰撞導致的能量交換更加頻繁,最終導致電子能量下降.圖3顯示了電子溫度與氣體放電壓力的對應關系,仿真計算的結果與實驗值保持相同的變化趨勢,但電子溫度的絕對值略高.

圖3 Toki[14]等離子體源實驗電子溫度與計算值對比(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 等離子體鞘層理論驗證

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

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

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

圖5顯示了離子漂移速度、當地Bohm速度和Debye長度隨半徑的變化曲線.從圖中可以看出,離子漂移速度在離壁面約0.4 mm處達到Bohm速度,約為數個當地德拜長度,與經典鞘層理論相符.

圖6分別顯示了天線電流為11 A時電子溫度和電勢在放電室中達到穩態時的分布.電子溫度分布相對較為均勻,但在靠近天線壁面的區域內較高,對應此處的能量沉積較多,電子被加熱到更高溫度.等離子體電勢在放電室呈現中間高邊緣低的分布,符合由鞘層穩定的等離子體電勢特性.

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

圖5 離子漂移速度和當地Bohm速度徑向分布以及Debye長度徑向分布Fig.5.The radial distribution of ion diffusion velocity and the local Bohm velocity,and the distribution of the local Debye length.

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

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

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

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

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

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

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

為了對低場放電磁場峰值現象的端面干涉機理進行驗證,本文設計了驗證算例：計算天線電流為3.3 A,放電室分別為200,350以及500 mm時電子密度與磁場強度對應關系.假設端面干涉是低場放電磁場峰值現象形成的主要機理,那么當放電室的長度越大時,干涉的現象越不明顯,所導致的磁場峰值應該越小.計算結果如圖10所示,可以看出放電室長度對磁場峰值的影響與假設高度符合.

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

圖11 軸向射頻磁場等高線圖 (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.

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

為了進一步研究電磁波在軸向上的干涉結構,圖12和圖13分別顯示了Ld=200 mm和Ld=350 mm條件下時間平均的磁場能量密度隨軸向位置的變化關系,可以看出當B0=185 G時干涉現象最為嚴重,此時干涉的區域包含了整個放電室軸向范圍,且在干涉區域中出現了數個磁場能量密度的峰值.

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

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

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

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

5 結 論

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

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

附錄A 電子碰撞關系

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

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

其中lnΛ為等離子體庫侖對數,

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

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

在本文的計算中采用文獻[10]的截面數據:電離截面σion=2.8×10-20m2,電子-中性粒子碰撞截面σen=15×10-20,Ar的一級電離能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)[成玉國,程謀森,王墨戈,李小康2014物理學報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)[成玉國2015博士學位論文(長沙：國防科技大學)]

[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

?國家自然科學基金(批準號：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