石 峰，王 昊

（河南理工大學 物理與電子信息學院，河南 焦作 454000）

平等板等離子體鞘層演化過程的PIC模擬研究

石 峰，王 昊

（河南理工大學 物理與電子信息學院，河南 焦作 454000）

為了詳細描述等離子體形成過程中鞘層的產(chǎn)生、演化過程，建立了平等板內部等離子體流動過程的PIC模型。在考慮壁面二次電子發(fā)射的前提下，模擬得到了放電開始直到放電穩(wěn)定后的平板等離子體內部和邊界處的電場強度、電子密度、離子密度、電勢等一系列參數(shù)分布，并分析了二次電子發(fā)射系數(shù)對鞘層厚度的影響分析。模擬結果表明：平等板內部電勢相對于邊界處較高，中部的電場強度為零，電子密度近似等于離子密度，而在邊界處電子密度略低于離子密度，電場方向指向邊界處，在等離子體邊界處有明顯的鞘層。

等離子體鞘層；粒子網(wǎng)格方法（PIC方法）；演化過程

0 引言

在有界的等離子體邊緣存在一個電位變化區(qū)域，電子的遷移率遠大于離子的遷移率，因此相對于接地的空間物體，等離子體是帶正電的，這個位于等離子體和物體之間的非電中性區(qū)域，就是鞘層[1]。等離子體鞘層在平面磁控濺射、刻蝕、集成電路、材料的表面改性等領域具有重要的意義和價值[2-3]。

1929年Langmuir提出了等離子體鞘層的概念，Boeuf采用二維流體模型模擬了陰極鞘層區(qū)的電子和離子動力學過程[4-5]。采用PIC/MCC方法研究了氬氣輝光放電陰極鞘層區(qū)域的微觀特性，研究了陰極鞘層區(qū)域內的離子能量分布、角向分布、電子能量分布和電子碰撞散射角分布等微觀特性。趙曉云等[5]研究了伴有二次電子發(fā)射弱等離子體鞘層中離子密度的影響[6]。

郝美蘭[7]用二維流體動力學理論對射頻偏壓等離子體鞘層進行了數(shù)值模擬，結果表明：鞘層電勢隨著低頻源功率的變化而降低，隨著壓力的增加，鞘層變厚，電勢增加[8]。

Toki等[9]采用動力學模型研究了在磁化等離子體鞘層中二次電子發(fā)射，結果表明發(fā)射系數(shù)增加[10]到達器壁的電子能量流隨著增加。

盡管上述文獻已經(jīng)對等離子體鞘層做了一定的研究，但對等離子體鞘層的具體形成過程分析較少。為此，將建立一個包含的等離子體的自洽的PIC（粒子網(wǎng)格）模型分析帶電粒子的運動，揭示等離子體鞘層的形成過程。并分析了二次電子發(fā)射系數(shù)對鞘層厚度的影響分析。

1 PIC方法

1.1 PIC方法的基本思路

等離子體粒子模擬的基本思路，是假設大量的帶電粒子具有初始的位置和速度，對其統(tǒng)計平均等離子體空間中格點處的電荷和電流密度分布，然后用Maxwell方程組求出網(wǎng)格處的電場和磁場。

真空中Maxwell方程組如式（1）：

除了方程式（1）外，還需要媒質的本構關系式才能最終解決場量的求解問題，如式（2）。

式中：ε為媒質的介電常數(shù)；μ為媒質的磁導率；σ為媒質的導電率。

再利用插值或權重法得出粒子所受到的力，通過求解運動方程式（3）推動粒子得到粒子新的位置和速度。以此循環(huán)下去，通過跟蹤大量帶電粒子的運動，得到宏觀等離子體的物質特性和運動方程。

1.2 模擬初始化

模擬開始前，首先需要對粒子的位置和速度初始化，位置在整個放電空間內隨機分布，速度為Maxwell分布，其中位置空間為一維，速度空間為三維。

1.3 電荷密度統(tǒng)計

對于一維情況下，采用粒子權重方法分配電荷，具體為式（4）：

對宏觀粒子統(tǒng)計平均得出網(wǎng)格點上的電荷密度。

1.4 求解電勢和電場

