PIC-MCC程序手冊X.doc

PIC-MCC 程序手冊1. PIC_MCC的模擬方法和數(shù)值計算1.1 PIC_MCC的模擬原理1.2 PIC模型 1.3 MCC模型1.4 模擬中所涉及的放電粒子1.5 模擬中的碰撞2. 數(shù)據(jù)分析2.1程序正常運行所需文件及文件意義2.2主要輸出文件的意義1 PIC_MCC的模擬方法和數(shù)值計算1.1PIC_MCC的模擬原理 目前在計算機模擬中大量采用低壓射頻放電模型來模擬材料的加工及改性。但在低壓情況下，粒子和中性氣體的碰撞不足，無法使其達到平衡，在這種情況下我們認為放電所產(chǎn)生的負離子及電子的速度已經(jīng)偏離了Maxwellian分布。因此，流體力學的模擬方法已經(jīng)無法準確的解決此類問題。我們選用一種新的方法，即用運動分析的方法來解決低壓等離子反應器中的物理和化學過程，用包含大量粒子的模型來解決Boltzmann方程。粒子模型不僅可以解決刻蝕在基板上粒子的能量問題，而且還可以很好的解決粒子刻蝕率和各向異性的問題。Monte Carlo算法與PIC模擬方法的有機結(jié)合就形成了PIC_MCC模擬方法。在PIC_MCC中，我們假設中性氣體是在時間與空間位置上的一種特定分布。PIC_MCC模擬中采用的是一維空間（Zn），速度方向為三維（Vx，Vy，Vz）。1.2 PIC模型在PIC模型中（如圖a），粒子在電場力的作用下運動。粒子模擬只能解決少量粒子存在的模型，這個模型中的粒子數(shù)量遠遠小于真實情況下等離子體中的粒子數(shù)目，模擬中的每個粒子即超粒子表示106109個粒子。在模擬中我們必須有足夠多的超粒子，以減少粒子離散及噪聲擾動。超粒子與所劃分的網(wǎng)格點數(shù)之比必須大于1。在模擬中，我們主要解決Maxwell方程以及F=ma=q（E+v*B）。電場可以通過Maxwall方程解出。粒子在電場和磁場中所受的力可通過Newton_Lorentz方程解出。圖1.2.1 模擬中所用的系統(tǒng)模型（a），系統(tǒng)模型z軸方向上網(wǎng)格點的劃分。0點為接射頻電極極板端，ZN 為接地電極極板端。1.2.1 PIC模型中一個時間步循環(huán)的數(shù)值計算 圖1.2.2 PIC一個時間步的循環(huán)運算（1） 電荷密度j 被分配到每個網(wǎng)格點j上，這個過程稱之為電荷分配。因此，先從連續(xù)的節(jié)點Zpi 然后再到離散的節(jié)點Zj來計算電荷密度j 。電荷分配函數(shù)可用S（Zj-Zpi）來表示，包括零節(jié)點，第一節(jié)點和最后節(jié)點。圖1.2.3中所描述的是第一節(jié)點的電荷分配方式。這種分配方式把Zj 節(jié)點上的j單胞和Zj+1 節(jié)點上的（j+1）單胞這部分電荷看做帶電粒子或電荷云。這種帶電粒子看做是一種有限度的剛性電荷云，他們可以在通過彼此時不受束縛而自由運動，這種模型我們稱之為cloud_in_cell或者CIC.圖1.2.3如果電荷粒子的密度是一定的，j和j+1之間的距離為Z，那么電荷粒子qpi 分配給節(jié)點j的電量為： 分配給節(jié)點j+1的電量為： 因此，在Zpi 上的電荷粒子qpi 分配給j節(jié)點的電荷密度j 為： （2） 電荷密度可以用來計算網(wǎng)格點上的電場E。在靜電場模擬中，*E=- B/t0,所以E=-，由一維條件下的Poisson方程可以得到： 電場可以由以下公式計算得出：（3） E又按照函數(shù)S（Zj-Zpi）分配給網(wǎng)格點上的粒子。在一維的靜電場模型中，電場分配各網(wǎng)格點的電場為： 帶電粒子所受的電場力為F=qE，一維靜電場模型中： （4） 運動方程可以計算出帶電粒子新的位置和速度。 在一維靜電場模型中，可用以下方程代替上述運動方程： 因為帶電粒子的速度v和位置x是不能同時確定的，所以leap_frog算法要采用不同模式原則。 圖1.2.2 leap_frog算法的網(wǎng)格點劃分示意圖。 應注意到初始條件下帶電粒子在時間t=0時的速度是需要改變的，把v（0）處的速度V退回到v（-t/2）處，然后通過帶電粒子所受的電場力還計算t=0時的速度。（5） 檢查邊界條件，檢查粒子是否附著在極板上。初始時間t=0時粒子的位置和初始-t/2時的速度已經(jīng)給出，帶電粒子的密度也可通過計算得出。圖1.2.2中所描述的（1）到（5）只是重復的循環(huán)，直到等離子體達到收斂。1.2.2 邊界處等離子體粒子的模擬 在射頻放電產(chǎn)生等離子體的模擬中，不僅要考慮中心等離子體處的粒子行為，也要模擬邊界處，即鞘層處的粒子行為。 位勢方程的邊界條件可通過Gauss法則得出： S等離子體區(qū)域和上下兩極板的總面積 A0下極板（接射頻電壓的極板）的表面區(qū)域面積 AN上極板（接地電極極板）的表面區(qū)域面積 0下極板（接射頻電壓的極板）的帶電粒子密度 N上極板（接地電極極板）的帶電粒子密度 網(wǎng)格點的電勢可通過以下方程計算得出: j=1，2，N-1，N為所劃分的網(wǎng)點。 一維系統(tǒng)的邊界條件為： 1.3 MCC model1.3.1 無碰撞的模擬方法 PIC模擬方法是一種碰撞模型。即使在低電壓情況下帶電粒子和中性氣體的碰撞也對維持放電起著非常重要的作用。碰撞可以將PIC和MC兩種方法結(jié)合起來進行模擬運算。PIC模擬的是所有粒子在同一時間步的運動。而MC方法模擬的是在碰撞中一些隨機粒子的行為，只對每個時間步中的部分粒子的行為進行模擬，我們稱之為MCC模擬方法。1.41.5 電子和中性粒子的碰撞1.5.1 碰撞截面的數(shù)據(jù)我們認為中性氣體（Ar，CF4和N2）的分布是均勻的，它們的速度分配在室溫情況下（Tg=0.026 eV或者300K）是服從麥克斯韋分布的。因此，中性氣體粒子和電子（平均Tg2 eV）相比所具有的能量很小，我們認為它們是靜止的。碰撞截面我們用（g）來表示，g是粒子在碰撞之前的相對速度。在模型中所有中性粒子的碰撞截面數(shù)據(jù)和相應的閾值在表1.1中給出。Type of collisionReactionThreshold（eV）ArElastic scatteringe+Are+ArTotalelectronic excitatione+Are+Ar*11.5Ionizatione+Ar2e+Ar+15.8CF4Momentum transfere+CF4e+CF4Vibrational excitatione+CF4e+CF40.108Vibrational excitatione+CF4e+CF40.168Vibrational excitatione+CF4e+CF40.077Electronic excitatione+CF4e+CF4*7.54Electron attachmente+CF4F- +CF35Electron attachmente+CF4F+CF3-5Dissociatione+CF4e+F- +CF3+12Dissociative ionizatione+CF42e+F+CF3+16Neutral dissociatione+CF4e+F+CF312Neutral dissociatione+CF4e+2F+CF217Neutral dissociatione+CF4e+3F+CF18N2Momentum transfere+ N2e+ N2excitatione+ N2e+ N2* (Y)aIonizatione+ N22+ N2+(Y)b15.6Ionizatione+ N22+ N2+(Y)+ (B2)18.1.5.2 粒子碰撞后的速度計算方法 一般情況下，我們所考慮的是兩個均勻球體粒子的之間的彈性碰撞。粒子在碰撞之前，我們設它們的質(zhì)量和速度分別為，m和M，v和V，它們之間的相對速度為g=v-V。在不考慮一般情況下的損失，我們認為在碰撞之前的系統(tǒng)條件為：V=0，v=g。碰撞之后的速度為v ，V，g= v V。因為碰撞前后系統(tǒng)的總動量守恒，在這種條件下，可以計算出質(zhì)心的速度VCM （圖1.11）。 (2.31)(2.32) (2.34)因為在質(zhì)心坐標系中兩個粒子的初始動量大小相等，方向相反，所以， (2.36)在系統(tǒng)中，兩個粒子所受的力大小相等，方向相反。所以，碰撞后的動量也是大小相等，方向相反。碰撞之后粒子的速度方向仍然是平行的，但是卻偏離了一定角度，如圖1.11圖1.11 粒子在質(zhì)心坐標系下的運動軌跡。為散射角要求出粒子碰撞之后的速度，就要先求出散射角。如果可以求出散射角，那么可以通過以下公式來求出粒子碰撞之后的速度，(2.38) (2.39)， (1.17)在質(zhì)心坐標系中粒子運動軌跡所在的平面叫做碰撞平面。碰撞平面和參考平面之間的夾角為，參考平面是任意選取的，因此夾角為： (2.40)R是在特定分布0,1之間的一個隨機數(shù)。兩個粒子碰撞之后的速度有兩個主要影響因素，第一個是角度，第二個是b。圖1.1.2 碰撞參數(shù)和b系統(tǒng)要維持通量守恒，所以通過環(huán)面2db和立體角2sind的總粒子數(shù)相等。(2.41)上式中負號表示如果增加b，會相應的減小。 (2.42)db/d在這里去絕對值，因為其為負值。(2.43)散射角與兩粒子之間的電勢和速度有關系，不同的散射角，碰撞截面不同。（2.44）（2.45）不同的碰撞截面取決于粒子之間的相互作用力和電勢。例如，在電子和Ar的碰撞中，在計算電勢時可以忽略掉庫倫電勢的影響（彈性碰撞，激發(fā)，電離）。不同情況下的碰撞截面由以下公式給出，（2.46）=E/E0 是無量綱的能量，E是電子的相對能量，E0 原子的單位能量（E0 =27.21 eV）.散射角與隨機數(shù)R及之間的關系是： （2.47） 這樣就可以求出散射角，而粒子彈性碰撞之后的速度也可通過公式（1.15）和（1.16）求出。角度可通過公式（1.18）求出。多數(shù)情況下的電子和中性粒子碰撞中，公式（1.15）和（1.16）中的M+mM，gv。a) e_Ar excitatione+Are+Ar*根據(jù)動量和能量守恒公式，（2.48）（2.49）Eth 是非彈性碰撞的閾值。M+mM，gv(2.50)v為的標量，（2.51）E=mv2/2是電子碰撞之前的能量。激發(fā)過程如果看做是一個彈性碰撞過程，碰撞前的速度為何V。碰撞之后的速度可以由公式（1.15）和（1.16）求出。公式中所有的v用來代替。b) electron_A ionizationA是質(zhì)量為M的中性粒子。（2.52）A+ 為碰撞后的離子，e1 是初始電子，e2 是發(fā)射出的電子。能量守恒公式，（2.53）Eth 為電離反應的閾值。因為入射離子與電子的質(zhì)量比很大，我們可以認為電子的動量遠小于中性粒子的動量，電子與中性粒子碰撞之后，電子偏離，中性粒子變成離子，碰撞之后的軌跡仍然是未被干擾的。這種假設一方面使得在碰撞之前中性粒子的速度V=V；另一方面，能量守恒公式表示為，（2.54）等式的左邊我們看做是電離能量的改變量E，需要找到一種方法把散射電子和發(fā)射電子區(qū)分開。（2.55）（2.56）（2.57）Einc =mv2/2為初始電子的能量。1和0 單位是電子伏特。當E區(qū)分開后，我們可以從公式（1.33）中來求得v。（2.58）類似于激發(fā)情況，（2.59）碰撞之后散射電子的速度v可有公式（1.15）求出。式中v用來代替。碰撞之前，發(fā)射出的電子的能量實際是不存在的，我們假設為，（2.60）碰撞之后的速度可由公式（1.15）可算出，用ej 來代替v，v ej 來代替v。由于很多反應的碰撞截面數(shù)據(jù)是無法在文獻中查詢到的。如果沒有碰撞截面的數(shù)據(jù)，我們就無法求出相應的散射角?？偟呐鲎步孛婵捎晒剑?.21）得出T =4（g），這樣可以得到sindd/4。意味著碰撞之后的速度方向隨機，我們稱這種性質(zhì)的散射為各向同性。因為electron_CF4和electron_N2碰撞的碰撞截面數(shù)據(jù)是沒有的，因此我們認為它們之間的碰撞是各向同性的。各向同性的散射角在區(qū)間0, ，可表示為，（2.61）角度可從公式（1.18）求出，而electron_CF4和electron_N2 彈性碰撞之后的速度可由公式（1.15）和（1.16）求出。與electron_Ar的彈性碰撞相似，M+mM ，gv。1.6 Ion_neutral collisions1.6.1 Cross_section data(a) Collisions of Ar+ with neutrals因為所有Ar+ 參與的反應的碰撞截面數(shù)據(jù)均已給出，我們用null_collision方法來處理此類問題。1) 主程序讀取inp1d_mc.dat文件，并計算gas densitygas density=pertorr*pressure/gtemp2) 主程序打開data/e_Ar_eltot.dat文件，讀取e和Ar的碰撞截面數(shù)據(jù)3) 根據(jù)語句“species.eq.XX”來判斷主程序讀取哪些碰撞截面數(shù)據(jù)文件。(species.eq.argon/CF4)4) 主程序讀取CF4ion與CF4反應所需能量的文件（70，file=“data/dif_engy_cf3p.dat”）（70，file=“data/dif_engy_fm.dat”）（70，file=“data/dif_engy_cf3m.dat”）5) Spec(i)所代表的粒子分別為：spec（1）-e spec（2）-Ar spec（3）-CF3+ spec（4）-F- spec（5）-CF3-a) 判斷語句if (nflag_col_spec(1).eq.1)，由于此時只有e一種離子存在，因此調(diào)用子程序prob_e_ngas(),來模擬計算電子和所有中性氣體粒子的碰撞反應。if (nflag_col_spec(1).eq.1) thencall prob_e_ngas()end ifb) 對模擬中氣體種類進行判斷（(species.eq.XX），然后對不同的species以及spec（i）選擇調(diào)用不同的子程序來進行碰撞模擬。if(species.eq.argon.or.species.eq.argon/CF4.or.species.eq.argon/CF4/* N2) thenif (nflag_col_spec(2).eq.1) then call prob_Ar_Ar() end if end ifcc*cc probability of various ion-molecule collisions; beta_inf=3. polar_cf4=2.9d-30 polar_Ar =1.642d-30 if (species.eq.argon/CF4.or.species.eq.argon/CF4/N2) then do is=3,nspeccc*cc probability of CF4 ions + CF4 reactive and ellastic collisions if (nflag_col_spec(is).eq.1) then prob_col_spec(is)=e*sqrt(polar_cf4*PI*(mass(is)+ * mass(0)/mass(is)/mass(0)/EPS0)*beta_inf*beta_inf* * ratio_gas(2)*gden*dtis(is)cc*/cc probability of CF4 ions - Ar ellastic collisions prob_col_CFAr(is)=e*sqrt(polar_Ar*PI*(mass(is)+ * mass(2)/mass(is)/mass(2)/EPS0)*beta_inf*beta_inf* * ratio_gas(1)*gden*dtis(is) end if end docc*cc probability of Ar+ + CF4 ellastic collisions if (nflag_col_ij_2(3).eq.1) then prob_col_ArCF=e*sqrt(polar_cf4*PI*(mass(2)+mass(0)/mass(2) * /mass(0)/EPS0)*beta_inf*beta_inf*ratio_gas(2)*gden*dtis(2) end ifcc*cc probability of CF4 ions + CF4 reactive and ellastic collisions cc in case of pure CF4 discharge else if (species.eq.CF4) then do is=2,nspec if (nflag_col_spec(is).eq.1) then prob_col_spec(is)=e*sqrt(polar_cf4*PI*(mass(is)+ * mass(0)/mass(is)/mass(0)/EPS0)*beta_inf*beta_inf* * ratio_gas(2)*gden*dtis(is) end if end do end ifcc*/cc probability of (N2+) + N2 and (N2+) + Ar charge exchange if (species.eq.argon/CF4/N2) then call prob_Ar_N2 call prob_N2_N2 end ifcc*cc write all probabilities write(*,*)prob,(prob_col_spec(i),i=1,nspec) write(*,*) if (species.eq.argon/CF4.or.species.eq.argon/CF4/N2) then write(*,*)prob CF4 ions - Ar,(prob_col_CFAr(i),i=3,nspec) write(*,*)prob Ar+ - CF4,prob_col_ArCF write(*,*)prob Ar+ - N2 charge exchange,prob_col_ArN2 write(*,*)prob N2+ - Ar ch ex and el col,prob_col_N2Ar write(*,*) end ifcc*/cc end of initialization end ifcc*cc ELECTRON IMPACT COLLISIONS if (nflag_col_spec(1).eq.1) then if (species.eq.argon) then call mc_e_ngas(-1,2,-1,-1,-1,-1) else if (species.eq.argon/CF4) then call mc_e_ngas(0,2,3,4,5,-1) else if (species.eq.argon/CF4/N2) then call mc_e_ngas(0,2,3,4,5,6) else if (species.eq.CF4) then call mc_e_ngas(0,-1,2,3,4,-1) end if end ifcc*/cc ION IMPACT COLLISIONS if (nt_ion_i.eq.nt_ion0.or.nt_ion_i.eq.0) then if (species.eq.argon.or.species.eq.argon/CF4.or. * species.eq.argon/CF4/N2) thencc*/cc Argon ion impact collisions with Ar, CF4 or N2 if (nflag_col_spec(2).eq.1) then call mc_Ar_Ar() end if end ifcc*/cc N2+ + N2 and N2+ + Ar charge exchange if (species.eq.argon/CF4/N2) then call mc_N2_N2() end ifcc*/cc if Ar/CF4 discharge if (species.eq.argon/CF4.or.species.eq.argon/CF4/N2) thencc*/cc CF3+ ion impact collisions with CF4 and Ar if (nflag_col_spec(3).eq.1) then call mc_ion_neut(3,0,dif_engy_cf3p,n_reac_cf3p, * num_col_ij_3,extra_col_ij_3,rate_col_ij_3) call mc_cf_Ar(3,2,num_col_ij_3,rate_col_ij_3) end ifcc*/cc F- ion impact collisions with CF4 and Ar if (nflag_col_spec(4).eq.1) then call mc_ion_neut(4,0,dif_engy_fm,n_reac_fm, * num_col_ij_4,extra_col_ij_4,rate_col_ij_4) call mc_cf_Ar(4,2,num_col_ij_4,rate_col_ij_4) end ifcc*/cc CF3- ion impact collisions with CF4 and Ar if (nflag_col_spec(5).eq.1) then call mc_ion_neut(5,0,dif_engy_cf3m,n_reac_cf3m, * num_col_ij_5,extra_col_ij_5,rate_col_ij_5) call mc_cf_Ar(5,2,num_col_ij_5,rate_col_ij_5) end ifcc*/cc if species are only CF4 else if (species.eq.CF4) thencc*/cc CF3+ ion impact collisions if (nflag_col_spec(2).eq.1) then call mc_ion_neut(2,0,dif_engy_cf3p,n_reac_cf3p, * num_col_ij_3,extra_col_ij_3,rate_col_ij_3) end ifcc*/cc F- ion impact collisions if (nflag_col_spec(3).eq.1) then call mc_ion_neut(3,0,dif_engy_fm,n_reac_fm, * num_col_ij_4,extra_col_ij_4,rate_col_ij_4) end ifcc*/cc CF3- ion impact collisions if (nflag_col_spec(4).eq.1) then call mc_ion_neut(4,0,dif_engy_cf3m,n_reac_cf3m, * num_col_ij_5,extra_col_ij_5,rate_col_ij_5) end ifcc*/ end if end ifcc*/ return endcc*/cc end of MCMC1、Electron_neutral collisions1）Elastic collisione+Are+Are+CF4e+CF4e+ N2e+ N22)Inelastic collisiona) Excitationi. The mian reactionTotalelectronic excitation e+Are+Ar* 11.5Vibrational excitation e+CF4e+CF4 0.108Vibrational excitation e+CF4e+CF4 0.168Vibrational excitation e+CF4e+CF4 0.077Electronic excitation e+CF4e+CF4* 7.54Ratational excitation e+ N2e+N2 0.02Vibrational excitation e+ N2e+N2 0.291Vibrational excitation e+ N2e+N2 0.59Vibrational excitation e+ N2e+N2 0.88Vibrational excitation e+ N2e+N2 1.17Vibrational excitation e+ N2e+N2 1.47Vibrational excitation e+ N2e+N2 1.76Vibrational excitation e+ N2e+N2 2.06Vibrational excitation e+ N2e+N2 2.35Vibrational excitation e+ N2e+N2 0.290 e+ N2e+N2 * 6.17 6.785 7.35 7.36 7.80 8.16 8.40 8.55 8.89 11.03 11.88 12.35 13.00 ii. The main formula (2.48) (2.49)Where Eth is the threshold energy of the inelastic collision, i.e.the excitation in this case . Taking into consideration that M+mM and gv we obtain ( 2.50) (2.51)Where E=mv2/2 is the electron energy before collision. The excitation process is treated as if it were an elastic collision with pre_collision velocities and V. The post_collision velocities are given by Eqs.(2.37)-(2.39) in which all vs are replaced by . b) Ionization(electron_A)Where A denotes the neutral particle with a mass , which can be assumed equal to that of the ion M .The process is represented as (2.52)Where A+ denotes an ion,e1 is the incident electron , e2 is the ejected electron , and the symbols in the parentheses denote the velocities .The energy balance equation is (2.53)Where Eth is the threshold energy of the ionization . Because of the large ion_to_electron mass ratio ,we can assume that the momentum of the incident electron is much lower than the momentum of the neutral particle,i.e. the incident electron removes an electron from the neutral ,and the neutral becomes an ion,continuing on its trajectron undisturbed . On one hand ,this assumption means that the created ion takes the velocity and direction of the neutral particle before collision ,i.e. V=V .On the other hand , the energy balance can be rewritten (2.54)The left_hand side is known as the excess energyE after ionization and we need to find an algorithm how to divide it into the scattered and ejected electrons. In general, when there is no published work on the division of the excess energy , it is divided randomly into two (2.55) (2.56)For the electron_Ar ionzation ,however,we use the expression for sampling the Eej by a random number R,(2.57)Where Einc =mv2/2 is the energy of the incident electron. The units of 0 and1 are electronvolts .Once the excess energy is divided, v is calculated from Eq.(2.54) (2.58) (2.59)The post_collision velocity of the scattered electron v is obtained by Eq(2.37) in which v is replaced by .The pre_collision velocity for the ejected electron , which does not exist in reality, is assumed (2.60)and the post_collision velocity is obtained again from Eq.(2.37) by replacing v with and v with vej . For many types of collision no data on differential cross_section is assumed not to depend to the deflection angle x .The total cross_section calculated from Eq.(2.43)gives T=4（g）.Therefore , the probabilityDefined by (2.44)becomes sinxdxd/4.This means that the post_collision velocities take random direction. The scattering with this property is called isotropie.Since there are no sufficient data for the differential cross_sections of electron_CF4 and electron_N2 collisions we assume that the scattering is isotropic .The approximation o