1、Global M在大氣壓射頻放電等離子體數值模擬中的應用東華大學dingke低溫等離子體數值模擬培訓班 SDU 2018第四屆內容一、低溫等離子體基本特性、粒子和能量平衡二、Global M三、 Global M: A Simple Example四、He、Ar大氣壓射頻放電的模擬討論五、He/H2O大氣壓射頻放電的模擬討論六、Ar/O2/TiCl4沉積TiO2的模擬討論七、He/O2/TiCl4射頻放電中的鞘層變化2低溫等離子體基本特性參考資料：1. Chapters 1-9, Lieberman M A and Lichtenberg A J, Principles of Plasma D

2、ischarges and Materials Processing, 2nd edn, (New York: Wiley) 2005;2. Lieberman M A, A Mini-Course on the principles of Low-Pressure Discharges andMaterials Processing, 2010;3Plasma Density Versus Pressure4Non-EquilibriumEnergy coupling between electrons and heavy particles is weakElectrons are not

3、 in thermal equilibrium with ions or neutrals:bulkHigh temperature processing in low temperature plasma For example:A. Electrons produce free radicals - chemistrysurfaceB. Electrons produce electron-ion pairs ion bombardment5ei >> eeTe >> TiDischarge BehaviorUniform density of electrons

4、and ions at time t = 0eiLow mass warm electrons quickly drain to the wall, forming sheath Ion bombarding energy ei = plasma-wall potential VpSeparation into bulk plasma and sheath occurs for ALL discharges6Thermal EquilibriumElectrons generally near thermal equilibriumIons generally not in thermal e

