非平衡态等离子体的仿真研究现状与新进展

荣命哲 刘定新 李 美 王伟宗

（西安交通大学电力设备电气绝缘国家重点实验室 西安 710049）

1 引言

等离子体是由大量带电粒子组成的非束缚态宏观体系[1]，因整体上呈电中性而得名。它是与固体、液体、气体3 种物质形态相并列的“物质第四态”。宇宙空间中99%的物质处于等离子体态，其温度、压力等物性参数范围跨越10 个数量级以上，相互之间具有巨大的特性差异。这种差异使得人们对等离子体往往进行分类研究，分类方式包括温度、压力、电离度、产生形式等。从仿真研究角度，热力学平衡条件是最常用的分类标准，据此可将等离子体分为3 类:

（1）热力学平衡态等离子体（Complete Thermal Equilibrium，CTE）:等离子体的宏观性质（温度、压力、体积、密度等）不随时间变化，但组成系统的大量微观粒子仍在不停运动，只是微观粒子所表现出的整体效果不变[2]。

（2）局部热力学平衡态等离子体（Local Thermal Equilibrium，LTE）:等离子体中碰撞过程（而非辐射过程）占主导，且碰撞过程与逆过程满足细致平衡条件；宏观性质在局部场中的梯度很小，使得物理运动特征时间不小于化学平衡特征时间，粒子在运动中有足够时间获得平衡[2-5]。

（3）非热力学平衡态等离子体（Non Local Thermal Equilibrium，NLTE）:等离子体的宏观性质随时间与空间变化，微观粒子温度差异大，其分布状态不满足Maxwell-Boltzman 分布规律。此类等离子体往往简称为非平衡态等离子体（下文同）[6]。

非平衡态等离子体技术在微电子制造、材料表面改性、镀膜等领域已经实现大规模产业化应用，产生了巨大的经济效益和社会效益[7]；在新兴的纳米技术、生物医学与环境保护等应用领域，也表现出广阔的应用前景[8,9]。非平衡态等离子体技术已成为当前国际高科技竞争的热点领域，深化相关理论与应用研究具有重要意义。

等离子体中弹性碰撞非常频繁，这是粒子间能量转移的主要过程。驱动等离子体的电场和磁场能量大部分被电子吸收，而电子能量主要通过弹性碰撞传递给重粒子。因为电子质量与重粒子质量之间的巨大差异，单次弹性碰撞转移的电子能量非常少。比如，在SF6等离子体中，SF6分子质量mh是电子质量me的约2.7×105倍，单次弹性碰撞转移的能量只有电子能量的2me/mh=1/1.8×106。因此，如果电子密度不够高，重粒子就难以通过足够的弹性碰撞获得与电子相当的能量，无法达到热力学平衡。一般说来，等离子体中电离度小于0.1%时，重粒子温度Th与电子温度Te之间就有明显差异，可认为等离子体处于非平衡状态[10-12]。

与弹性碰撞相对应，非弹性碰撞是等离子体中物质转化的关键过程。通过包括激发、电离、复合、电子吸附和解吸附等非弹性碰撞过程，等离子体中往往产生种类繁多的带电粒子、激发态粒子、自由基等。特别是多原子气体或混合气体形成的等离子体，其中粒子组分与化学过程更加复杂。比如:空气等离子体中常见的粒子种类有数百种，化学反应有数千个[13,14]。除了化学过程复杂外，这些粒子还受电场、磁场、温度场、气流场等多场耦合作用，物理过程也不易描述。

非平衡态等离子体中，如果粒子化学反应的特征时间比物理运动特征时间短，则认为等离子体处于局部化学平衡状态（Local Chemical Equilibrium，LCE）；反之则为非化学平衡状态（Non Local Chemical Equilibrium,NLCE）[15]。严格说来，LCE状态不易存在于非平衡态等离子体中，因为等离子体鞘层中物理场的梯度非常大，使得某些低能态粒子（如:振动激发态粒子）在存活时间内会产生大于德拜长度的位移。但是，针对特定的研究和应用需求，当这些粒子的作用不很重要时，可近似按LCE状态对等离子体进行仿真分析[16]。

在LCE 状态下，等离子体中的粒子组分及分布规律主要由非弹性碰撞过程决定。由于该状态下非弹性碰撞是平衡的可逆过程，所以各种粒子的密度是自身势能的函数，等离子体系统中所有粒子的吉布斯自由焓最小。基于这样的原理，在计算粒子密度与分布规律时，可以不考虑等离子体中复杂的微观物理与化学过程，而是通过宏观的热力学统计规律来得到。在NLCE 状态下，由于物理运动（漂移、对流与扩散）带来的影响不能忽略，统计物理的方法无法准确描述等离子体中粒子的密度与分布规律，深入考察等离子体中微观的物理与化学过程在所难免。具体说来，每一种粒子都需要一个质量守恒方程来描述，方程中既包括复杂的化学过程，又包括对流与扩散等物理过程。由于等离子体中粒子种类繁多，化学过程复杂，且这些过程的特征时间尺度差异巨大，所以相比而言，NLCE 状态等离子体的建模与计算难度远大于LCE 状态等离子体[15]。

过去的数十年内，人们对非平衡态等离子体开展了大量仿真与实验研究，但主要针对低气压等离子体[17]。在大气压及更高气压条件下，等离子体中的物理与化学过程极为复杂，研究工作相对滞后。本文主要针对高气压、非平衡态等离子体，对近十年来的仿真模型进行综述，分别介绍LCE 和NLCE状态下的建模与仿真方法，并重点介绍仿真研究的新进展。

2 局部化学平衡态等离子体的仿真模型

