气体电子崩发展的蒙特卡罗模拟模型

在气体放电的发展过程中，电子崩的发展主要受到电场力的作用。从宏观上看，电子在气体中的运动轨迹是曲线；从微观上来讲，电子在气体中的运动本身具有随机的性质。一个在电极表面或放电空间中产生的自由电子，如何及何时与气体分子碰撞都是偶然的，但是有一定的概率分布；而且电子同气体分子的碰撞过程本身也是一个随机事件，可能发生的碰撞有弹性碰撞、激发碰撞、电离碰撞、附着碰撞以及和多元分子还可能发生中性分解碰撞等，每种碰撞过程的发生都具有一定的概率分布。

这里以电子在平板均匀电场中运动的过程，来描述蒙特卡罗模拟法。模拟过程大致有三个过程，即源分布抽样过程（产生粒子的初始状态），空间、能量和方向的随机游动过程和记录分析结果的过程。首先介绍状态参数和状态序列的概念。

状态参数，电子在气体中运动的状态，可用一组参数来描述，包括电子的空间位置、能量、运动方向和时间：

S=（r，ε，Ω，t）（3-4）

式中，r是电子空间位置矢量，ε是电子能量，Ω是电子的运动方向，t是电子的运动时间，则碰撞后的状态参数如下所示：

Sn=（rn，εn，Ωn，tn）（3-5）

式中，rn是电子第n次碰撞点的位置，εn是电子第n次碰撞后的能量，Ωn是电子第n次碰撞后的运动方向，tn是电子在第n次碰撞时经历的时间。

状态序列，一个由源出发的电子，在气体中运动，一般要经历若干次碰撞（弹性、激发、电离、附着以及中性分解碰撞），电子在气体中的运动过程，可以用以下碰撞点的状态序列来描述：

S0，S1，…，Sn-1，Sn（3-6）

这里S0是从源发出的电子初始状态，Sn是终止状态。可以用一个矩阵来表示：

式（3-6）和式（3-7）可以称为电子在气体中的随机运动历史。

3.1.2.1 单一气体电子崩发展的蒙特卡罗模拟模型

运动过程的模拟和记录结果包括：①初始状态S0的确定；②确定下一碰撞点；③确定碰撞性质；④确定碰撞后的方向Ωn+1；⑤确定碰撞后的能量εn+1；⑥记录结果。

1）初始状态S0的确定

设电子的空间、能量和方向分布为

S（r0，ε0，cosα0）=S1（r0）S2（ε0）S3（cosα0）（3-8）

源分布常常归一到单位源强度，即

∭S（r0，ε0，cosα0）dr0dε0dcosα0=1（3-9）

表示源分布S（r0，ε0，cosα0）是一个概率密度函数，从源发射的一个电子，就是由源分布抽样得到的。

2）确定下一碰撞点

由于两次碰撞之间电子运动路径是弯向电场方向的抛物线，在电场作用下加速获得能量。这里对它的运动轨迹用折线来近似描述，用L表示折线的长度。已知状态参数Sn来确定下次碰撞点Sn+1，则第n+1次碰撞点的位置rn+1近似为

rn+1=rn+LΩn

xn+1=xn+Lun

yn+1=yn+Lvn

zn+1=zn+Lwn （3-10）

式中，un、vn、wn是Ωn的方向余弦，L为两次碰撞点的距离，它由分布

抽样确定。通常先由自由程分布

e-ρ，ρ≥0（3-12）

抽样ρ，再由关系式

解出L。ρ的抽样结果是

ρ=-ln R（3-14）

式中，R是均匀分布在[0，1]中的随机数。对于均匀介质，可以得出：

式中，∑t（εn）是电子能量为εn时，气体分子的宏观总截面。

3）确定碰撞性质

对于单一气体，电子与气体分子可能会发生弹性、振动、激发、附着、电离和中性分解等碰撞，这时需要判断发生哪种碰撞，设它的微观截面为Qm、Qev、Qex、Qat、Qion和Qnd，则气体分子的宏观总截面为

Qt（εn）=Qm（εn）+Qev（εn）+Qex（εn）+Qat（εn）+Qion（εn）+Qnd（εn） （3-16）

电子与气体发生各种性质碰撞的概率为

Pj=Qj（εn）/Qt（εn）

∑Pj=1（3-17）

式中，j代表气体的各种微观碰撞截面。用标准抽样法，确定气体发生碰撞的性质。

4）确定碰撞后的方向Ωn+1