5、quilibriummv2mMaxwellian distribution of electrons= n)exp(-) 2kT3/ 2f (v)(eep kT2eep = nkTPressureFor neutral gas at room temperature, 300 K7n cm-3 » 3.3´1016 pTorrgAverages Over Maxwellian Distribution1、Average energy1 mv22113ò d v 2 mv fe (v) = 2 kTe=32ne2、 Average speed3、 Average e

6、lectron flux lost to a wall4、 Average kinetic energy lost per electron flux lost to a wall8ee = 2Te¥¥¥1eò-¥x ò-¥y ò0zze4eeG=dvdvdv v f (v) =n v m-2s-1v = 1d 3vv2 f (v) =(8kTe )1/2en òep meForces on ParticlesFor a unit volume of electrons (or ions)ruepen m

7、E: electron flow velocitype = nekTe: electron pressure: collision frequency of electrons with neutrals: electric field9electric field forcefriction force (gas drag)m n due = qn r - Ñp - m n nreedte Eeeemuepressure gradient forceBoltzmann Factor for ElectronIf electric field and pressure gradien

8、t forces almost balance-ene E - Ñpe » 0Q E = -ÑF, pe = nekTeÑneÑF = kTeÞenekTe / e = Te (Volts) and integrate to obtainrr10ne (r ) = ne0 exp(F(r ) / Te (V )Plasma Dielectric ConstantRF discharge are driven at a frequency wE(t) = Re(E%e jwt )Define e p from the total cur

9、rent in Maxwells equationsÑ´ H% = J%+ jwe E% º jwe E%C0pJ%= -en u%Conduction currentCeejwmu%e = -eE% - mn mu%eNewtons lawSolve for u%e and evaluate J% to obtainCe2nw= (e )1/ 2pee m0e, eFor w ?nis mainly realpm11w2e= e k= e 1-pep0p0w(w - jn)mPlasma ConductivityIts useful to introduce R

10、F plasma conductivityJ%Cº s p E%E%J%Sinceis a linear function ofCe2ns= epm(n+ jw)mDC plasma conductivity ( w =n m)RF current flowing through the plasma heats electronsJust like a resistor12s= e2n e DCmnm粒子和能量平衡（Argon Discharge）13Electron Collisions With Argone + Ar ® Ar + + 2ee + Ar ®

11、 Ar* + e ® e + Ar + photone + Ar ® Ar + e14Electron Collisions With ArgonMaxwellian electrons collide with Ar atoms (density ng )# CollisionType =n n= Kn ns × m3n : collision frequency s-1 egeK : rate coefficient m3/ sRate coefficient K (T ) is average of cross section s (v)m2 overReM

12、axwellian distributionK (Te ) = s vRMaxwellianvR : relative velocity of collision particles15Ion Collisions With ArgonArgon ions collide with Ar atomsAr + + Ar ® Ar + + ArAr + + Ar ® Ar + Ar +(elastic scattering)(charge transfer)Total cross section for room temperature ions s i» 10-14

13、 cm2Ions-neutral mean free pathPractical formulap(Torr)16l (cm) =1i300 pl » 1in sgiThree Energy Loss ProcessesCollisional energy ec lost per electron-ion pair createdKizec = Kizeiz + KexeexÞ ec (Te )+ Kel (2m / M )(3Te / 2)Te (Volts)Electron kinetic energy lost to wallsee = 2TeIon kinetic

14、energy lost to walls is mainly due to the DC potential V%sacross the sheathei » V%sTotal energy lost per electron-pair lost to walls17eT = ec + ee + eiCollisional Energy Loss18Discharge Behavior: SheathUniform density of electrons and ions at time t = 0eiLow mass warm electrons quickly drain to

15、 the wall, forming sheath Ion bombarding energy ei = plasma-wall potential VpSeparation into bulk plasma and sheath occurs for ALL discharges 19Bohm (Ion Loss) VelocityDue to formation of a presheath, ions arrive at the plasma-sheath edgewith directed energy kTe/ 21 Mu2= kTei22Electron-ion pairs are

16、 lost at the Bohm velocity at the plasma-sheathedge (density ns)20u = u= ( kTe )1/2iBMSteady State DiffusionParticle balance: losses to walls = creation in volumeÑ ×Ge,i =n e,ineAmbipolar: equal fluxes (and densities) of electrons and ionsG = -DaÑnDa = kTe / Mn iBoundary condition: Gw

17、all= nsuBat plasma-sheath edgeÞedge-to-center density ratio21hl º ns / n0Plasma Diffusion at Low PressurePlasma density profile is relatively flat in the centerand sharply falls near the sheath edgeIon and electron loss flux to the walls:Gwall= nsuB = hl n0uBThe edge-to-center density rati

18、o :i = ion-neutral mean free path P 16Applies for low pressure in Argon< 100 mTorr22h º ns»0.86ln(3 + l / 2l )1/20iPlasma Diffusion inLow Pressure Cylindrical DischargeLoss fluxes to the axial and radial walls:The edge-to-center density ratio :Applies for low pressure in Argon< 100 m

19、Torr23h»0.80R(4 + R / l )1/2ih »0.86l(3 + l / 2l )1/2iGaxial = hl n0uB , Gradial = hRn0uBPlasma Diffusion at High PressurePlasma density profile is relatively flat in the centerand sharply falls near the sheath edge ?Ion and electron loss flux to the walls:Gwall= nsuBThe edge-to-center den

20、sity ratio :( Kawamura, PSST 23, 035014, 2014 )24h º ns= 1+ ( l uB )-1/2lnp D0a在一個半徑為R、長度為 l 的圓柱形等離子體中，面積上的粒子損失率是由流向器壁的粒子通量密度決定的。這些nsuB的乘積表示。鞘層邊界的離子密度 ns 與等離通量密度可用子體中心離子密度 n0 之間的關系非常復雜：對不同的氣壓和R、lil 值，離子和電子的雙極性輸運尺度可小至可以將等離子體分成三種狀態：、大到R、l。1）低氣壓狀態，li ³ (R, l) ，即Langmuir狀態。離子的輸運是無碰撞型的，主等離子體區的密

21、度分布可用離子自由落體模型 得到，這個分布曲線在靠近等離子體中心處相對平坦，而在靠近鞘層邊界處則迅速下降。» 0.5 ；無窮長圓柱形體形狀時，25R ? l ，ns/ n0在平板形狀時，當/ n0 » 0.4 。當 l ? R ，ns2）中等氣壓狀態，(R, l) ³ li³ (Ti/ Te )(R, l) 。在這種狀態下，離子輸運由擴散過雙極性電場來決定。但因為在主等離子體區的大部分區域里，離子的漂移速度遠高于熱運動速度，這種條件下的擴散方程是非線性的。可用平板模型得到密度分布： 在中心處也較平坦，在靠近鞘層邊界處陡峭起來。將無碰撞（低氣壓）狀態和碰撞

22、（中等氣壓）狀態下得到的結果擬合，可以得到近似表在軸向鞘層邊界：º nsl0.86»hl(3 + l / 2l )1/ 2n0iº nsR在徑向鞘層邊界：0.80»h(4 + R / l )R1/ 2n0i263）高氣壓狀態，li£ (Ti / Te )(R, l) 。此時，擴散過雙極性電場決定了離子輸運。離子密度的徑向分布可用Bessel函數J0描述，軸向分布可用余弦函數描述。º nsl= 1+ ( luB)-1/ 2hlp Dn0a因為離子碰撞頻率高，不能再假設離子密度在主等離子體區域均勻分布而僅在鞘層迅速下降，通常可假設鞘層邊界

23、處離子密; 0度 nsl，而離子流向器壁的通量密度由其密度梯度給出，對于平板模型，Gi參考資料：Chapters 10= -Dadni / dx 。Lieberman M A and Lichtenberg A J, Principles of Plasma Discharges and MaterialsProcessing, 2nd edn, (New York: Wiley) 200527Global M參考資料：Chapters 10Lieberman M A and Lichtenberg A J, Principles of Plasma Discharges and Materi

24、als Processing, 2nd edn, (New York: Wiley) 200528Global M的引入在低氣壓等離子體中，電源能量主要耦合給了電子。在主等離子體區，離子和中性粒子可以通過很弱的碰撞過程從電子 那里獲得少許能量，此外還可以從較弱的雙極性電場中獲得一 些能量。一個電子通過彈性碰撞轉移給重離子或中性粒子的能量比例為。因此，主等離子體區中的電子溫度通常要遠高于離子溫度和中性粒子的溫度。不過，分解和激發也能夠產生少部分較高能量的重粒子。此外，雙極性電場能夠向鞘層邊界獲得的動能為離子，典型的主等離子子體區的離子到達鞘層時，。29在低氣壓下，能產生電離作用的電子的能量一般在

25、10-15V，其平均自由放電區差不多。因此，若沒有外加磁場，即使電源能量只注入到一個小體積內，各處的電子-中性粒子碰撞頻率仍然具有均勻的空間分布。在高氣壓下，電離電子的平均自由程通常小于放電區因此，當所加電功率的空間分布不均勻時，放電區各處的電離 頻率也不均勻。等離子體的電子分布函數并不一定具有Maxwell分布。不過，可以將分布函數近似于Maxwell分布，并假設主等離子體區有均勻分布的電子溫度和各種電子碰撞反應速率分布，可由 此估計放電參數。30。要準確地求得電子分布函數，需要了解在一定氣壓下的各種碰撞過程以及電子加熱機制。電子和中性粒子之間的碰撞過 程，對粒子的產生，以及一些能量損失過程

26、都很重要。離子和 中性粒子之間的碰撞過程，對于粒子的產生，輸運過程以及粒 子能量分布很重要。氣體放電中含有多種碰撞過程，所以難以確定等離子體中各種粒子數平衡和能量平衡的基本規律。因此，研究基本 的放電過程可以了解等離子體的性質，但對理解氣相化學和表 面化學卻幫助很少。而對大多數材料處理過程，這些化學過程恰恰非常關鍵。材料處理中用到的氣體大都是（并可能是電負性的），其等離子體性質和惰性氣體等離子體差異很大。31因此，可以考慮用Global M來研究氣體放電以及相關的化學反應、表面反應。在該模型中，等離子體的空間變化都是假設的，不是計算得到的。在Global M中，最簡單的是忽略所有的空間變化，稱

27、為零維模型：對復雜放電過程中各種等離子體參數和其變化關系的初步估計非常有用。在實際應用的氣體放電中，會有多種反應和多種正離子，零維模型對研究這種復雜放電很有用。Global M到具體放電參數。通過求解粒子平衡方能量平衡方程來得32Particles Balance and Electron TemperatureAssume uniform cylindrical plasma absorbing power PabsParticle balance :production due to ionization = loss to wallsSolve to obtain= 1Rldeff2 R

28、h + lheffective plasma sizelRGiven ngÞand deffelectron temperature TeTevaries over a narrow range of 2-5 volts33Kiz (Te ) =1uB (Te )ng deffK n n p R2l = (2p R2h n + 2p Rlh n )uizg0l0R0B34Ion Energy For Low Voltage Sheaths= energy entering sheath + energy gained traversing sheatheiIon energy ent

29、ering sheath = Te / 2 (Volts)Sheath voltage determined from particle conservationG= 1 n v e-VG= n u= (8eT (V ) / p m)/T (V )1/ 2vseisBeseee4The ion and electron fluxes at the wall must balance» 4.7T (V ) for ArgonFor example: Vseei » 5.2Te (V )Accounting for the initial ion energy35V = Te

30、(V ) ln( M)s22p mIon Energy For High Voltage SheathsLarge ion bombarding energies can be gained near RF-driven electrodesThe sheath width s is given by the Child lawEstimating ion energy is not simple as it depends on the type ofdischarge and the application of RF or DC bias voltages3642eV 3/2J = en

31、 u=e ()1/2 s isB90Ms2Power Balance and n0Assume low voltage sheaths at all surfaceseT (Te ) = ec (Te ) + 2TeÞ+ 5.2TePower balancepower in = power outSolve to obtaineffective area for particle lossDensity n0 is proportional to the absorbed power Pabspthrough hl, hR, andTeDensity n0 depends on pr

32、essure37n =Pabs=Pabs0(h 2p R2 + h n 2p Rl)u eeAu eelR0BTeffBTP= (h n 2p R2 + h n 2p Rl)u eeabsl0R0BTeT = ec + ee + eiSummaryParticle balance Þ electron temperatureindependent of plasma densityTeÞPower balanceplasma density n0Teonce electron temperatureis known38n =Pabs=Pabs0(h 2p R2 + h n

33、2p Rl)u eeAu eelR0BTeffBTKiz (Te ) =1uB (Te )ng deffExampleLet R = 0.15m , l = 0.3m , ng = 3.3´10m( p = 1mTorr ) at 298 K,19-3= 800WPabsandAssume low voltage sheaths at all surfaces1. Find li= 0.03m, then hl » 0.31, hR » 0.27 and» 0.18mdeffP 16,23,33P 34» 3.8Vng deff2. From

34、the Te versusfigure, Te3. From the ecec» 47V» 5.2TeP 18versus Tefigure,ee = 2Te » 7.6V, and eieT» 20V= 74VP 174. Addingyields» 3.0 ´103 m / s , and findA» 0.12m25. Find uP 20,37P 37Beffn » 1.9 ´1017 m-36. Power balance yields0J= eh n u» 2.8mA / cm27.

35、 Ion current densityill0Bei = Vs » 18VP 358. Ion bombarding energyExample 1, P335, Chapters 10, Lieberman M A and Lichtenberg A J, Principles ofPlasma Discharges and Materials Processing, 2nd edn, (New York: Wiley) 200539Reactive Neutral BalanceAdsorption and Desorption Kineticsin Discharge40Pl

36、ane-Parallel DischargeP 33P 2241Plane-Parallel DischargeP 37P 22P 22P 3542Reactive Neutral Balance43Loading EffectP 4344Adsorption45Sticking Coefficient46Desorption47Adsorption-Desorption KineticsP 46P 4748Adsorption-Desorption Kinetics49此部分略去：請具體查閱參考資料相應章節Chemical Reactions and Equilibrium Molecula

37、r Collisions（Reaction Rates）Chemical Kinetics and Surface Processes參考資料：Chapters 7-9Lieberman M A and Lichtenberg A J, Principles of Plasma Discharges and MaterialsProcessing, 2nd edn, (New York: Wiley) 200550Global MSummary對于大氣壓等離子體（如介質阻擋放電、等離子體射流等），是一種容性耦合等離子體。由于這時放電氣壓較高，對電子和離子成份，均可以采用 遷移-擴散近似。這樣，

38、等離子體的流體力學方程組為：¶ne+ Ñ× =nneione¶t¶ni=n+ Ñ× niione¶t¶¶ÑT ö + × u+ W÷eeeeè2 mèøeene1= -n E -Ñ(nT )em nem ne ee ene ene1 =n E -Ñ(nT )imnimni ii ini inÑ× E = e n (r, t) - n (r, t)eie051P51-56：年，低溫等離子體

39、物理基礎之整體模型部分可以看出，流體方程是在時間和空間上的一組偏微分方程組，需要借助于復雜的數值計算才能完成。如果對該方程進行空間化Global M（1）粒子數平衡方程，就可以得到一種簡為了簡單起見，首先考慮位于兩個無限大平板之間的電正性等離子體（一維幾何模型），兩個無限大平板分別位于一維坐標 x=-l/2處。和x=l/2平面壁平面壁兩個間距為l 的無限大平行板之間的體電離等離子體的一個區域，圖中灰色的強度表示等離子體密度，灰色箭頭表示粒子流的大小和方向。粒子最終到達平板并在平板表面復合。52將連續性方程中的每一項在位置空間中，可得到：¶¶¶(neue )dx =

40、l /2 l /2 l /2 òòòn dx +n n dxe¶xe lt-l /2 -l /2 -l /2 式中第二項可以分成兩部分，并化簡為指向兩個平板的粒子流。由于所用一維幾何模型的對稱性，中心處的粒子流為零：¶(neue )dx = 2l /2 l /2 dG = 2Gòòwall¶x-l /2 0所以，可以得出有效的全局（普適）粒子平衡方程：- 2Gwalldne= nenldtl可以將上面的粒子平衡方程推廣到體積為V 和總面積為A的三維腔室情形。具體方法是：用V/A 取代平板半間隔 l/2，另平板的面積

41、Asect>>l2，并定義積為lAsect 的塊狀等離子體，其面積近似地為2Asect。一53對于穩態等離子體，有æ A ö= GwallK (T ) =n (T ) / nç n V ÷lelegnèøeg上式中的反應系數 Kl(Te) 只是電子溫度Te 的指數函數。由于器壁上的粒子流量與平均電子密度成正比，所以右邊與電子密度無關。氣體數密度ng 與氣壓成正比，所以電子溫度與氣壓及系統尺度有關。（2）能量平衡方程通常在電能維持的放電中，電場能量總是耦合給電子，因此在確定放電參數時，沒有必要考慮離子對電場能量的吸收。對于電

42、場維持的穩態等離子體，電子吸收的能量Pabs 等于損失的能量平衡方程。Ploss，由此以得到整體模型下的能量54電子的能量損失途徑有兩類：（）與中性粒子碰撞，以非彈性碰撞的方式將電子能量轉化為中性粒子的電離能和激發能，或以彈性碰撞的 方式轉化為氣體熱能，（）將動能帶入等離子體界面區。對于惰性氣體 等離子體，電子能量的碰撞損耗可以表示如下：éT ù+ 3mee+ Ke= n ngPêëKiz izKel e úûeloss,collexc excM式中 eizeexc是Wm-3電離能和激發能的是焦耳，Ploss,coll的。在等離子體中

43、，電子還很多其他能量損失途徑，例如工作氣體的解離，振動激發等。在這種情況下，上式式應該包括這些能量損失機制。第二類電子能量損失途徑是電子攜帶能量進入等離子體界面區，電子 能量在界面區靜電場中或在器壁上損失。電子在界面區損失的能量可以表 示如下：55AV= (2kT + eDf)GPloss,bound .ewall式中A為界面表面積，eDf是電子克服界面鞘層勢壘所做的功。因此，在等離子體中體積電子的能量損失為兩種途徑之和：Ploss = Ploss,coll+ Ploss,boundset initial density nedne/dt = 0, obtain Teadvance time

44、stepobtain the other densityYes Nosteady stateoutput data56Global M：A Simple Example57Electron Collisions With Argone + Ar ® Ar + + 2ee + Ar ® Ar* + e ® e + Ar + photone + Ar ® Ar + e58A Simple Example Code59d æ 3 nT ö = Pabs- nn( K (T )e+ K (T )e+ K(T )e)dt ç 2e &

45、#247;VNizeizeleelexeexèø- nUB (Te ) æ eV + 5 T ö - a 3 T n2-dçs2e ÷2eeffèødn = nn K (T ) - nUB (Te ) - a n2dtNizedeffKelKexKizdn =+ dt × f (n,Te , t )f (n,T , t )nnew = noldedtæ 3 nT ö = g (n,T , t )æ 3 nèöd3è 2+ dt × g

46、(n,T , t )dt ç 2e ÷÷eeeèøøoldnew6061pi = 4.0*atan(1.0)N = 10000t = 0.0dt = 1e-7 Pn = 5.0Tg = 300.0ng = 3.25e19*Pn*297/Tg lamda = 1e18/ngalpha = 2e-16 R = 0.07L = 0.02hR = 0.80/sqrt(4.0+R/lamda)hL = 0.86/sqrt(3.0+0.5*L/lamda) ideff = 2.0*(R*hL+L*hR)/R*L Pabs0 = 50Volume

47、= pi*R*R*L e = 1.602e-19Eiz = 15.6Eex = 11.6Eel = 4.085e-5ne = 1e16 te = 0.162function Dte(ne,te)Dte = rPabs(t)/Volume - ne*ng*e*&(rKiz(te)*Eiz+rKex(te)*Eex+rKel(te)*Eel) &- ne*rUb(te)*ideff*(7.2*e*te)&- alpha*3.0/2.0*e*te*ne*ne returnendfunction Dn(ne,te)Dn = ne*ng*rKiz(te) - ne*rUb(te)

48、*ideff - alpha*ne*ne returnendsubroutine ode(ne,te,dt)real ne,te,dtreal new_ne,new_te,ryp,new_yp real Pabs0,Volume,e,Eiz,Eex,Eelcommon /input2/ Pabs0,Volume,e,Eiz,Eex,Eelryp = 3.0/2.0*e*te*nenew_ne = ne + dt*Dn(ne,te) new_yp = ryp + dt*Dte(ne,te) new_te = new_yp/3.0*2.0/new_ne/ene = ne + 0.5*dt*(Dn(

49、new_ne,new_te) + Dn(ne,te) ryp = ryp + 0.5*dt*(Dte(new_ne,new_te) + Dte(ne,te) te = ryp/3.0*2.0/ne/eenddo 220 i=1,Nwrite(2,*), t,ne,te print *,'p = ',rPabs(t) call ode(ne,te,dt)t = t+dtprint *,'t = ',t,'ne =',ne,'te =',te220continue63function rKex(te)real terKex = 0.3

50、71e-13*exp(-15.06/te)+ (0.06271e-13*exp(-14.27/te)*14.3 +0.0352e-13*exp(-15.07/te)*14.15+0.009237e-13*exp(-15.66/te)*14.1 +0.002501e-13*exp(-15.92/te)*14.3)/11.6return endfunction rPabs(t)real t,Tpreal Pabs0,Volume,e,Eiz,Eex,Eelcommon /input2/ Pabs0,Volume,e,Eiz,Eex,Eel if(t .lt. 1e-4) thenrPabs = P

51、abs0 elserPabs = 0 endifC Tp = 1e-2C pi = 4.0*atan(1.0)C rPabs = 2.0*Pabs0*sin(2.0*pi/Tp*t)*sin(2.0*pi/Tp*t)C Tp = 100e-6C ifl = t/TpC if( t-ifl*Tp .gt. 0.2*Tp) then C rPabs = 0.0C elseC rPabs = Pabs0/0.2 C endifC returnendfunction rKiz(te)real terKiz = 2.3e-14*te*0.68*exp(-15.76/te) returnendfuncti

52、on rKel(te)real te real tttt = te*teif (te .gt. 4) then rKel=1e-13*(0.399*te-0.0186*tt+0.004672*tt*te-6.429e-6*tt*tt) elserKel=1e-13*(5.14e-3+5.51e-2*te+0.229*tt-.0642*tt*te+.006*tt*tt) endifendDensity Temperature3.00E+01862.50E+0182.00E+01841.50E+01821.00E+0185.00E+01700.00E+0000.00000.00020.00040.

53、00060.00080.0010time secSimulation Results64Temperature eVDensitySummaryReaction Rate Coefficient I:H2 / Ar / He 65e-M Reaction Rate CoefficientsThe Boltzmann solver Bolsig+ with the Morgan cross section sets is used to calculate the eM (atom or molecule) rate coefficients, and theAT Be-C /Teresults are fit to the formeBoltzmann solver Bolsig+energetique/projets-en-cours/bolsig-resolution-de-l-equation-de/?lang=e