如引言中所述，在满足LCE 条件的非平衡等离子体中，粒子密度与分布规律可以通过宏观的热力学统计规律计算得到，无需深入考察等离子体中微观的物理和化学过程。研究表明，电子和重粒子的弹性碰撞使得电子的能量分布函数成平衡态的麦克斯韦分布，而非弹性碰撞和较高的电场则造成一个非平衡态的能量分布。一般认为，在等离子体的电离度小于0.1%后，电子和重粒子之间的弹性碰撞频率较低，能量交换不足，电子温度Te就会明显高于重粒子的温度Th。此时，等离子体可以用“双温度模型”描述，也就是仿真模型中有Te和Th两个温度变量。严格说来，等离子体中不同种类重粒子（离子、振动激发态粒子、基态粒子等）之间也有温度差异，在Th<9 000K 时需采用“多温度模型”[15]。但是温度变量越多模型越复杂，使得目前多温度模型应用很少，且通过双温度模型很容易推导得出多温度模型的控制方程。因此，本文重点介绍双温度模型。

双温度模型是在考察Te和Th两个温度变量条件下，采用热力学统计规律或化学动力学方法计算等离子体中任意位置的粒子组分，采用Navier-Stokes方程（简称N-S 方程）等流体方程计算等离子体中气体密度、温度与流速等物理特性，从而获取NLTE、LCE 态等离子体中微观粒子组分和宏观物理特性的时空演变规律。

2.1 双温度条件下的粒子组分计算模型

在双温度条件下计算粒子组分主要基于质量作用定律、Gibbs 自由焓最小化原理或化学动力学方法。其中，基于质量作用定律和Gibbs 自由焓最小化的两种方法是通过热力学统计规律计算LCE 态等离子体中粒子组分[18-21]，而化学动力学方法是基于反应动力学的规律得到粒子组分。

上述 3 种方法都需要与元素化学计量平衡条件、道尔顿分压定律、以及电荷准中性条件相结合，才能计算得出NLTE、LCE 态等离子体中任意位置的粒子组分[22]。

2.1.1 基于质量作用定律的计算模型

等离子体内部的电离和解离过程遵循质量作用定律原理，具体的表现形式为 Saha 电离方程和Guldberg-Waage 解离方程。应用质量作用定律计算粒子组分，不同学者推导出了各自不同的公式，形成了Potapov 方法[23-27]和Eindhoven 方法[23,28]等典型方法。在NLTE 条件下，这些方法都根据等离子体的双温度特点作了修正。下面将对较早出现的Potapov 方法作具体的阐述。

（1）Saha 电离方程。r+1 级电离方程式为:Ar↔Ar+1+e-EI,r+1。非平衡态时，电子温度Te和重粒子温度Th不再相等，用双温度模型来考虑，两者的比值Te/Th=θ。此时，配分函数Zr不再是电子和重粒子统一温度的函数，而是电子温度Te和重粒子温度Th的函数。

式中，ne、nr、nr+1分别是电子、粒子Ar、粒子Ar+1的数密度；Zr、Zr+1分别是粒子Ar、Ar+1的配分函数，电子的配分函数值是2；k、h、Te、Tex、me分别是玻耳兹曼常数、普朗克常量、电子温度、粒子激发温度以及电子的质量；EI,r+1为上述r+1 级电离反应的电离能；ΔEI,r+1为等离子体内部因粒子间相互作用而导致的电离能跌落。

Tex为粒子激发温度，其值与决定输运过程的主要粒子碰撞方式有关。对于电离反应涉及到的原子和离子，粒子的激发对配分函数的计算贡献最大，而粒子激发温度在双温度条件下近似认为等于电子温度。式（2）中，λd为等离子体的德拜长度，计算式为

式中，0ε 为真空介电常数；zt为第t 种粒子的带电电荷；nt为第t 种离子的密度；υ 为等离子体中粒子种类数。

（2）Guldberg-Waage 解离方程。解离反应方程式为:AB↔A+B-Ed。对于双温情况，考虑到解离反应主要是由于重粒子相互碰撞发生的，推出Zr(Te,Th)≈Zr(Th)(i=A,B,AB)。因此，对于非平衡态等离子体中解离反应涉及的各种粒子，其数密度关系用Guldberg-Waage 解离方程表示为

式中，ni、im、Zi（i=A, B,AB）分别代表粒子的数密度、质量和配分函数；Ed为该解离反应的解离能。

相比于Potapov 方法，Eindhoven 方法的Saha方程少了指数项1/θ，相应配分函数的指数项θ 也去掉了。研究发现，在热力学平衡条件下，Potapov和Eindhoven 方法有着很好的一致性[29]。但在非热力学平衡态时，Eindhoven 方法的结果更精确。对于Eindhoven 方法下的Saha 方程和Guldberg-Waage方程，在文献[30]中有详细表述。

2.1.2 基于Gibbs 自由焓最小化的计算模型

系统在温度和压力不变情况下，对于各种可能的变动，平衡态的Gibbs 函数最小，这就是Gibbs自由焓最小化计算模型的理论依据[30-32]。

双温条件下的Gibbs 自由焓定义由André 等人给出，公式如下[31]

采用 Gibbs 自由焓最小法计算等离子体组分时，需预先求出每种粒子在压力p0下的熵、焓和定压热容，计算比较繁琐复杂，增加了工作量，因而这种方法不常被采用。

2.1.3 化学动力学模型动力学模型的主要理论依据是化学反应的Arrchenius 定理[31-38]，Arrchenius 定理反映了一定浓度下化学反应速率随温度的变化关系。

通常的化学反应表达式为