電荷分配到網(wǎng)格點以后每個網(wǎng)格點的電荷密度，將所得的電荷密度帶入泊松方程，通過求解泊松方程，得到每個網(wǎng)格點處的電勢值。

式中：Φ為電位分布；e為電荷；ε0為空氣介電常數(shù)；ni和ne分別為正離子和電子的數(shù)密度，電場由激發(fā)的電場求出。

1.5帶電粒子在電場中的運動

求出電場后根據(jù)插值法求出每個帶電粒子所在位置處電場的大小，對于單個帶電粒子的運動，需要求解其在電磁場作用下的運動方程如式（6）：

式中：v→為速度矢量；M為帶電粒子質量；x→為位置矢量。采用蛙跳格式方法進行求解粒子運動方程。

采用蛙跳格式推進粒子，電場作用下帶電粒子的離散化運動方程如式（7）：

2 仿真結果及分析

模擬放電空間中假設只存在電場，位置為一維，速度空間為三維。模擬的放電空間的長度為z=2.5 cm、z=0 cm，2.5 cm為接地極板，電勢為零。初始時刻電子溫度為5 eV，離子的溫度為0.026 eV。初始時刻宏粒子個數(shù)為105個，每個宏粒子代表104個真實粒子。放電空間劃分為500個網(wǎng)格，電子時間步長取為10-12s，離子時間步長取為10-9s。

圖1 不同網(wǎng)格點處的電勢分布圖Fig.1 The potential distribution of different grid points

圖2為模擬達到穩(wěn)定狀態(tài)時的e-1密度分布?？梢钥闯觯趦蓚€極板處，電子數(shù)密度比較小，形成一個梯度，在中部等離子體區(qū)域，電子數(shù)密度比較高，達到10×1018cm，而在等離子體中部區(qū)域，電子數(shù)密度比較均勻，體現(xiàn)了等離子體的電中性。

圖2 放電空間內的電子數(shù)密度分布圖Fig.2 The electron density distribution inside the discharge space

圖3為電荷數(shù)密度分布?？梢钥闯觯谇蕦舆吔缣?，電荷密度為正值，表明離子數(shù)密度大于電子數(shù)密度，這是由于電子的運動速度快，消失的數(shù)量多，而離子的質量大，速度低，離子消失的速度慢。而鞘層形成以后，粒子數(shù)趨于平衡，不再變化。而在等離子體中部區(qū)域，電子和離子的數(shù)目趨于相等，因此總的電荷密度為零。

圖3 不同網(wǎng)格點處的電荷密度分布圖Fig.3 The charge density distribution of different grid points

圖4 放電空間內的Ar+密度分布圖Fig.4 TheAr+density distribution inside the discharge space

圖5 放電空間內的Ar數(shù)密度分布圖Fig.5 TheArdensity distribution inside the discharge space

圖6表明了極板之間的電場強度E→分布，在左側邊界處，電場強度的值為-100 V/m，其方向指向邊界處，這是由于電子的質量比較小，熱運動速度快，電子打向邊界后消失，而離子的質量比較大，熱運動速度慢，電子的消失速度大于離子，因此邊界處存在過剩的正電荷，因此電場強度E→的方向指向外部邊界。而在右側邊界處，存在過剩的離子，因此電場強度E→的方向指向右側邊界，因此為正值。

圖6 放電空間內的電場強度分布圖Fig.6 The electric field intensity distribution inside the discharge space

圖7表明了不同網(wǎng)格點處有無二次電子發(fā)射的電勢的對比圖，二次電子發(fā)射系數(shù)為0.8。可以看出由于電子達到壁面后消失，電子的損失數(shù)量比較多，導致中部等離子體處的電勢較高，同時鞘層的厚度較小。而有二次電子的情況，由于到達邊界處的電子會發(fā)射回來，導致電子的數(shù)目較多，因此鞘層的厚度比較大，同時由于電子數(shù)較多，導致中部等離子體區(qū)域電勢相對較低。

圖7 二次電子發(fā)射對等離子體鞘層的影響圖Fig.7 the effect of the two electron emission on the plasma sheath

3 結論

通過PIC方法利用一維模型模擬了極板之間等離子體的動力學行為，模擬結果結論為：

