global在大气压射频放电等离子体数值模拟中的应用

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