式中，Ai为参加化学反应的粒子i，υji为正向反应化学计量系数；υ 'ji为逆向反应化学计量系数；N 为粒子总数；L 为总的反应个数。当等离子体处于LCE状态时，粒子物理运动的影响可以忽略，化学反应满足可逆条件。如下，第j 个反应的速率方程可表示为

基于动力学的模型准确度最高，但该方法中正反应速率系数kf和逆反应速率系数kr很难获得。目前，已有很多学者通过大量实验获得了反应速率系数并做了较为完善的归纳[39,40]，然而这些数据都还局限在一定的温度范围之内，高温情况（如:Th>4000K）的数据还很缺乏。

2.2 双温度条件下的流体模型

基于粒子组分的计算模型，可以根据特定Te和Th组合求取等离子体的输运参数，这些参数是建立流体模型的基础。进一步，运用流体模型可以求取等离子体宏观物理特性（包括Te和Th）的时空分布。因此，粒子组分计算模型与流体模型相结合，才能全面求取等离子体中宏观的物理特性与微观的粒子组分，以及它们的时空演变规律。

在NLTE、LCE 条件下，需要对LTE 条件适用的流体模型（一般为磁流体动力学模型）做相应的修改。简单说来，质量守恒方程与动量守恒方程保持不变，但电子和重粒子分别采用各自的能量守恒方程。

在粒子组分已知的前提下，式（8）～式（12）的输运参数（通常指的是扩散系数，粘滞系数，热导率、电导率）可以由热力学表达式计算得到。等离子体的输运参数跟粒子质量、动量和能量的传输有关，而后者由分子的随机运动和碰撞完成，可用Boltzmann 方程来描述。该方程是一个复杂的多重积分微分方程，直接求解非常困难，可以基于局部热力学平衡态假设在流动平衡状态附近，对其作Chapman-Enskog 展开近似得到输运参数[41]。

虽然双温模型早在1981 年就由Chen 等人提出[34]，但是迄今的应用还偏少。一般认为电弧等离子体是LTE 状态，但是在冷空气注入[43]，等离子体边界[44-49]以及弧后初期阶段[50]，等离子体一般处于NLTE、LCE 状态。采用双温模型仿真电弧等离子体非常少，目前仅限于稳态情况下的氩气以及氧气电弧[51,52]。笔者采用双温模型仿真了SF6电弧等离子体，发现在强气吹条件下，电弧边缘和喷口下游有明显偏离热力学平衡的情况发生（如图1 所示），说明有必要采用双温模型进行仿真分析。

图1 600A 电流下喷口SF6电弧非平衡度θ=Te/Th分布Fig.1 Distribution of non-equilibrium degree θ=Te/ Thin the nozzle of SF6arc under 600A current

3 非化学平衡态等离子体的仿真模型

目前，对NLCE 态等离子体的仿真主要有流体模型[53-56]，化学动力学模型[57-59]、PIC 模型[60,61]等。在大气压甚至更高气压条件下，PIC 模型的计算量往往难以承受，因此应用非常少。本文主要讨论流体模型与动力学模型，以及这两种模型相结合的仿真方法。

相比于LCE 状态的等离子体，NLCE 状态等离子体的仿真要复杂许多，这主要是因为如下几个原因[15]:

（1）NLCE 状态等离子体中，粒子组分不满足Saha 方程和吉布斯最小自由焓原理，因此每一种粒子都必须采用一个守恒方程来描述。粒子守恒方程中不但包括等离子体中复杂的化学反应，而且包括漂移、扩散、对流等物理过程。等离子体中粒子种类往往很多，使得仿真模型需要对大量的守恒方程进行耦合求解，计算量非常大。

（2）与粒子生成/去除密切相关的化学反应时间尺度差异大，有的反应特征时间小于1ns，有的则在ms 量级甚至更长时间。为了能够对每一种粒子进行精确计算，仿真模型的时间步长必须小于化学反应的最短特征时间；但为了求取粒子平衡的稳态解，仿真模型的总体时间尺度必须大于化学反应的最长特征时间（或者物理运动过程的特征时间）。这进一步加大了仿真模型的计算量。

（3）模型中的输入参数，包括迁移率、扩散率、化学反应率系数等，都是电子温度或气体温度的函数，在仿真计算过程中不断变化。尤其是化学反应率系数的数据资料还很不全面，大量的化学反应率系数只能通过估计和推算得来[62]。对于电子碰撞反应率，还常常假定电子能量符合Maxwell-Boltzman分布来获取反应率系数[63,64]，这都增加了模型的复杂性和不确定性。

3.1 基于漂移-扩散方程的流体模型

在流体模型中，NLCE 态等离子体中的粒子守恒方程可以表示为

式中，Γi表示粒子i 的通量，Si表示因非弹性碰撞（化学反应）引起的粒子生成/去除源项。带电粒子的通量主要包括漂移项与扩散项，Γi=μni/E-D ∇ni；中性粒子的通量主要包括对流项和扩散项，表示为Γi=niv -D∇ ni。

这个方程一般需要联立电子能量方程、气体温度方程和麦克斯韦方程进行求解。电子能量守恒方程表示如下[65]

式中，右边第一项表示电子从电场中获得的功率，第二项表示非弹性碰撞能量损失，第三项表示弹性碰撞（动量转移）能量损失。Γε表示电子能量的通量，在假定电子能量满足Boltzman 分布条件下[66]，

气体温度方程表示为[67]

式（15）中表明，除了对流与传导（右边第一项）带来的温度变化外，速度与压力变化（右边第二项），化学反应带来的热焓变化（右边第三项）和离子焦耳热（右边第四项）对气体温度变化起着关键作用。