（1）由于電子的質量遠小于離子，因此電子的速度遠大于離子，電子在邊界上消失的速度遠大于離子，在邊界處的電子數(shù)量小于離子的數(shù)量，而在中部的等離子體區(qū)域，電子和離子的數(shù)量近似相等，體現(xiàn)了等離子體的電中性性質；

（2）在等離子體邊界處很短時間內形成鞘層，粒子數(shù)在很短時間內達到平衡，邊界處的電場阻礙電子和離子向邊界運動；

（3）粒子的運動區(qū)域平衡以后，在邊界處形成一個電場強度，方向指向邊界，從而阻礙電子向邊界運動，粒子數(shù)達到平衡。而在中部等離子體區(qū)域，電場強度幾乎為零，并且其電勢大于邊界區(qū)域，從而體現(xiàn)了等離子體的電中性性質。

[1]Horisawa H，Sawada F，Hagiwara S，et al.Micro-multiplasmajet array thruster for space propulsion applications[J].Vacuum，2010，85（5）：574-578.

[2]Shiina T，Shima T，Suzuki Y，et al.Inhibitory actions of a local neural reflex on propulsive activity of the esophageal striated muscle portion in rats[J].Research in Veterinary Science，2013，94（2）：331-335.

[3]吳靜，劉國，姚列明，等.等離子體鞘層附近塵埃顆粒特性的數(shù)值模擬[J].物理學報，2012，61（7）：301-307.

[4]Sch?nherr T，Komurasaki K，Romano F，et al.Analysis of At?mosphere-Breathing Electric Propulsion[J].IEEE Transac?tionson Plasma Science，2015，43（1）：287-294.

[5]趙曉云，劉金遠，段萍，等.伴有二次電子發(fā)射的磁化等離子體鞘層結構特性[J].真空科學與技術學報，2012，32（4）：279-284.

[6]Toki K，Kuninaka H，Shimizu Y，et al.Accomplishment and prospect of isas electric propulsion[J].Acta Astronautica，2002，50（5）：285-294.

[7]郝美蘭.射頻偏壓等離子體鞘層特性的二維流體力學模擬[D].大連：大連理工大學，2009.

[8]Shiina T，Shima T，Suzuki Y，et al.Inhibitory actions of a local neural reflex on propulsive activity of the esophageal striated muscle portion in rats[J].Research in Veterinary Science，2013，94（2）：331-335.

[9]Toki K，Shinohara S，Tanikawa T，et al.Small helicon plasma source for electric propulsion[J].Thin Solid Films，2006，507（13）：597-600.

[10]Okawa Y，Kitamura S，Kawamoto S，et al.An experimental study on carbon nanotube cathodes for electrodynamic tether propulsion[J].Acta Astronautica，2007，61（11）：989-994.

PIC SIMULATION OF THE SHEATH EVOLUTION PROCESS IN THE PARALLEL PLATE PLASMA

SHI Feng，WANG Hao
（School of Physics and Electronic Information Engineering，Henan Polytechnic University，Jiaozuo Henan 454000，China）

In order to describe in detail the information and evolution of sheath during plasma formation，the PIC model of internal plasma flow in an equal plate is established.On the premise of considering the two electron emission from the wall，a series of parameters such as electric filed intensity，electron density，ion density and electric potential are obtained at the inner and boundary of the plate plasma at the beginning of the discharge until the discharge is stable.The influence of two electron emission coefficient on sheath thickness was analyzed.Simulation results show that：the internal potential of an equal plate is higher than that at the boundary，the strength of the electric filed in the middle is zero，and the electron density is approximately equal to the ion density.At the boundary，the electron density is slightly lower than the ion density，having an approximate positive electric field distribution，its direction is directed to the border.There is obvious sheath at the plasma boundary.

plasma sheath；PIC method；evolution process

V439

A

1006-7086（2017）06-0345-04

10.3969/j.issn.1006-7086.2017.06.007

2017-07-08

國家自然科學基金（61501175）、河南理工大學博士基金（B2017-56）

石峰（1984-），男，河南人，博士，講師，主要從事電推進技術。E-mail：shf19841009@163.com。