近年来，科学计算的验证和确认逐渐成为一个成熟的研究领域。美国航空航天协会AIAA关于计算流体力学的验证和确认指南被视为这一领域的“金标准”^[1]。美国桑迪亚国家实验室进一步指出，为满足工程物理领域的应用需求，需要专注于通过验证和确认的方法来提高科学计算的预测能力成熟度^[2]。验证是确认过程的前提：前者指是否正确求解方程组，后者指模拟结果是否描述了正确的物理问题^[1]。相比计算流体力学领域，计算放电等离子体领域的验证研究仍处在不断发展和完善的阶段，研究者对于这一环节的具体实施过程呈现出明显的个性化^[3-8]。一方面，不同严格程度的方法被用于代码验证的实践，包括简单测试，代码间比较，离散误差评估，收敛测试，精度阶测试等。能量守恒性检验作为一种简单测试方法，通常被用来验证粒子模拟算法在长时间演化时误差传播的性质^[9]。代码间比较在数值模拟中用于验证时较容易实现，因而也是最常用的方法^[10-14]。例如Smith等人对激光等离子体中的离子加速问题比较了3款PIC程序：EPOCH、LSP和WarpX并提供了相应的输入文件以鼓励更多同行从事代码间交叉比较^[10]。B. Bagheri等人针对空气中的正极性流注放电模型，采用6种不同模拟方法的计算并比较结果^[11]。Riva等人^[4]和Tranquilli等人^[6]分别构造了1D1V和2D2V的Vlasov-Poisson系统的精确解。另一方面，在解验证中由于粒子离散过程通常伴随着统计误差，传统的里查德森外推技术并不能直接用于粒子模拟的验证。Radtke等人^[15]因此提出了针对随机模拟方法的“鲁棒”外推技术，其核心思想是通过系综平均给出均值和方差，基于多元回归拟合并进行统计检验。本文提出粒子云网格（PIC）蒙特卡罗（MC）模拟的分层级验证方法，基于分层级验证方法可以准确定位PIC程序中错误模块，提升PIC模块的验证速度并对PIC程序的研制具有指导作用。

1 粒子云网格蒙特卡罗方法的特点

当前，科学计算的验证主要面向的是具有确定性偏微分方程组的网格离散方法。但PIC方法是一种粒子-网格离散的耦合方法，其验证过程呈现出明显不同的特点，主要体现在以下几点。

第一，验证过程的离散域。PIC方法的离散网格除了与一般数学模型相同的空间t和时间x的离散，还包括粒子的离散：离散量通常称为宏粒子的权重w。作为对比，等离子体直接动力学方法的离散网格为相空间和时间，广义离散网格的无量纲形式h表示如下

$ {h^p} = {\Bigg(\dfrac{{\Delta {x_{{\text{coarse}}}}}}{{\Delta {x_{{\text{fine}}}}}}\Bigg)^{{p_x}}} = {\Bigg(\dfrac{{\Delta {v_{{\text{coarse}}}}}}{{\Delta {v_{{\text{fine}}}}}}\Bigg)^{{p_v}}} = {\Bigg(\dfrac{{\Delta {t_{{\text{coarse}}}}}}{{\Delta {t_{{\text{fine}}}}}}\Bigg)^{{p_t}}} $ (1)

式中：下标coarse表示粗化网格，而fine表示细化网格，p为对应不同网格的精度阶，v代表速度。而PIC方法是一种统计动理学方法，离散网格和无量纲形式h为

$ {h^p} = {\Bigg(\dfrac{{\Delta {x_{{\text{coarse}}}}}}{{\Delta {x_{{\text{fine}}}}}}\Bigg)^{{p_x}}} = {\Bigg(\dfrac{{\Delta {t_{{\text{coarse}}}}}}{{\Delta {t_{{\text{fine}}}}}}\Bigg)^{{p_t}}} = {\Bigg(\dfrac{{{w_{{\text{coarse}}}}}}{{{w_{{\text{fine}}}}}}\Bigg)^{{p_w}}} $ (2)

有时采用无量纲的网格步长更为方便