联立式（13）～式（15）及麦克斯韦方程，就可以求解 NLCE 态等离子体中粒子密度及分布规律。因为每一种粒子就需要一个质量守恒方程，因此模型计算量非常大。现有的模型大多数只针对原子气体（如:He 和Ar），或者原子气体中含极少量分子气体（如:He+N2）的等离子体进行仿真，以避免复杂的粒子成分和化学过程带来的计算量过大、程序收敛性差的问题[68-72]。同时，大多数模型仅限于1 维，只有少量2 维模型报道[73-75]。

上述模型对LCE 态等离子体同样适用，但是由于LCE 态等离子体的气体温度比较高（接近或高于1eV），上述方程中的参数难以全面获得，无法实现准确求解。比如:某些气体分子的振动温度与气体温度相当，等离子体中存在大量振动激发态粒子，而这些粒子的化学反应率系数在现有数据资料中严重不全。

3.2 化学动力学模型

相比于流体模型而言，化学动力学模型无需计算等离子体内部的物理场与粒子运动过程，因此计算量大大降低，程序易收敛。对比NLCE 和LCE 态等离子体，化学动力学模型在本质上是一致的，只是前者中化学反应不满足可逆条件。也就是，式（6）也是 NLCE 态等离子体化学动力学模型的基本方程，但是式（7）不适用。虽然本质上几乎没有差别，但是化学动力学模型在NLCE 态等离子体仿真中应用广泛，而在LCE 态等离子体仿真中应用较少。这主要是因为如下两个原因:

（1）由于实验水平的限制，化学反应率系数往往只能在4 000K 温度以下测量得到[76,77]，这一温度范围内等离子体一般处于NLCE 态。对于更高温度的LCE 态等离子体，往往需要通过数据拟合求取化学反应率系数。由于化学反应率系数往往是温度的指数函数，拟合得到的数据误差很大，最终使得仿真结果准确度不高[15]。

（2）对于流体模型而言，LCE 态等离子体比NLCE 态等离子体仿真难度小很多。前者可对成分复杂的等离子体进行3 维仿真，而后者一般对简单成分等离子体进行1 维仿真。因此，对LCE 态等离子体，人们更趋向于采用流体模型；而对于成分复杂的NLCE 态等离子体，化学动力学模型是当前主要的选择。

普通的化学动力学模型忽略空间变化，且不考虑等离子体边界的影响，因此计算结果误差较大，往往只能用于定性判断和趋势分析。对于成分复杂的NLCE 态等离子体，当前除了化学动力学模型外，其他模型的应用还有困难。因此，人们期望通过改进化学动力学模型，获得更加准确的计算结果。加州大学伯克利分校Lieberman 教授率先开发了全局模型（global model）[78]。全局模型是一种优化的化学动力学模型，它通过考察等离子体特性，预设粒子的时空分布状态，并计算边界条件的影响。

全局模型在低气压NLCE 态等离子体仿真中大量应用[79-81]。研究发现，低气压RF 放电等离子体中带电粒子分布状态与电负性（α）大小有关[82]:①α 很小时，一般中心区域是电负性的，呈抛物线形状，边界区域是电正性的；电子密度在中心区域均匀分布，在边界区域逐步减小。②随着α 值的增加，电正性边界逐步消失。③如果α 继续增加，中心区域离子分布逐步平滑。在电极边界，一般采用Bohm 判据和热通量分别估计带电粒子和中性粒子的通量。通过设定带电粒子分布状态并引入边界条件，可以更加准确的求取等离子体中粒子密度。

全局模型应用到高气压（如:大气压）等离子体仿真是近几年才开始的[63,57-59]。相比于低气压等离子体，高气压等离子体有很大差别，具体说来主要包括如下几点[83]:①粒子空间分布状态不同；②电子平均自由程往往小于电极间距，因此在电负性α 值很大的时候，离子的焦耳热比较重要；③等离子体内部碰撞频繁，正离子在鞘层中的速率远低于Bohm 速率，不能简单使用Bohm 速率估计正离子的电极通量；④高气压下三体反应比较重要；⑤某些存活时间较长的粒子（如:O3）会累积到很高的密度，并通过等离子体侧边界对流或扩散出去，因而等离子体侧边界损失不可忽略。考虑高气压条件下等离子体的特点，建立的全局模型具有较高的计算准确度。

对射流等离子体，近年来研究人员还开发了1维全局模型[84,85]，通过在射流横截面上假定粒子的分布状态，在射流方向上应用活塞流假设，求取射流方向上粒子密度分布规律，并揭示其内部机理。

3.3 全局模型与流体模型的联合仿真方法

对比全局模型和流体模型，前者计算量小，但准确度低，且不能求取粒子的空间分布状态；后者可用于求取粒子的时空演变规律，准确度高，但计算量大。全局模型往往用于仿真化学过程复杂的等离子体，而这类等离子体目前难以直接应用流体模型进行仿真。比如:He+O2+H2O 等离子体中常见粒子达数十种，常见化学反应有近千个[59]。复杂的粒子组分和化学反应相互耦合，且粒子生命周期与化学反应特征时间变化范围很大（快到ns 以下，慢到ms 以上），会给流体模型带来难以承受的计算量。但是，在众多粒子和化学反应中，对等离子体特性及应用效果影响显著的粒子和化学反应数量却要少很多。如果只针对关键的粒子和化学反应进行分析，计算量将大大减小，流体模型就能够得以应用，这将大幅度提升人们对复杂等离子体的认知水平。前提条件是，如何准确提取关键的粒子和化学反应。