如果能确定碰撞前后的方向夹角，即散射角θL和散射方位角φ，就可以确定Ωn+1的方向余弦。φ是碰撞后的方向Ωn+1偏离碰撞前Ωn所在的固定平面的角度，均匀分布在[0，2π]，则Ωn+1的方向余弦按如下确定：

式中，a=cosθL，b=sinθL=，c=cosφ，d=sinφ。

当1-→0时，不能应用上式，可应用下面简单公式：

un+1=sinθLcosφ

vn+1=sinθLsinφ

wn+1=cosθLwn （3-19）

5）确定碰撞后的能量εn+1

碰撞后的能量εn+1，一般与碰撞前后粒子运动方向间的夹角有关。电子与气体分子可能发生弹性、振动、激发、附着、电离和分解等碰撞。如果发生附着碰撞，则此电子的模拟结束；如果为弹性碰撞，碰撞后的能量εn+1为

式中，A是气体分子质量与电子质量之比，将式（3-21）代入式（3-20），得

电子质量me为9.109534×10-31kg，而气体分子质量M的数量级最小为10-27kg，由于A=M/me≈104≫1，上式可简化为

对于发生电离碰撞，则旧电子能量为

εn+1=（εn-Δεion）（1-R）（3-24）

新产生的二次电子的能量为

ε′n+1=（εn-Δεion）R（3-25）

式中，Δεion为发生电离碰撞的阈值能量，R是均匀分布在[0，1]中的随机数。而对于其他的碰撞：

εn+1=εn-Δεi（3-26）

式中，Δεi分别为振动、激发和分解等碰撞的阈值能量。

以上模拟电子的运动过程，可分为两大步：第一步是由源分布来确定电子的初始状态S0；第二步是由Sn来确定Sn+1。第二步又可分为两个过程：第一是在空间位置上的变化，称为空间输运过程；第二是由εn、cosαn来确定εn+1、cos αn+1，称为碰撞过程。以后就重复这两种过程，直到模拟结束。

6）记录结果

对于每一个电子运动的历史过程，如式（3-6）和式（3-7）所示，它可能运动到模拟的截止时间，也可能被附着。对模拟过程中需要得到的量，按其状态的空间、能量、方向和时间特征分别记录下来，将这些记录结果按要求处理即可。

3.1.2.2 混合气体电子崩发展的蒙特卡罗模拟模型

对于混合气体电子崩发展的蒙特卡罗模拟模型，与单一气体模型相比基本相同，不同之处是：首先判断与哪种气体分子发生碰撞，然后再判断发生哪种性质的碰撞；其余部分一致。假如是由气体1和2组成的二元混合气体，有

和分别是气体分子1和2的宏观总截面，定义如下：

式中，N是气体i（i=1，2）的分子数密度，是气体i（i=1，2）的微观总截面，即弹性、振动、激发、电离、附着和中性分解等碰撞截面的总和。这样对于二元混合气体，电子与气体分子1和2发生碰撞的概率分别是：

用标准抽样法，确定与哪种气体发生碰撞。

电子与气体分子可能会发生弹性、振动、激发、附着、电离和中性分解等碰撞，这时需要判断碰撞的性质，设它们的微观截面为，则

式中，i（i=1，2）分别代表气体1和2。电子与气体i发生各种性质碰撞的概率为

式中，j代表气体的各种碰撞截面。用标准抽样法，确定气体发生碰撞的性质。

3.1.2.3 单一气体电子崩发展的蒙特卡罗模型的建立

本节结合图3-1所示的蒙特卡罗模拟气体放电的流程图，具体说明此方法程序的编写。采用C语言编写，运行速度快。

图3-1 单一气体计算流程图

在主程序中包含在[0，1]中产生均匀分布的随机数子程序和插值子程序。以下所用的Rm（Rm，m=1，2，3，…）是均匀分布在[0，1]上的随机数，在下面各公式中不再说明。

1）模拟电子的初始化

模拟过程假定气体中带电粒子的密度都很低，从而可以不考虑带电粒子之间的库仑作用，分子之间不发生化学作用。模拟气体放电条件：压强p为133.322Pa，温度为20℃，气体分子数密度N=3.29×1022m-3[由式（3-32）计算得到]，电场E方向与z轴相反

式中，p是压强（单位为Pa），k B是玻耳兹曼常数，k B=1.38×10-23J/K，T是温度（单位为K）。

首先判断此电子是初始电子还是二次电子。如果是初始电子，假设在t=0时，能量为ε0=1.0e V的n0个初始电子以各向同性分布从阴极表面出发，出发点设为坐标原点（x0=0，y0=0，z0=0）。这些初始电子进入电场的方向角α和φ随机决定于以下余弦分布：

