PICMCC程序手册X.docx

PICMCC程序手册X.docx

《PICMCC程序手册X.docx》由会员分享，可在线阅读，更多相关《PICMCC程序手册X.docx（57页珍藏版）》请在冰点文库上搜索。

PICMCC程序手册X

PIC-MCC程序手册

1.PIC_MCC的模拟方法和数值计算

1.1PIC_MCC的模拟原理

1.2PIC模型

1.3MCC模型

1.4模拟中所涉及的放电粒子

1.5模拟中的碰撞

2.数据分析

2.1程序正常运行所需文件及文件意义

2.2主要输出文件的意义

1PIC_MCC的模拟方法和数值计算

1.1PIC_MCC的模拟原理

目前在计算机模拟中大量采用低压射频放电模型来模拟材料的加工及改性。

但在低压情况下，粒子和中性气体的碰撞不足，无法使其达到平衡，在这种情况下我们认为放电所产生的负离子及电子的速度已经偏离了Maxwellian分布。

因此，流体力学的模拟方法已经无法准确的解决此类问题。

我们选用一种新的方法，即用运动分析的方法来解决低压等离子反应器中的物理和化学过程，用包含大量粒子的模型来解决Boltzmann方程。

粒子模型不仅可以解决刻蚀在基板上粒子的能量问题，而且还可以很好的解决粒子刻蚀率和各向异性的问题。

MonteCarlo算法与PIC模拟方法的有机结合就形成了PIC_MCC模拟方法。

在PIC_MCC中，我们假设中性气体是在时间与空间位置上的一种特定分布。

PIC_MCC模拟中采用的是一维空间（Zn），速度方向为三维（Vx，Vy，Vz）。

1.2PIC模型

在PIC模型中（如图a），粒子在电场力的作用下运动。

粒子模拟只能解决少量粒子存在的模型，这个模型中的粒子数量远远小于真实情况下等离子体中的粒子数目，模拟中的每个粒子即超粒子表示106—109个粒子。

在模拟中我们必须有足够多的超粒子，以减少粒子离散及噪声扰动。

超粒子与所划分的网格点数之比必须大于1。

在模拟中，我们主要解决Maxwell方程以及F=ma=q（E+v*B）。

电场可以通过Maxwall方程解出。

粒子在电场和磁场中所受的力可通过Newton_Lorentz方程解出。

图1.2.1模拟中所用的系统模型（a），系统模型z轴方向上网格点的划分。

0点为接射频电极极板端，ZN为接地电极极板端。

1.2.1PIC模型中一个时间步循环的数值计算

图1.2.2PIC一个时间步的循环运算

（1）电荷密度ρj被分配到每个网格点j上，这个过程称之为电荷分配。

因此，先从连续的节点Zpi然后再到离散的节点Zj来计算电荷密度ρj。

电荷分配函数可用S（Zj-Zpi）来表示，包括零节点，第一节点和最后节点。

图1.2.3中所描述的是第一节点的电荷分配方式。

这种分配方式把Zj节点上的j单胞和Zj+1节点上的（j+1）单胞这部分电荷看做带电粒子或电荷云。

这种带电粒子看做是一种有限度的刚性电荷云，他们可以在通过彼此时不受束缚而自由运动，这种模型我们称之为cloud_in_cell或者CIC.

图1.2.3

如果电荷粒子的密度是一定的，j和j+1之间的距离为△Z，那么电荷粒子qpi分配给节点j的电量为：

分配给节点j+1的电量为：

因此，在Zpi上的电荷粒子qpi分配给j节点的电荷密度ρj为：

（2）电荷密度可以用来计算网格点上的电场E。

在静电场模拟中，▽*E=-əB/ət≈0,所以E=-▽Φ，由一维条件下的Poisson方程可以得到：

电场可以由以下公式计算得出：

（3）E又按照函数S（Zj-Zpi）分配给网格点上的粒子。

在一维的静电场模型中，电场分配各网格点的电场为：

带电粒子所受的电场力为F=qE，一维静电场模型中：

（4）运动方程可以计算出带电粒子新的位置和速度。

在一维静电场模型中，可用以下方程代替上述运动方程：

因为带电粒子的速度v和位置x是不能同时确定的，所以leap_frog算法要采用不同模式原则。