一般说来，人们通过简单比较化学反应率系数来提取关键的粒子和化学反应。这样的方法对于化学过程简单的等离子体有效，而复杂等离子体中大量粒子与化学反应相互耦合，难以实现准确提取。全局模型可以对等离子体特性进行整体化描述，且由于考察了等离子体空间分布状态与边界条件，其计算结果相比普通的化学动力学模型更加准确。因此，采用全局模型提取关键的粒子和化学反应，是一种相对准确、可行的方法。比如:应用全局模型对He+H2O（1～30×10-6）等离子体进行分析，从46 种粒子和577 个化学反应中提取出了21 种重要粒子和28 个主要的化学反应[57]。将这些提取出来的粒子和化学反应建立了简化的全局模型，相比于原模型的计算误差不超过60%，但计算速度提高了2 个数量级以上。考虑到流体模型中粒子、化学反应、物理场与粒子运动的耦合更加复杂，提取关键粒子和化学反应，计算量减小幅度会更大。

全局模型和流体模型在特性上具有互补性，因此采用全局模型和流体模型的联合仿真方法，对深入揭示复杂等离子体中的微观机理非常有效[86]。具体说来，就是针对特定的研究和应用需求，应用全局模型全面考察等离子体中的理化过程，提取关键粒子和化学反应，进而建立流体模型对等离子体进行时空分辨的仿真研究。这种新的仿真方法已经在He+O2、He+H2O 大气压冷等离子体仿真中得以较好应用[53,55-58]。图2 给出了大气压He+O2等离子体的仿真结果，其中全局模型与流体模型得到的粒子密度都在同一数量级[86]。

图2 大气压He+O2等离子体的全局模型和流体模型仿真结果Fig.2 Simulation results of global model and fluid model for atmospheric-pressure He+O2plasma

4 结论

非平衡态等离子体是当前国际高科技竞争的热点领域，对其开展深入研究具有重要意义。由于等离子体自身的复杂性，目前的实验技术需要与仿真研究相互配合，才能深入揭示此类等离子体的微观机理，阐释其宏观特性变化规律，进而推动非平衡等离子体的产业化应用发展。

目前，大规模产业化应用的非平衡态等离子体主要是在低气压条件下产生的，因此仿真研究也主要集中于低气压等离子体。对大气压甚至更高气压的非平衡态等离子体，仿真研究的报道还偏少。为此，本文综述了近年来高气压非平衡态等离子体的仿真模型，并重点强调了仿真方法的新进展，以飨读者。

一般说来，热力学平衡特征时间长于化学平衡特征时间，因而非平衡等离子体中往往包括局部化学平衡（LCE）和非化学平衡（NLCE）两种状态。在LCE 状态下，等离子体中某一位置的粒子组分主要由该位置处的非弹性碰撞过程（化学反应）决定，输运过程对粒子组分的影响可忽略。因此，可通过宏观的热力学统计原理（主要是质量作用定律和Gibbs 自由焓最小化原理）计算等离子体中任一位置的粒子组分，再结合N-S 方程等流体方程获取等离子体特性的时空分布规律。通过化学动力学方法深入考察微观粒子的理化过程，也可以求取等离子体中的粒子组分，且原理上比热力学统计方法更加准确。但是，LCE 态等离子体的温度比较高（一般不低于数千开），化学过程往往非常复杂，化学反应率系数非常缺乏，使得化学动力学方法的复杂度高、准确度低。目前，对LCE 态等离子体中粒子组分的求取主要是采用热力学统计方法。

非平衡态等离子体的仿真模型必需考察粒子温度的差异。针对满足LCE 条件的非平衡态等离子体，可近似认为重粒子具有相似的温度且远低于电子温度，因而将电子温度Te和重粒子温度Th分别作为模型的关键参数，由此建立的模型称为“双温度模型”。双温度模型的相关理论体系已经建立，贯穿粒子组分、物性参数和等离子体特性时空分布规律整个链条的仿真计算。尽管如此，由于粒子组分的计算方法还存有争议，对电子和重粒子之间能量交换的机理还有待于进一步研究，目前该模型的应用还偏少，人们还是习惯假定局部热力学平衡对非平衡态等离子体进行仿真研究。

相比于LCE 态等离子体，NLCE 态等离子体的仿真计算复杂得多。具体说来，每一种粒子需要一个质量守恒方程，方程中需要同时考虑化学反应和物理运动对粒子密度的影响。等离子体中粒子成分和化学过程比较复杂，尤其是多原子气体或混合气体形成的等离子体，常见粒子种类可达上百种，化学反应可达数千个。这些粒子的生命周期和化学反应的特征时间尺度可跨越6 个数量级以上。因此，现有的模型大多数只针对原子气体（如:He 和Ar），或者原子气体中含极少量分子气体（如:He+N2）的NLCE 态等离子体。对于复杂的NLCE 态等离子体（如:空气），应用流体模型进行仿真计算量过大，程序难以收敛。

近几年来，全局模型与流体模型的联合仿真方法得以提出，并针对NLCE 态等离子体成功实践。全局模型是一种优化的化学动力学模型，可以在假定空间分布的前提下对等离子体进行较为准确的定量分析，计算量小。应用全局模型系统的提取出关键的粒子和化学反应，并用之建立流体模型，可以在保证较高准确度的基础上，大大降低流体模型的计算量、提高收敛性。全局模型与流体模型的联合仿真方法，突破了化学过程复杂的NLCE 态等离子体仿真瓶颈，对于新兴的等离子体纳米和医学应用（这些应用往往采用化学过程复杂的等离子体）具有重要价值。

[1]刘万东.等离子体物理导论[M].合肥:中国科学技术大学出版社,2002.

[2]Boulos M I,Fauchais P,Pfender E.Thermal plasmas-fundamental and applications[M].New York:Plenum Press,1994.