cosα=2R1-1（3-33）

φ=2πR2（3-34）

φ均匀分布在[0，2π]上。在实验情况下，由于初始电子只能进入两电极间隙中，对于在z轴的分速度与z轴方向相反的初始电子可以舍去，则初始电子的初始速度为

式中，ε0是电子的初始能量，e/m是荷质比。

如果是二次电子，则将产生电离碰撞后生成的二次新电子位置、能量、方向和时间记录下来，直接作为其初始参数。

对于选定每种气体的碰撞截面数据，以矩阵或拟合多项式的形式出现在程序里，对点截面根据能量ε用插值法来得到相应的截面数据。

2）空碰撞技术（null technique）

电子在电场中飞行时间有四种计算方法：①自由飞行时间；②平均碰撞距离模型；③平均碰撞时间模型；④空碰撞技术。这里采用了空碰撞技术，此方法是对前三种方法的改进，提高计算速度。在电子的运动过程中，不仅发生真实碰撞，还发生空碰撞。在确定了气体的碰撞截面组后，由电子的能量ε就可得到对应的截面数据和速度，按照下式就可以找到碰撞频率的最大值υmax：

υmax=max[NQt（ε）|V（ε）|]（3-36）

式中，N是气体总的分子数密度，Qt（ε）是气体分子的宏观总截面，V（ε）是电子速度，ε是电子能量。那么自由飞行时间可由下式得到：

Δt=-（ln R3/υmax）（3-37）

电子在z轴方向受到电场力的作用，在z轴方向做匀加速运动，从电场中获得能量。电子在自由飞行时间Δt内遵守能量守恒定律，运动到第n+1个Δt时间点的位移、速度和能量由以下方程决定：

式中，a为电子的加速度，E是电场强度，电子电荷量e=1.60207×10-19C，ε是电子能量（单位e V，1e V=1.602189×10-19J）。

3）确定碰撞类型

在电子每运动完Δt时间，这时判断电子与气体分子发生真实或空碰撞。对于单一气体，即由下式所表示的概率判断：

υreal=N|V（ε）|Qt（ε）

Preal=υreal/υmax （3-39）

式中，υreal是真实碰撞频率，Preal是发生真实碰撞的概率。如果Preal＞R4，则发生真实碰撞，否则为空碰撞。如果发生空碰撞，电子将保持运动方向不发生改变，继续下一个Δt时间的运动，再重复上面碰撞过程的判断。

电子与气体分子可以发生多种碰撞，如弹性、振动、激发、电离、附着和分解等碰撞，气体分子的宏观总截面为

Qt（ε）=Qm（ε）+Qev（ε）+Qex（ε）+Qat（ε）+Qion（ε）+Qnd（ε）（3-40）

由下式来判断碰撞的性质：

Pj=Qj（ε）/Qt（ε）

∑Pj=1 （3-41）

式中，Qj（ε）分别代表气体的第j个微观碰撞截面（j=1，2，3，…，包括弹性、激发、振动、附着、电离和中性分解等碰撞）。如果

P1+P2+…+Pj-1≤R5≤P1+P2+…+Pj（3-42）

则第j个碰撞过程发生。

4）碰撞后的运动方向和能量

电子发生碰撞后，电子运动方向和能量会发生改变，如果发生附着碰撞，此电子的模拟结束。对于碰撞后方向，在各向同性的情况下，各向同性散射角余弦x=cosα遵从如下分布：

对于各向同性散射角余弦分布的直接抽样方法如下：

各向同性散射方位角余弦x=cosφ遵从如下分布：

对此分布的直接抽样方法如下：

cosφ=cos2πR7

φ=2πR7 （3-46）

这时可以确定各向同性散射的方向Ω=ui+vj+wk，根据下面公式：

u=sinαcosφ

v=sinαsinφ

w=cosα （3-47）

碰撞之后再按照相应的算法确定带电粒子的能量和速度。如果发生附着碰撞（电负性气体），则电子模拟过程结束。如果发生弹性碰撞，电子的能量几乎不变，则按下式计算碰撞后电子的能量：

如果发生非弹性碰撞，则按下式计算碰撞后电子的能量：

εl=ε-Δεj（3-49）

式中，ε是碰撞前电子能量，Δεj是每一非弹性碰撞的阈值能量，εl是碰撞后的能量。对于发生电离碰撞，则旧电子能量为