图1.2.2leap_frog算法的网格点划分示意图。

应注意到初始条件下带电粒子在时间t=0时的速度是需要改变的，把v（0）处的速度V退回到v（-△t/2）处，然后通过带电粒子所受的电场力还计算t=0时的速度。

（5）检查边界条件，检查粒子是否附着在极板上。

初始时间t=0时粒子的位置和初始-△t/2时的速度已经给出，带电粒子的密度也可通过计算得出。

图1.2.2中所描述的

（1）到（5）只是重复的循环，直到等离子体达到收敛。

1.2.2边界处等离子体粒子的模拟

在射频放电产生等离子体的模拟中，不仅要考虑中心等离子体处的粒子行为，也要模拟边界处，即鞘层处的粒子行为。

位势方程的边界条件可通过Gauss法则得出：

S——等离子体区域和上下两极板的总面积

A0——下极板（接射频电压的极板）的表面区域面积

AN——上极板（接地电极极板）的表面区域面积

σ0——下极板（接射频电压的极板）的带电粒子密度

σN——上极板（接地电极极板）的带电粒子密度

网格点的电势可通过以下方程计算得出:

j=1，2，……，N-1，N为所划分的网点。

一维系统的边界条件为：

1.3MCCmodel

1.3.1无碰撞的模拟方法

PIC模拟方法是一种碰撞模型。

即使在低电压情况下带电粒子和中性气体的碰撞也对维持放电起着非常重要的作用。

碰撞可以将PIC和MC两种方法结合起来进行模拟运算。

PIC模拟的是所有粒子在同一时间步的运动。

而MC方法模拟的是在碰撞中一些随机粒子的行为，只对每个时间步中的部分粒子的行为进行模拟，我们称之为MCC模拟方法。

1.4电子和中性粒子的碰撞

1.4.1碰撞截面的数据

我们认为中性气体（Ar，CF4和N2）的分布是均匀的，它们的速度分配在室温情况下（Tg=0.026eV或者300K）是服从麦克斯韦分布的。

因此，中性气体粒子和电子（平均Tg>2eV）相比所具有的能量很小，我们认为它们是静止的。

碰撞截面我们用σ（g）来表示，g是粒子在碰撞之前的相对速度。

在模型中所有中性粒子的碰撞截面数据和相应的阈值在表1.1中给出。

Typeofcollision

Reaction

Threshold（eV）

Ar

Elasticscattering

e+Ar→e+Ar

Totalelectronicexcitation

e+Ar→e+Ar*

11.5

Ionization

e+Ar→2e+Ar+

15.8

CF4

Momentumtransfer

e+CF4→e+CF4

Vibrationalexcitation

e+CF4→e+CF4

0.108

Vibrationalexcitation

e+CF4→e+CF4

0.168

Vibrationalexcitation

e+CF4→e+CF4

0.077

Electronicexcitation

e+CF4→e+CF4*

7.54

Electronattachment

e+CF4→F-+CF3

5

Electronattachment

e+CF4→F+CF3-

5

Dissociation

e+CF4→e+F-+CF3+

12

Dissociativeionization

e+CF4→2e+F+CF3+

16

Neutraldissociation

e+CF4→e+F+CF3

12

Neutraldissociation

e+CF4→e+2F+CF2

17

Neutraldissociation

e+CF4→e+3F+CF

18

N2

Momentumtransfer

e+N2→e+N2

excitation

e+N2→e+N2*（Y）a

…

Ionization

e+N2→2+N2+（Y）b

15.6

Ionization

e+N2→2+N2+（Y）+（B2Σ）

18.

1.4.2粒子碰撞后的速度计算方法

一般情况下，我们所考虑的是两个均匀球体粒子的之间的弹性碰撞。

粒子在碰撞之前，我们设它们的质量和速度分别为，m和M，v和V，它们之间的相对速度为g=v-V。

在不考虑一般情况下的损失，我们认为在碰撞之前的系统条件为：

V=0，v=g。

碰撞之后的速度为v’，V’，g’=v’–V’。

因为碰撞前后系统的总动量守恒，在这种条件下，可以计算出质心的速度VCM（图1.11）。

（2.31）

（2.32）

（2.34）

因为在质心坐标系中两个粒子的初始动量大小相等，方向相反，所以，