[3]Freton P,Gonzalez J J,Gleizes A.Comparison between a two-and a three-dimensional arc plasma configuration[J].Journal of Physics D:Applied Physics,2010,43(19):2442-2452.

[4]Trelles J P,Pfender E,Heberlein J V R.Modelling of the arc reattachment process in plasma torches[J].Journal of Physics D:Applied Physics,2007,40(18):5635-5648.

[5]Bernardi D,Colombo V,Ghedini E,et al.3-D turbulent modelling of an ICPT with detailed gas injection section[J].IEEE Transactions on Plasma Science,2005,33(2):426-427.

[6]Murphy A B.Thermal plasmas in gas mixtures[J].Journal of Physics D:Applied Physics,2001,34(20):151-157.

[7]Plasma science:Advancing knowledge in the national interest.[P/OL].http://sites.nationalacademies.org/DE PS/ DEPS_037579.

[8]Samukawa S,Hori M,Rauf S,et al.The 2012 plasma roadmap[J].Journal of Physics D:Applied Physics,2012,45(25):434-441.

[9]Kong M G,Kroesen G,Morfill G,et al.Plasma medicine:an introductory review[J].New Journal of Physics,2009,11(11):1367-2360.

[10]Murphy A B.Electron heating in the measurement of electron temperature by Thomson scattering:Are thermal plasmas thermal [J].Physical Review Letters,2002,89(4):1-4.

[11]Gregori G,Kortshagen U,Heberlein J,et al.Analysis of Thomson scattering results from an arc plasma jet[J].Physical Review E,2002,65(046411):1-8.

[12]Dzierzega K,Zawadzki W,Pokrzywka B,et al.Experimental investigations of plasma perturbation in Thomson scattering applied to thermal plasma diagnostics[J].Physical Review E:Statistical Nonlinear and Soft Matter Physics,2006,74(2):30-45.

[13]NIST chemical kinetics database[DB].http://kinetics.nist.gov/kinetics/index.jsp.

[14]Millar T J,Farquhar P R A,Willacy K.The UMIST database for astrochemistry 1995[J].Astronomy and Astrophysics Supplement Series,1997,121:139-185.

[15]Rat V,Murphy A B,Aubreton J,et al.Treatment of non-equilibrium phenomena in thermal plasma flows[J].Journal of Physics D:Applied Physics,2008,41(3):183-190.

[16]Yan J D,Fang M T C,Liu Q S.Dielectric breakdown of a residual SF6plasma at 3000K under diatomic equilibrium [J].IEEE Transactions on Plasma Science,1997,4(1):114-119.

[17]Fauchais P,BoulosI M I,Pfender E.Thermal plasmas fundamentals and applications[M].New York:Plenum,1994.

[18]高尚惠,孙煜.热力学与统计物理学[M].北京:高等教育出版社,1980.

[19]André P,Abbaoui M,Lefort A,et al.Numerical method and composition in multi-temperature plasmas:application to an Ar-H2mixture[J].Plasma Chemistry and Plasma Processing,1996,16(3):379-381.

[20]André P.Partition functions and concentrations in plasmas out of thermal equilibrium[J].IEEE Transactions on Plasma Science,1995,23(3):453-457.

[21]Coufal O.Composition and thermodynamic properties of thermal plasma up to 50kK[J].Journal of Physics D:Applied Physics,2007,40(11):3371-3385.

[22]王伟宗,荣命哲,吴翊,等.平衡态与非平衡态等离子体物性参数计算模型研究[J].高压电器,2010,46(7):41-45.Wang Weizong,Rong Mingzhe,Wu Yi,et al.Studies on properties calculation models of equilibrium and non-equilibrium plasma[J].High Voltage Apparatus,2010,46(7):41-45.

[23]Gleizes A,Chervy,Gonzale J J.Calculation of a two-temperature plasma composition:bases and application to SF6[J].Journal of Physics D:Applied Physics,1999,32(16):2060-2067.

[24]Wang W Z,Rong M Z,Yan J D.The reactive thermal conductivity of thermal equilibrium and non equilibrium plasmas application to nitrogen [J].IEEE Transactions on Plasma Science,2012,40(4):980-990.

[25]Wang W Z,Rong M Z,Yan J D.Thermophysical properties of nitrogen plasmas under thermal equilibrium and non-equilibrium conditions [J].Physics of Plasmas,2011,18(10):2-19.

[26]Ghorui S,Heberlein J V R,Pfender E.Thermodynamic and transport properties of two temperature nitrogen-oxygen plasma[J].Plasma Chemistry and Plasma Processing,2008,7(28):553-582.

[27]Ramshaw J D.Simple approximation for thermal diffusion in gas mixtures[J].Journal of Nonequilibrium Thermodynamics,1996,21(3):233-238.

[28]Tanaka Y,Yokomizu Y,Ishikawa M,et al.Particle composition of high-pressure SF6plasma with electron temperature greater than gas temperature [J].IEEE Transactions on Plasma Science,1997,25(5):991-995.

[29]Chervy B,Gleizes A,Razafinimanana M.Thermodynamic properties and transport coefficients in SF6-Cu mixtures at temperatures of 300～3 0000 K and pressures of 0.1～1 MPa[J].Journal of Physics D:Applied Physics,1994,27(6):1193-1206.

[30]臧春艳,何俊佳,程礼椿.平衡态和非平衡态等离子体的微观模型研究[J].高压电器,2005,41(6):416-419.Zang Chunyan,He Junjia,Cheng Lichun.Study on microcosmic models of equilibrium and nonequilibrium plasma[J].High Voltage Apparatus,2005,41(6):416-419.

