PIC-MCC方法的计算步骤详解

前面两节分别介绍了PIC方法和MCC方法，下面将两种方法在完整的程序框架中进行整合。图5-38呈现了完整的PIC-MCC模拟仿真流程图。首先必须定义模型中的网格并且预置宏粒子，然后执行循环操作：

1）将网格中的宏粒子所带电荷分配到网格的节点处，获得电荷或者电流的分布状况。

2）基于网格节点电荷或者电流的分布求解电磁场方程，求解宏粒子之间的电磁场力。

3）在得到电磁场与宏粒子之间的相互作用关系后，利用经典或者相对论力学原理求解宏粒子的运动状态。

4）基于蒙特卡洛法模拟宏粒子之间的碰撞。

图5-38 PIC-MCC模拟仿真流程图

1.仿真设置

此处针对鞘层的仿真采用1d3v的PIC-MCC模型，即在相空间中坐标空间为一维而速度空间为三维。这种仿真模型可以较好地反映出鞘层的运动情况，最大限度地保留各种粒子的信息，同时由于采用一维位置空间，可以有效地控制计算量，节省计算时间，对于定性分析平板触头间的弧后等离子体现象具有较好的效果。

2.PIC模拟中参数的选取

为了保证PIC模拟过程中的数值计算稳定性，保证结果的收敛和可信，同时能够分辨出物理过程的细微变化，需要时间步长、网格尺寸满足如下条件：

Δx＜λ_D

Δt＜2/ω_p

Δt＜Δx/c

式中 Δx——网格宽度（m）；

Δt——仿真时间步长（s）；

λ_D——等离子体德拜长度（m）；

ω_p——等离子体振荡频率（rad/s）；

c——真空中的光速（m/s）。

对于静电模拟问题，可以消除上述时间步长和网格尺寸的约束，但是由于存在等离子体扩散问题，需要增加一个限制条件：

Δt＜Δx/v_thermal

式中 v_thermal——电子的平均热速度（m/s）。

此外，在使用PIC-MCC方法进行仿真，用蒙特卡洛碰撞方法模拟粒子间的碰撞时，还应控制时间步长以保证每个时间步长Δt内不发生太多次碰撞，即

νΔt1

式中 ν——碰撞发生的概率。(www.xing528.com)

网格宽度、时间步长等仿真参数的选择对仿真结果的收敛性、计算速度等有很大的影响。例如，如果Δx的设置远远大于德拜长度，则可能出现因为数值振荡导致所谓的“自加热”现象。在本节的仿真中，我们取网格宽度为以满足相应的约束条件，并留有一定的裕度。经过一定的尝试，根据仿真模型尺寸的不同，我们取时间步长为10^-12s，这样设定满足上文中的限制条件，同时也可以满足我们对运算速度的要求。同时，我们设定仿真的总时间为5μs，此时间长度可以完整地观察真空电弧弧后等离子体鞘层发展过程。

3.粒子初始条件

在仿真开始时，粒子被初始化，仿真中均设定铜原子、铜离子、电子在初始时刻均匀地分布于触头间隙之间，并具有一定的热运动速度。电子温度和离子温度选取标准的讨论参见5.2节，金属蒸气温度认为接近于触头表面温度，假设Cu金属蒸气温度为0.2eV。为了简化起见，仿真中认为全部铜离子均为正二价。

初始时刻的铜原子、铜离子、电子密度分别用n₀、n_i、n_e表示，为了满足初始时刻的等离子体电中性条件，粒子密度存在如下关系

n_e=2n_i （5-54）

另外我们将粒子边界条件设置为：当粒子运动至电极区域时即被电极吸收。

4.电磁场求解设置

由于我们所关心的物理问题中，电磁场的演化速度远小于等离子体特征时间，因此可以忽略感应磁场，采用静电场模型计算电场力，在每一时间步长内只需要求解Possion方程。采用此方法可以节省计算时间，提高计算效率。

仿真中电场的边界条件为Dirichlet型边界条件，即给定边界处的电势值。其中，左侧为阳极，给定电势为0；右侧为阴极，其电势为按固定斜率下降的瞬态恢复电压TRV。

5.MCC碰撞设置

在建立鞘层发展模型时，仅考虑铜离子与中性铜原子的碰撞，而不考虑电子的影响。因为电子相对于离子质量较小，鞘层的发展取决于离子的扩散速度，所以此处将电子的影响忽略。但是如果需要讨论鞘层发展阶段击穿等现象则需要把电子可能引起的激发和电离加以考虑。在我们的仿真模型中，考虑了铜离子与铜原子之间的电荷交换及弹性碰撞。

电荷交换碰撞：

Cu+Cu²⁺→Cu²⁺+Cu （5-55）

弹性碰撞：

Cu+Cu²⁺→Cu+Cu²⁺ （5-56）

电荷交换碰撞和弹性碰撞的碰撞截面如图5-39所示。

铜离子与铜原子弹性碰撞截面σ_ex为

σ_ex=（7.0-0.38lnε）²×10^-20 （5-57）

铜离子与铜原子电荷交换碰撞截面σ_mx为

σ_mx=（6.45-0.365lnε）²×10^-20 （5-58）

图5-39 Cu离子碰撞截面数据

6.后处理——鞘层边界位置的判定

通过仿真，我们可以得到触头间隙间的电势分布，电子、离子数量，电子、离子密度分布等参数随时间变化情况。然而要得到鞘层厚度随时间的变化情况还需根据仿真得到的电子、离子密度分布情况进行一定的后处理。在后处理过程中，我们认为，如果某一网格后的所有网格中的离子数均多于电子数，而此网格前一网格中离子数少于电子数，则此网格所在位置即为离子鞘层的边界，此网格到弧后阴极的距离即为鞘层厚度。通过后处理过程可以得到鞘层厚度随时间变化的情况，便于我们对不同参数下的鞘层厚度变化情况进行比较，并得到相应的鞘层变化规律。