$ \left\{\begin{array}{l}\delta_x=\dfrac{\Delta x}{\lambda_{\text{D}}} \\ \delta_t=\dfrac{\lambda_{\text{D}}\omega_{\text{pe}}\Delta t}{\Delta x} \\ \delta_w=\dfrac{w}{n_{\text{ref}}\lambda_{\text{D}}^3}\end{array}\right. $ (3)

式中：$ {\lambda _{\text{D}}} $为德拜长度，$ {\omega _{{\text{pe}}}} $为等离子体振荡频率，$ {n_{{\text{ref}}}} $为参考数密度。相比直接动理学方法，PIC方法并没有直接在速度空间进行离散，而是在粒子空间进行离散。因此，为了获得光滑的粒子相图，需要减小宏粒子的权重以增加其数目。

第二，验证过程的观测量。PIC方法应该选取粒子量或者/并且网格量。粒子量是定义在拉格朗日网格的离散粒子分布函数$ f_{\mathrm{p}} $，它是Boltzmann方程的解。网格量包括定义在欧拉网格上的电势$ \varphi $、电场E、电荷密度ρ等，直接采用$ f_{\mathrm{p}} $难以与精确解进行对比。在一维情形下用$ f_{\mathrm{p}} $的积分即累积概率分布函数CDF来验证在理论上是可行的，但在2D2V或3D3V系统中定义和计算CDF存在困难。最近，Tranquili等人^[6]证明了网格量的误差范数与粒子分布函数的误差范数在同一量级。

$ \left\| f-\tilde{f} \right\| _2=\left\{\begin{array}{l}O( \left\| E-\tilde{E} \right\| _2) \\ O( \left\| \rho-\tilde{\rho} \right\| _{\infty})\end{array}\right. $ (4)

式中：下标代表误差范数的阶数。因此两种验证量的选择可以认为是等价的，而选择场量例如电场和电荷密度在实践中更加简单。

第三，广义网格的全局收敛精度。PIC方法中粒子和网格分别离散之后需要耦合，表现在电场插值与电荷分配的过程中。因此，对耦合方法的验证也是PIC方法验证中必要的部分。粒子模拟程序即使对场求解器采用相同的中心差分二阶精度收敛算法，对粒子时间积分算法也采用相同的Boris二阶精度收敛格式：粒子-网格耦合在采用二阶cloud in cell 格式时整个PIC程序是二阶精度收敛的；而采用一阶最近邻格点NGP格式则整个PIC程序是一阶精度收敛的^[4]。

对于有碰撞过程的放电等离子体，需要将PIC方法与蒙特卡罗（MC）方法进行混合。MC方法包括MCC方法和DSMC，前者在碰撞过程中入射物作为宏粒子追踪，而靶物质一般作为背景处理；后者不管是入射物还是靶粒子均应该作为宏粒子处理。MC方法在验证过程中有以下特点。

首先，将碰撞和运动过程进行解耦时通常忽略时间步长二阶以上的项，Boltzmann方程的时间离散写成

$ f(v, r, t+\Delta t)=(1+\Delta t \mathcal{C})(1+\Delta t \mathcal{F}) f(v, r, t) $ (5)

式中：$ \mathcal{F} $和$ \mathcal{C} $分别为表征运动和碰撞的算子。因此即使在运动过程采用了时间收敛二阶精度以上的求解格式，上述解耦过程会使得整个程序为时间收敛一阶精度^[16]。

其次，对于存在随机过程的模拟，同一程序采用不同的伪随机数发生器计算的结果会不一致^[17]；采用同一个随机数发生器，但是不同的随机数种子，结果也不一致^[18]。随机过程MC方法与基于确定论方法的本质不同，一般情形下会使得随机误差和离散误差混合在一起。采用系综平均的方法可以将随机误差和离散误差进行分离^[15]。

2 分层级验证方法

PIC-MC方法的验证需要考虑到第一部分中所描述的诸多特点，但直接对程序整体进行验证仍然存在一些不足：当模拟结果例如精度阶测试出现错误时难以快速定位是哪一个程序模块的问题；整体验证需要大量的计算资源；验证过程可能是代码侵入式的。由于PIC-MC程序通常是模块化的，我们建议采用分层级验证的策略。第一层级，分别对粒子推进、电磁场求解、粒子碰撞进行单元模块的收敛精度阶测试，这一层级能够识别代码模块中的微小缺陷，可以直接验证数学方程的求解过程，因此2.1节中的结果均是无量纲的；第二层级，将单元模块与粒子-场耦合，以及边界效应集成进行PIC、MCC、DSMC功能验证，在这一层级可以利用经典的物理问题构造验证算例，一般存在解析精确解；第三层级，采用基准算例进行PIC-MCC代码间比较，此时物理过程一般不存在解析解。下面将具体介绍静电PIC-DSMC程序的分层级验证过程。

2.1 单元模块验证

首先，粒子推进的验证通常以单个带电粒子在电磁场中的速度和位置与解析精确解进行比较。验证算例假设只存在径向电场$ {E_r} $和轴向磁场$ {B_{\textit{z}}} $，分别依赖于空间的径向和轴向位置^[19]，即

$ \left\{\begin{array}{*{20}{l}}E_r=E_{r_0}+kr， & E_{r_0}=0.01,\; k=0.01 \\ B_{{\textit{z}}}=B_{{\textit{z}}_0}+h{\textit{z}}， & B_{{\textit{z}}_0}=1.0,\; h=1.0\end{array}\right. $ (6)

采用的粒子的电荷数和质量都是1，其初始位置(r, θ, z)=(0.01, 0, 0)，初始速度为(v_r, v_θ, v_z)=(0.1, 0, 0)。根据牛顿第二定律，位置和速度各分量满足的常微分方程如下

$ \left\{\begin{gathered}\dfrac{\text{d}r}{\text{d}t}=v_r \\ \dfrac{\text{d}\theta}{\text{d}t}=\dfrac{v_{\theta}}{r} \\ \dfrac{\text{d}{\textit{z}}}{\text{d}t}=v_{{\textit{z}}} \\ \dfrac{\text{d}v_r}{\text{d}t}=\dfrac{q}{m}(E_r+v_{\theta}B_{{\textit{z}}}-v_{{\textit{z}}}B_{\theta})+\dfrac{v_{\theta}^2}{r} \\ \dfrac{\text{d}v_{\theta}}{\text{d}t}=\dfrac{q}{m}(E_{\theta}-v_rB_{{\textit{z}}}+v_{{\textit{z}}}B_r)-\dfrac{v_rv_{\theta}}{r} \\ \dfrac{\text{d}v_{{\textit{z}}}}{\text{d}t}=\dfrac{q}{m}(E_{{\textit{z}}}+v_rB_{\theta}-v_{\theta}B_r) \\ \end{gathered}\right. $ (7)

上述常微分方程的理论解可以通过Python中scipy库的odeint函数得到。图1（a）所示为时间步长取0.01时的粒子推进Boris算法和odeint的计算结果径向位置和径向速度的比较。为了定量描述两者的误差，图1（b）所示为全局误差随时间步长减小的趋势，虚线所示为二阶收敛的参考线，在达到收敛域之后，粒子推进模块的收敛速度与二阶收敛参考线的斜率一致。因此，粒子推进模块满足时间二阶精度收敛。

图 1. 自编程序的粒子推进Boris算法和odeint的计算结果比较

Fig. 1. Comparison of calculation results for self-implemented particle push using Boris algorithm and odeint

下载图片 查看所有图片

其次，对于二维柱坐标泊松方程

$ \nabla^2\phi=f，\quad\mathit{\Omega}=r\in[0,20]\times {\textit{z}}\in[0,20] $ (8)

人为构造以下精确解和源项

$ \left\{ \begin{gathered} \phi = \dfrac{{\textit{z}}({\textit{z}} - 20)({r^3}/3 - 10{r^2})}{10\;000} + 0.004\;8{\textit{z}} \\ f = \dfrac{{\textit{z}}({\textit{z}} - 20)(3r - 40)}{10\;000} + \dfrac{2({r^3}/3 - 10{r^2})}{10\;000} \\ \end{gathered} \right. $ (9)

满足边界条件

$ \left\{\begin{gathered}\phi|_{{\textit{z}}=0}=0 \\ \phi|_{{\textit{z}}=20}=0.096 \\ \dfrac{\partial\phi}{\partial r}\Bigg|_{r=0}=\dfrac{\partial\phi}{\partial r}\Bigg|_{r=20}=0 \\ \end{gathered}\right. $ (10)

静电泊松方程基于五点中心有限差分法进行离散，并通过直接LU分解求解矩阵方程，因此不存在迭代误差。数值模拟结果显示有限差分方法的计算结果收敛速率与空间二阶收敛的斜率一致。对于电磁PIC方法可以人为构造类似的解析解。另外，在求解过程中还需关注Maxwell方程组中的两个散度方程特别是磁场散度的相对误差。

$ {\varepsilon _{\nabla \cdot {\boldsymbol{B}}}} = \dfrac{{\displaystyle\int {{{\left| {\nabla \cdot {\boldsymbol{B}}} \right|}^2}{\text{d}}V} }}{{\displaystyle\int {{{\left| {\nabla \times {\boldsymbol{B}}} \right|}^2}{\text{d}}V} }} $ (11)

接着对于碰撞模块进行验证，气体之间的弹性碰撞是稀薄气体动力学中重要的碰撞过程之一。Maxwell气体的碰撞截面反比于相对速度，Bobylev-Krook-Wu （BKW）模型给出了Maxwell气体弛豫的一个特例^[20-21]，无量纲的速度分布函数f具有解析解，另外，在求解过程中还需关注Maxwell方程组中的两个散度方程特别是磁场散度的相对误差。

$ \left\{ \begin{gathered} f(u,\tau ) = \dfrac{1}{{{{(\pi C)}^{3/2}}}}\dfrac{1}{{2C}}\Bigg[5C - 3 + \dfrac{{2(1 - C){u^2}}}{C}\Bigg]\exp\Bigg( - \dfrac{{{u^2}}}{C}\Bigg) \\ C = 1 - \dfrac{2}{5}\exp \Bigg( - \dfrac{\tau }{6}\Bigg) \\ \tau = 4{\text{π}}nkt \\ u = \dfrac{v}{{\sqrt {{k_{\text{B}}}T/m} }} \\ \end{gathered} \right. $ (12)

通过分布函数对速度幂次的积分得到其偶数阶速度矩满足

$ {M_{2n}}(t) = {C^{n - 1}}[n - (n - 1)C] $ (13)

图2（a）所示为BKW弛豫下Maxwell气体的四阶和六阶速度矩MC模拟与解析精确解的比较，数据点是使用蒙特卡罗碰撞方法计算的，而曲线是理论值。类似的，对于不同组分A-B二体碰撞^[22]，BKW弛豫模型速度的四阶矩也存在解析精确解，图2（b）所示为MC模拟与解析精确解的比较。

图 2. BKW气体的速度矩与时间的关系：数据点是使用蒙特卡罗碰撞方法计算的，而曲线是理论值

Fig. 2. Velocity moments of BKW gas versus time: data points are calculated with Monte Carlo collision method and line is from theory

下载图片 查看所有图片

2.2 集成模块验证

集成模块功能验证如表1所示，分别包括DSMC、MCC和PIC验证。对于DSMC功能的验证将结合系综平均和理查德森外推详细介绍，评估随机误差和离散误差随时间步长的变化。对于MCC和PIC验证的过程类似，因此不再赘述，仅给出简单的算例。

DSMC功能的验证包含了中性粒子在无加速情形下的推进和粒子之间的直接碰撞。稀薄气体动力学领域的经典物理问题，包括温度差异导致的傅里叶流动，自由分子的Couette流动，自由边界激波结构，圆柱绕流等可以用来作为验证算例^[23]。以气体的傅里叶流动模型为例^[24]：两端壁面的温度恒定，壁面之间气体的初始温度介于两端的壁面温度大小之间，则稳恒条件下传导的热通量可以用傅里叶热传导公式计算。验证算例中，壁面两端距离为L=1×10⁻³ m，左边壁面温度T_c=223.15 K，右边壁面温度T_h=323.15 K，壁面之间初始均匀分布Ar气，气体温度273.15 K，压强266.64 Pa。气体碰撞采用硬球模型，硬球直径为3.66×10⁻¹⁰ m。气体输运至壁面的边界为完全反射热适应边界，即入射粒子被吸收，发射粒子的温度完全由壁面温度确定。

第一组模拟保持空间和时间步长不变，仅改变离散宏粒子的权重，因此宏粒子数目不同。由于DSMC模拟涉及随机过程，在宏粒子权重确定后又改变随机数发生器的种子，每组计算共变化32个种子。在系统模拟达到稳恒状态后，取1 ms的时间段对热通量进行时间平均：只有稳恒时间较长时抽样过程的统计误差可以忽略；如果抽样数过少，则此过程的统计误差可能与随机数种子变化产生的误差可比拟。对这32个不同随机数种子得到的模拟结果采用系综平均进行统计分析，随机统计误差用标准差来评价。图3所示为宏粒子权重不同时，DSMC模拟的热通量及其随机统计误差的变化。红色圆点所示为基于理查德森外推法通过两个粗化网格的解外推出细化网格的精确解的估值，与细化网格的直接DSMC模拟在误差范围内吻合。图中红色线性拟合直线证实了热通量满足宏粒子权重离散的一阶精度收敛。

图 3. 空间和时间步长固定时，热通量随宏粒子权重的变化

Fig. 3. Variation of heat flux with macro-particle weight when spatial and temporal steps are fixed

下载图片 查看所有图片

第二组模拟保持空间网格数和宏粒子权重不变，仅改变离散时间步长，保证随机误差小于离散误差。如图4（a）所示为热通量随时间步长的变化，图中黑色方框数据为DSMC模拟的结果，红色圆圈数据为通过两个粗化网格的Richardson外推得到的精确解估值，可以看到在时间步长小于1时与细化网格的DSMC模拟结果符合较好。图4（b）所示为离散误差和统计误差随时间步长的变化。可见，统计误差基本与时间步长无关而离散误差随时间步长满足一阶精度收敛，与公式（5）理论预计的一致。

图 4. 空间步长和宏粒子权重固定时，热通量、离散误差、统计误差随无量纲时间步长的变化

Fig. 4. Variation of heat flux, discretization error, and statistical error with dimensionless time step when spatial step and macro-particle weight are fixed

下载图片 查看所有图片

第三组模拟对空间步长、时间步长和宏粒子权重采用广义离散网格的一致细化，即空间步长和时间步长减小的同时，宏粒子的权重也变小。图5所示为不同细化网格下的DSMC模拟的热通量结果。图中黑色方框数据为DSMC模拟的结果，红色圆圈为基于理查德森外推的数据，而黑色虚线为根据傅里叶热传导公式给出的理论解。当时间步长和空间步长趋于零，和宏粒子数趋于无穷大时，通过数值解的理查德森外推方法得到的结果与理论精确解一致符合。

图 5. 离散网格一致细化时热通量随广义网格间距的变化

Fig. 5. Variation of heat flux with generalized grid spacing during grid refinement with discrete grid uniformity

下载图片 查看所有图片

程序的MCC功能采用带电粒子在等离子体中的输运特性来验证，可以同时考察带电粒子推进和测试粒子蒙特卡罗碰撞两个模块。上述物理过程中空间电荷效应不显著，因此场和粒子的耦合可以忽略。以电子的swarm模型为例，包含电子在弱电离等离子体中存在均匀电场情形下的运动以及与背景气体的碰撞过程。当背景气体为Maxwell气体时，电子平均能量存在解析精确解^[25]。电子与Maxwell气体之间仅有弹性碰撞，取碰撞截面为$ {\sigma}_{\text{el}}[\text{m}^2]={\sigma}_0/{g}= 6\times10^{-20}{\varepsilon}[\mathrm{eV}]^{-1/2} $，假设电场沿着y方向，则可以证明电子在给定电场下达到稳态时的能量为

$ \left\langle\boldsymbol{\mathit{\varepsilon}}\right\rangle=\boldsymbol{\varepsilon}_0+\dfrac{1}{2}\boldsymbol{\mathbf{\mathit{m}}}_{\text{gas}}\boldsymbol{\mathbf{\mathit{\mathit{v}}}}_y^2 $ (14)

$\mathit{\boldsymbol{\mathbf{\mathit{v}}}_y}=\dfrac{\boldsymbol{\mathit{e}\mathbf{\mathit{\mathit{E}}}}}{\boldsymbol{\mathbf{\mathit{m}}}_{\text{e}}\boldsymbol{\mathbf{\mathit{\nu}}}_{\text{coll}}}$ (15)

$\boldsymbol{\mathbf{\mathit{\nu}}}_{\text{coll}}=\mathit{\mathit{\boldsymbol{N}\boldsymbol{\mathbf{\mathit{v}}}}}_{\text{rel}}\boldsymbol{\mathbf{\mathit{\sigma}}}=\boldsymbol{\mathit{N}}(2/\mathbf{\boldsymbol{\mathit{\mu}}})^{1/2}\boldsymbol{\mathbf{\mathit{\sigma}}}_0 $ (16)

除此之外，还可采用硬球模型来验证弹性碰撞，在给定均匀场下电子的速度分布为Druyvesteyn分布。

真实气体相比Maxwell气体存在非弹性碰撞过程，因此电子swarm模型的输运系数不存在解析精确解。但是电子在弱电离等离子体中输运的数值解许多文献都有记载，已经形成了一套标准验证算例，包括Reid气体模型^[26]、Lucas-Salee气体模型^[27]，以及Ness-Robson气体模型^[28]，由此可以依次验证激发、电离、粘附等碰撞反应过程。表2还总结了离子在弱电离等离子体中的输运^[29]，以及完全电离等离子体中的输运问题^[30]。

PIC功能的验证主要包括以下模块：带电粒子推进、场的求解、场的插值与电荷电流的分配，以及需要的边界条件。等离子体物理中有许多经典问题存在（近似）解析解，可以用于PIC功能的验证，包括鞘层、空间电荷限制流、线性与非线性朗道阻尼、双流不稳定性、等离子体的自由膨胀^[31]、等离子体中的波^[32]等现象。以空间电荷限制流SCLC（space-charge limited current）为例，根据高压鞘层理论，一维平板之间由于空间电荷限制效应所能传输的最大电流密度$ {j_{{\text{SCLC}}}} $由极板之间压降$ V_{\text{g}} $和间距d决定。

$ j_{\text{SCLC}} = \dfrac{{4\varepsilon _{\text{0}}}}{9}{\left( {\dfrac{{2e}}{{m_{\text{e}}}}} \right)^{1/2}}\dfrac{{{V_{{\text{g}}}}^{3/2}}}{{{d^2}}} $ (17)

研究者发展了SCLC的模型，发现某些几何构型^[33-34]或电压波形^[35-36]下的引出电流密度比公式（17）的SCLC更强。例如两个球形电极之间所能传输的最大电流密度$ {j_{{\text{SCLC}}}} $除了决定于电极间距离和压降，还与阴极半径R_c、阳极半径R_a相关。我们采用一维球坐标PIC计算了SCLC并与Zhu等人^[34]的渡越时间理论进行了比较。图6（a）所示为阴极半径固定为1×10⁻³ m、阳极半径为2×10⁻³ m时，SCLC随间隙电压的变化，图6（b）所示为间隙电压固定时，SCLC随阴阳极半径比值的变化。在不同的参数变化情况下，PIC的计算结果均与理论解符合一致。

图 6. PIC计算的SCLC密度与朱氏模型的比较

Fig. 6. Comparison of SCLC density calculated by PIC and Zhu’s model

下载图片 查看所有图片

2.3 代码间比较

放电等离子体的数值模拟通常是PIC方法和MC方法的耦合。不同于前两节中经典的物理模型，难以获得解析精确解，或者通过人为方法构造精确解。因此代码间比较成为验证数值模拟程序的常用手段。我们将自研的PIC程序与某商业软件进行了对比。碰撞过程只考虑电子与假想原子的弹性和电离碰撞，采用二维轴对称平行平板电极结构，电极间距为3×10⁻⁴ m。两电极间为均匀分布的常温（293 K）背景气体，其压强为1×10⁴ Pa，初始存在均匀分布的等离子体，数密度为10¹³ m⁻³。离子轰击电极产生二次电子的产额设定为常数。采用与商业软件相同的内置参数和输入条件，输出不同时刻的宏粒子数，如图7所示。可以看到无论是电子和离子的宏粒子数目，两个程序的模拟结果是基本一致的。

图 7. 简化反应通道模型气体击穿的粒子模拟基准校验：与商业软件对比宏粒子数随时间的变化

Fig. 7. Benchmark validation of particle simulation for simplified reaction channel model of gas breakdown: comparison of macro-particle count vs time with commercial software

下载图片 查看所有图片

3 结　论

本文首先介绍了放电等离子体PIC模拟的验证方法区别于常见的计算流体力学验证的特点，主要体现在离散域、观测量和全局收敛精度。而MC方法中随机误差和离散误差的叠加也需要注意。随后，根据PIC-MC混合模拟的特点提出了一种分层级验证的思路，按照分层验证思想可以快速对任何PIC程序的子模块进行验证，也可作为PIC程序研制中的指导思想。第一层级验证对于粒子推进、场求解和粒子碰撞进行严格的收敛精度阶测试，从而快速定位代码中的微小缺陷；第二层级验证对于PIC、MCC和DSMC功能进行集成验证，基于经典物理问题例如空间电荷限制流、带电粒子在等离子体中的输运以及气体的傅里叶热流动等构造具有理论精确解的算例，从而快速检查模块与模块之间的耦合性和程序的准确性；第三层级对于不存在解析精确解的实际放电过程，设计简化的基准放电模型，采用代码间比较进行验证。本论文未涉及数值模拟的确认工作，后续需与实验方共同迭代讨论，给出适用于放电等离子体粒子模拟的确认方案。

致　谢　感谢哈尔滨工业大学（深圳）曹勇在静电场求解方面给予的指导和帮助。