[31]Coufal O.Composition and thermodynamic properties of thermal plasma up to 50kK[J].Journal of Physics D:Applied Physics,2007,40 (13):3371–3385.

[32]Tanaka Yasunori,Michishita T,Uesugi Y.Hydrodynamic chemical non-equilibrium model of a pulsed arc discharge in dry air at atmospheric pressure[J].Plasma Sources Science and Technology,2005,14(5):134–151.

[33]Fauchais P,Boulos M I,Pfender E.Thermal plasmas-fundamentals and applications[M].New York:Plenum Press,1994.

[34]Chen D M,Hsu K C,Pfender E.Two-temperature modeling of an arc plasma reactor[J].Plasma Chemistry and Plasma Processing,1981,4(3):295-314.

[35]Li H P,Heberlein J,Pfender E.Three-dimensional,nonequilibrium effects in a high-intensity blown arc[J].IEEE Transactions on Plasma Science,2005,33(2):402-405.

[36]Dinulescu H A,Pfender E.Analysis of the anode boundary layer of high intensity arcs[J].Journal of Physics D:Applied Physics,1980,51(6):3149-3157.

[37]Morrow R,Lowke J J.A one dimensional theory for the electric sheaths of electric arcs[J].Journal of Physics D:Applied Physics,1993,26(3):634-642.

[38]Amakawa T,Jenista J,Heberlein J,et al.Anode-boundary-layer behaviour in a transferred high-intensity arc[J].Journal of Physics D:Applied Physics,1998,31(20):2826-2834.

[39]Tanaka M,Ushio M.Observations of the anode boundary layer in free-burning argon arcs[J].Journal of Physics D:Applied Physics,1999,32(8):906-912.

[40]Benilov M S.Analysis of ionization non-equilibrium in the near-cathode region of atmospheric-pressure arcs[J].Journal of Physics D:Applied Physics,1999,32(3):257-262.

[41]Yang G,Heberlein J.Anode attachment modes and their formation in a high intensity argon arc [J].Plasma Sources Science and Technology,2007,16(3):529-542.

[42]Yang G,Heberlein J.The anode region of high intensity arcs with cold cross flow [J].Journal of Physics D:Applied Physics,2007,40(18):5649-5662.

[43]Colombo V,Ghedini E,M Boselli,et al.3D static and time-dependent modelling of a dc transferred arc twin torch system[J].Journal of Physics D:Applied Physics,2011,44(9):194-199.

[44]Trelles J P,Heberlein J V R,Pfender E.Non-equilibrium modelling of arc plasma torches[J].Journal of Physics D:Applied Physics,2007,40(6):5937–5952.

[45]Adamec L,Coufal O.On kinetics of reactions in HV circuit breakers after current zero[J].Journal of Physics D:Applied Physics,1999,32(6):1702–1710.

[46]Girard R,Belhaouari J B,Gonzalez J J,et al.A two-temperature kinetic model of SF6plasma[J].Journal of Physics D:Applied Physics,1999,32(5):2890–2901.

[47]Gleizes A,Mbolidi F,Habibt A A M.Kinetic model of a decaying SF6plasma over the temperature range 12000 K to 3000 K[J].Plasma Sources Science and Technology,1993,2 (10):173-179.

[48]Richley E,Tuma D T.On the determination of particle concentrations in multi-temperature plasmas[J].Journal of Physics D:Applied Physics,1982,53(12):8537-8542.

[49]Teulet P H,Cressault Y,Hingana H,et al.Calculation of two-temperature thermodynamic and transport properties for argon,oxygen,nitrogen and air plasmas at atmospheric pressure[C].International Symposium on Plasma Chemistry,2011.

[50]Hossain M M,Tanaka Y,Sakuta T.Particle concentrations and transport properties of a partially ionized Ar plasma in a two-temperature reaction kinetic approach[J].Journal of Physics D:Applied Physics,2002,35(7):529-535.

[51]Phelps A V,Van Brunt R.Electron-transport,ionization,attachment and dissociation coefficients in SF6and its mixtures[J].Journal of Physics D:Applied Physics,1988,64(9):4269-4277.

[52]Atsushi Komuro,Ryo Ono,Tetsuji Oda.Kinetic model of vibrational relaxation in a humid-air pulsed corona discharge[J].Plasma Sources Science and Technology,2010,19(5):055004.

[53]Liu D X,Yang A J,Wang X H,et al.Wall fluxes of reactive oxygen species of an rf atmospheric-pressure plasma and their dependence on sheath dynamics[J].Journal of Physics D:Applied Physics,2012,45(30):305205.

[54]Mckay K,Liu D X,Rong M Z,et al.Generation and loss of reactive oxygen species in low-temperature atmospheric-pressure RF He+O2+H2O plasmas[J].Journal of Physics D:Applied Physics,2012,45(17):172001.

[55]Yang A J,Wang X H,Rong M Z,et al.1-D fluid model of atmospheric-pressure rf He+O2cold plasmas:Parametric study and critical evaluation[J].Journal of Physics D:Applied Physics,2011,18(11):113503.

[56]McKay K,Liu D X,Iza F,et al.Double-layer structures in low-temperature atmospheric-pressure electronegative RF microplasmas:separation of electrons and anions[J].IEEE Transactions on Plasma Science,2011,39(11):2138-2139.

[57]Liu D X,Bruggeman P,Iza F,et al.Global model of low-temperature atmospheric-pressure He+H2O plasmas[J].Plasma Sources Science and Technology,2010,19(2):025018.

[58]Liu D X,Rong M Z,Wang X H,et al.Main species and physicochemical processes in cold atmospheric pressure He+O2plasmas[J].Plasma Processes and Polymers,2010,7(9):846-865.