（2.36）

在系统中，两个粒子所受的力大小相等，方向相反。

所以，碰撞后的动量也是大小相等，方向相反。

碰撞之后粒子的速度方向仍然是平行的，但是却偏离了一定角度χ，如图1.11

图1.11粒子在质心坐标系下的运动轨迹。

χ为散射角

要求出粒子碰撞之后的速度，就要先求出散射角χ。

如果可以求出散射角χ，那么可以通过以下公式来求出粒子碰撞之后的速度，

（2.38）

（2.39）

，

（1.17）

在质心坐标系中粒子运动轨迹所在的平面叫做碰撞平面。

碰撞平面和参考平面之间的夹角为ψ，参考平面是任意选取的，因此夹角ψ为：

（2.40）

R是在特定分布[0,1]之间的一个随机数。

两个粒子碰撞之后的速度有两个主要影响因素，第一个是角度ψ，第二个是b。

图1.1.2碰撞参数ψ和b

系统要维持通量守恒，所以通过环面2πdb和立体角2πsinχdχ的总粒子数相等。

（2.41）

上式中负号表示如果增加b，χ会相应的减小。

（2.42）

db/dχ在这里去绝对值，因为其为负值。

（2.43）

散射角χ与两粒子之间的电势和速度有关系，不同的散射角，碰撞截面不同。

（2.44）

（2.45）

不同的碰撞截面取决于粒子之间的相互作用力和电势。

例如，在电子和Ar的碰撞中，在计算电势时可以忽略掉库伦电势的影响（弹性碰撞，激发，电离）。

不同情况下的碰撞截面由以下公式给出，

（2.46）

ε=E/E0是无量纲的能量，E是电子的相对能量，E0原子的单位能量（E0=27.21eV）.散射角χ与随机数R及ε之间的关系是：

（2.47）

这样就可以求出散射角χ，而粒子弹性碰撞之后的速度也可通过公式（1.15）和（1.16）求出。

角度ψ可通过公式（1.18）求出。

多数情况下的电子和中性粒子碰撞中，公式（1.15）和（1.16）中的M+m≈M，g≈v。

a）e_Arexcitation

e+Ar→e+Ar*

根据动量和能量守恒公式，

（2.48）

（2.49）

Eth是非弹性碰撞的阈值。

M+m≈M，g≈v

（2.50）

v’为

的标量，

（2.51）

E=mv2/2是电子碰撞之前的能量。

激发过程如果看做是一个弹性碰撞过程，碰撞前的速度为

何V。

碰撞之后的速度可以由公式（1.15）和（1.16）求出。

公式中所有的v’用

来代替。

b）electron_Aionization

A是质量为M的中性粒子。

（2.52）

A+为碰撞后的离子，e1是初始电子，e2是发射出的电子。

能量守恒公式，

（2.53）

Eth为电离反应的阈值。

因为入射离子与电子的质量比很大，我们可以认为电子的动量远小于中性粒子的动量，电子与中性粒子碰撞之后，电子偏离，中性粒子变成离子，碰撞之后的轨迹仍然是未被干扰的。

这种假设一方面使得在碰撞之前中性粒子的速度V’=V；另一方面，能量守恒公式表示为，

（2.54）

等式的左边我们看做是电离能量的改变量△E，需要找到一种方法把散射电子和发射电子区分开。

（2.55）

（2.56）

（2.57）

Einc=mv2/2为初始电子的能量。

α1和α0单位是电子伏特。

当△E区分开后，我们可以从公式（1.33）中来求得v’。

（2.58）

类似于激发情况，

（2.59）

碰撞之后散射电子的速度v’可有公式（1.15）求出。

式中v用

来代替。

碰撞之前，发射出的电子的能量实际是不存在的，我们假设为，

（2.60）

碰撞之后的速度可由公式（1.15）可算出，用

ej来代替v，v’ej来代替v’。

由于很多反应的碰撞截面数据是无法在文献中查询到的。

如果没有碰撞截面的数据，我们就无法求出相应的散射角Χ。

总的碰撞截面可由公式（1.21）得出σT=4πσ（g），这样可以得到sinΧdΧdψ/4π。

意味着碰撞之后的速度方向随机，我们称这种性质的散射为各向同性。