εl=（εn-Δεion）（1-R8）（3-50）

新产生的二次电子的能量为

ε′l=（εn-Δεion）R8（3-51）

式中，Δεion为发生电离碰撞的阈值能量。将二次电子此时的位置、能量、速度和时间用存储数组记录下来。旧电子模拟结束后，接着模拟其产生的二次电子，这样可以减少数据存储空间。电子碰撞后的速度及在x、y、z轴上的速度分量由下式得出：

在每一步的模拟过程中，将所需的变量记录下来。最后将电子的运动时间累加起来，达到预定的模拟截止时间，则此电子的模拟过程结束。当所有旧电子和新产生的二次电子的运动时间达到预定的模拟截止时间，整个的模拟过程结束。

5）记录数据的抽样和计算

所有的电子模拟结束后，要对记录的结果抽样处理。不同的实验方法对应不同的抽样方法。常用方法有SST实验、TOF实验和PT实验的处理方法，而这里是针对PT实验的模拟，所以给出PT实验对结果的处理方法，电子的特性和它在电子崩中的位置有关，按下式抽样：

式中，ξi是t时刻对电子崩中的第i个电子的抽样，Nt是t时刻电子崩中的电子总数。电离系数密度比α/N和附着系数密度比η/N按下式计算：

有效电离系数密度比为（α-η）/N。

受到电场平衡区域的影响，虽然电子崩中的电子以自己的速度运动，但整个电子崩作为一个整体以漂移速度Ve沿电场方向运动。漂移速度Ve定义为电子崩中所有电子的平均速度。按下式计算：

上述几个式子中，α和η是电离和附着系数，N是气体分子数密度，是在取样时间t内所有电子的平均位移，n+和n-是在时刻t产生的正离子和负离子数。在实际计算过程中，n+、n-和的取值如下：

（1）n+，当电子与气体分子发生电离碰撞时，气体分子电离出一个正离子和一个电子，则将产生的正离子数累加，就是n+的值。

（2）n-，当电子与气体分子发生附着碰撞时，电子被气体分子吸附形成负离子，则将发生附着碰撞的次数累加，就是n-的值。

（3），将所有达到预定模拟截止时间的初始电子和新产生的二次电子的总数记录下来，将它们的位移累加，那么就等于位移和除以电子总数。

3.1.2.4 混合气体电子崩发展的蒙特卡罗模型的建立

在单一气体电子崩发展的蒙特卡罗模型的基础上，阐述混合气体电子崩发展的蒙特卡罗模型的建立，流程图如图3-2所示。以下所用的Rm（Rm，m=1， 2，3，…）是均匀分布在[0，1]上的随机数，在下面各公式中不再说明。

从图3-2可以看出在判断是否发生真实碰撞之后，首先判断与哪种气体分子发生碰撞，然后再判断发生哪种性质的碰撞，其他模拟流程与单一气体的相同，本节只给出不同之处，其余相同部分不再重复。在确定了混合气体的碰撞截面组后，由电子的能量ε就可得到对应的截面数据和速度，按照下式就可以找到碰撞频率的最大值υmax：

图3-2 混合气体计算流程图

Δt=-（ln R3/υmax）（3-58）

在电子每运动完Δt时间，这时判断电子与气体分子发生真实或空碰撞。对于二元混合气体，在模拟过程中对碰撞过程的处理分三步进行，首先由式（3-58）所表示的概率判断是否发生真实碰撞：

式中，υreal是真实碰撞频率，Preal是发生真实碰撞的概率。如果Preal＞R4，则发生真实碰撞，否则为空碰撞。如果发生空碰撞，电子将保持运动方向不发生改变，继续下一个Δt时间的运动。

如果发生真实碰撞，则第二步确定和哪种气体分子发生碰撞：

由下面式子来判断和哪种气体分子发生碰撞的概率：

P1=υm，1/υmax

P2=υm，2/υmax

P1+P2=1 （3-61）

式中，υm，1和υm，2分别是气体1和2的最大碰撞频率，υmax为υm，1和υm，2之和。如果P1＞R5，则跟气体1发生碰撞，否则跟气体2发生碰撞。

第三步确定电子与气体分子发生哪种碰撞，如弹性、振动、激发、电离、附着和中性分解等碰撞，每种气体分子的宏观总截面为

式中，分别为气体1和2的宏观总截面。由下式来判断发生碰撞的性质：

Pi，1+Pi，2+…+Pi，j-1≤R6≤Pi，1+Pi，2+…+Pi，j（3-64）

则第j个碰撞过程发生。