[59]Liu D X,Iza F,Wang X H,et al.He+O2+H2O plasmas as a source of reactive oxygen species[J].Applied Physics Letters,2011,98(22):221501.

[60]Kim H C,Iza F,Yang S S,et al.Paricle and fluid simulations of low-temperature plasma discharges:benchmarks and kinetic effects[J].Journal of Physics D:Applied Physics.2005,38(10):283-301.

[61]Iza F,Lee J K,Kong M G.Electron Kinetics in radio-frequency atmospheric-pressure microplasmas[J].Physical Review Letters,2007,99(7):075004.

[62]Kossyi I A,Kostinsky A Y,Matveyev A A,et al.Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures[J].Plasma Sources Science and Technology,1992,1(5):207-220.

[63]Park G,Lee H,Kim G,et al.Global model of He/O2and Ar/O2atmospheric pressure glow discharges[J].Plasma Processes and Polymers,2008,5(3):569-576.

[64]Rauf S,Kushner M J.Dynamics of a coplanarelectrode plasma display panel cell[J].Journal of Physics D:Applied Physics,1999,85(14):3460-3469.

[65]Sakiyama Y,Graves D B,Stoffels E.Influence of electrical properties of treated surface on RF-excited plasma needle at atmospheric pressure[J].Journal of Physics D:Applied Physics,2008,41(9):095204.

[66]Hagelaar G J M,Pitchford L C.Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[J].Plasma Sources Science and Technology,2005,14(7):722-733.

[67]Kang W S,Kim H S,Hong S H.Gas temperature effect on discharge-mode characteristics of atmosphericpressure dielectric barrier discharge in a helium-oxygen mixture[J].IEEE Transactions on Plasma Science,2010,38(8):1982-1990.

[68]Shi J J,Kong M G.Mechanisms of the α and γ modes in radio-frequency atmospheric glow discharges [J].Journal of Physics D:Applied Physics,2005,97(2):023306.

[69]Martens T,Bogaerts A,Brok W,et al.Computer simulations of a dielectric barrier discharge used for analytical spectrometry[J].Analytical and Bioanalytical Chemistry,2007,388(45):1583-1594.

[70]Balcon N,Hagelaar G J M,Boeuf J P.Numerical model of an argon atmospheric pressure RF discharge[J].IEEE Transactions on Plasma Science,2008,36(5):2782-2787.

[71]张增辉,邵先军,张冠军.大气压氩气介质阻挡辉光放电的一维仿真研究[J].物理学报,2012,61(4):452-464.Zhang Zenghui,Shao Xianjun,Zhang Guanjun.One-dimensional simulation of argon dielectric barrier glow discharge at atmospheric pressure[J].Journal of Physics,2012,61(4):452-464.

[72]张远涛.大气压介质阻挡放电时空演化行为理论研究[D].大连:大连理工大学,2006.

[73]Sakiyama Y,Graves D B.Neutral gas flow and ring-shaped emission profile in non-thermal RF-excited plasma needle discharge at atmospheric pressure[J].Plasma Sources Science and Technology,2009,18(5):342-347.

[74]Al-Mamun S A,Tanaka Y,Uesugi Y.Twotemperature two-dimensional non chemical modeling of Ar-CO2-H2induction thermal plasmas at atmospheric pressure[J].Plasma Chem.Plasma Process,2010,30(4):141-172.

[75]Trelles J P,Heberlein J V R,Pfender E.Nonequilibrium modelling of arc plasma torches[J].Journal of Physics D:Applied Physics,2007,40(3):5937-5952.

[76]Gordillo-V´azquez F J.Air plasma kinetics under the influence of sprites[J].Journal of Physics D:Applied Physics,2008,41(23):3-32.

[77]Baulch D L,Cox R A,Hampson R F.Evaluated kinetic and photochemical data for atmospheric chemistry:Supplement II[J].Journal of Physical and Chemical Reference Data,1984,13(4):1259-1378.

[78]Lieberman M A,Lichtenberg A J.Principles of plasma discharges and materials processing[M].2nd ed.NewYork:Wiley,2005.

[79]Monahan D D,Turner M M.Global models of electronegative discharges:critical evaluation and practical recommendations[J].Plasma Sources Science and Technology,2008,17(3):45-53.

[80]Gudmundsson J T.On the effect of the electron energy distribution on the plasma parameters of an argon discharge:a global (volume-averaged) model study[J].Plasma Sources Science and Technology,2001,10(5):76-81.

[81]Lee C, Lieberman M A.Global model of Ar,O2,Cl2and Ar/O2high-density plasma discharges[J].Journal of Vacuum Science and Technology A Science,1995,13(2):368-380.

[82]Kim S,Lieberman M A,Lichtenberg A J.Improved volume-averaged model for steady and pulsed-power electronegative discharges[J].Journal of Vacuum Science and Tecfunology A,2006,24(6):2025-2039.

[83]刘定新.大气压复杂冷等离子体全局模型的建模与仿真研究[D].西安:西安交通大学,2010.

[84]Gonzalez E,Barankin M D,Guschl P C,et al.Surface activation of poly(methyl methacrylate) via remote atmospheric pressure plasma[J].Plasma Processes and Polymers,2010,7(13):482-493.

[85]Stafford D S,Kushner M J.O2(1Δ) production in He/O2mixtures in flowing low pressure plasmas [J].Journal of Physics D:Applied Physics,2004,96(5):2451-2465.

[86]Li M,Liu D X,Yang A J,et al.Synergy of global model and fluid model for numerical study on plasmas with complex chemistry[C].Proceedings of the 19th International Conference on Gas Discharges and Their Applications,2012:844-847.