因为electron_CF4和electron_N2碰撞的碰撞截面数据是没有的，因此我们认为它们之间的碰撞是各向同性的。

各向同性的散射角在区间[0,π]，可表示为，

（2.61）

角度ψ可从公式（1.18）求出，而electron_CF4和electron_N2弹性碰撞之后的速度可由公式（1.15）和（1.16）求出。

与electron_Ar的弹性碰撞相似，M+m≈M，g≈v。

1.5Ion_neutralcollisions

1.5.1Cross_sectiondata

（a）CollisionsofAr+withneutrals

因为所有Ar+参与的反应的碰撞截面数据均已给出，我们用null_collision方法来处理此类问题。

1）主程序读取inp1d_mc.dat文件，并计算gasdensity

gasdensity=pertorr*pressure/gtemp

2）主程序打开data/e_Ar_eltot.dat文件，读取e和Ar的碰撞截面数据

3）根据语句“species.eq.XX”来判断主程序读取哪些碰撞截面数据文件。

（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一种离子存在，因此调用子程序prob_e_ngas（）,来模拟计算电子和所有中性气体粒子的碰撞反应。

if（nflag_col_spec

（1）.eq.1）then

callprob_e_ngas（）

endif

b）对模拟中气体种类进行判断（（species.eq.XX），然后对不同的species以及spec（i）选择调用不同的子程序来进行碰撞模拟。

if（species.eq.'argon'.or.species.eq.'argon/CF4'.or.species.eq.'argon/CF4/*N2'）then

if（nflag_col_spec

（2）.eq.1）then

callprob_Ar_Ar（）

endif

endif

cc********************************************************\\\

ccprobabilityofvariousion-moleculecollisions;

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

dois=3,nspec

cc********************************************************\\\

ccprobabilityofCF4ions+CF4reactiveandellasticcollisions

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********************************************************///

ccprobabilityofCF4ions-Arellasticcollisions

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）

endif

enddo

cc********************************************************\\\

ccprobabilityofAr++CF4ellasticcollisions

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）

endif

cc********************************************************\\\

ccprobabilityofCF4ions+CF4reactiveandellasticcollisions

ccincaseofpureCF4discharge

elseif（species.eq.'CF4'）then

dois=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）

endif

enddo

endif

cc********************************************************///

ccprobabilityof（N2+）+N2and（N2+）+Archargeexchange

if（species.eq.'argon/CF4/N2'）then

callprob_Ar_N2

callprob_N2_N2

endif

cc********************************************************\\\

ccwriteallprobabilities

write（*,*）'prob',（prob_col_spec（i）,i=1,nspec）

write（*,*）

if（species.eq.'argon/CF4'.or.species.eq.'argon/CF4/N2'）then

write（*,*）'probCF4ions-Ar',（prob_col_CFAr（i）,i=3,nspec）

write（*,*）'probAr+-CF4',prob_col_ArCF

write（*,*）'probAr+-N2chargeexchange',prob_col_ArN2

write（*,*）'probN2+-Archexandelcol',prob_col_N2Ar

write（*,*）

endif

cc********************************************************///

ccendofinitialization

endif

cc********************************************************\\\

ccELECTRONIMPACTCOLLISIONS

if（nflag_col_spec

（1）.eq.1）then

if（species.eq.'argon'）then

callmc_e_ngas（-1,2,-1,-1,-1,-1）

elseif（species.eq.'argon/CF4'）then

callmc_e_ngas（0,2,3,4,5,-1）

elseif（species.eq.'argon/CF4/N2'）then

callmc_e_ngas（0,2,3,4,5,6）

elseif（species.eq.'CF4'）then

callmc_e_ngas（0,-1,2,3,4,-1）

endif

endif

cc********************************************************///

ccIONIMPACTCOLLISIONS

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'）then

cc********************************************************///

ccArgonionimpactcollisionswithAr,CF4orN2

if（nflag_col_spec

（2）.eq.1）then

callmc_Ar_Ar（）

endif

endif

cc********************************************************///

ccN2++N2andN2++Archargeexchange

if（species.eq.'argon/CF4/N2'）then

callmc_N2_N2（）

endif

cc********************************************************///

ccifAr/CF4discharge

if（species.eq.'argon/CF4'.or.species.eq.'argon/CF4/N2'）then

cc*