榮命哲 劉定新 李 美 王偉宗

（西安交通大學電力設備電氣絕緣國家重點實驗室 西安 710049）

1 引言

等離子體是由大量帶電粒子組成的非束縛態(tài)宏觀體系[1]，因整體上呈電中性而得名。它是與固體、液體、氣體3 種物質形態(tài)相并列的“物質第四態(tài)”。宇宙空間中99%的物質處于等離子體態(tài)，其溫度、壓力等物性參數(shù)范圍跨越10 個數(shù)量級以上，相互之間具有巨大的特性差異。這種差異使得人們對等離子體往往進行分類研究，分類方式包括溫度、壓力、電離度、產生形式等。從仿真研究角度，熱力學平衡條件是最常用的分類標準，據(jù)此可將等離子體分為3 類:

（1）熱力學平衡態(tài)等離子體（Complete Thermal Equilibrium，CTE）:等離子體的宏觀性質（溫度、壓力、體積、密度等）不隨時間變化，但組成系統(tǒng)的大量微觀粒子仍在不停運動，只是微觀粒子所表現(xiàn)出的整體效果不變[2]。

（2）局部熱力學平衡態(tài)等離子體（Local Thermal Equilibrium，LTE）:等離子體中碰撞過程（而非輻射過程）占主導，且碰撞過程與逆過程滿足細致平衡條件；宏觀性質在局部場中的梯度很小，使得物理運動特征時間不小于化學平衡特征時間，粒子在運動中有足夠時間獲得平衡[2-5]。

（3）非熱力學平衡態(tài)等離子體（Non Local Thermal Equilibrium，NLTE）:等離子體的宏觀性質隨時間與空間變化，微觀粒子溫度差異大，其分布狀態(tài)不滿足Maxwell-Boltzman 分布規(guī)律。此類等離子體往往簡稱為非平衡態(tài)等離子體（下文同）[6]。

非平衡態(tài)等離子體技術在微電子制造、材料表面改性、鍍膜等領域已經(jīng)實現(xiàn)大規(guī)模產業(yè)化應用，產生了巨大的經(jīng)濟效益和社會效益[7]；在新興的納米技術、生物醫(yī)學與環(huán)境保護等應用領域，也表現(xiàn)出廣闊的應用前景[8,9]。非平衡態(tài)等離子體技術已成為當前國際高科技競爭的熱點領域，深化相關理論與應用研究具有重要意義。

等離子體中彈性碰撞非常頻繁，這是粒子間能量轉移的主要過程。驅動等離子體的電場和磁場能量大部分被電子吸收，而電子能量主要通過彈性碰撞傳遞給重粒子。因為電子質量與重粒子質量之間的巨大差異，單次彈性碰撞轉移的電子能量非常少。比如，在SF6等離子體中，SF6分子質量mh是電子質量me的約2.7×105倍，單次彈性碰撞轉移的能量只有電子能量的2me/mh=1/1.8×106。因此，如果電子密度不夠高，重粒子就難以通過足夠的彈性碰撞獲得與電子相當?shù)哪芰?，無法達到熱力學平衡。一般說來，等離子體中電離度小于0.1%時，重粒子溫度Th與電子溫度Te之間就有明顯差異，可認為等離子體處于非平衡狀態(tài)[10-12]。

與彈性碰撞相對應，非彈性碰撞是等離子體中物質轉化的關鍵過程。通過包括激發(fā)、電離、復合、電子吸附和解吸附等非彈性碰撞過程，等離子體中往往產生種類繁多的帶電粒子、激發(fā)態(tài)粒子、自由基等。特別是多原子氣體或混合氣體形成的等離子體，其中粒子組分與化學過程更加復雜。比如:空氣等離子體中常見的粒子種類有數(shù)百種，化學反應有數(shù)千個[13,14]。除了化學過程復雜外，這些粒子還受電場、磁場、溫度場、氣流場等多場耦合作用，物理過程也不易描述。

非平衡態(tài)等離子體中，如果粒子化學反應的特征時間比物理運動特征時間短，則認為等離子體處于局部化學平衡狀態(tài)（Local Chemical Equilibrium，LCE）；反之則為非化學平衡狀態(tài)（Non Local Chemical Equilibrium,NLCE）[15]。嚴格說來，LCE狀態(tài)不易存在于非平衡態(tài)等離子體中，因為等離子體鞘層中物理場的梯度非常大，使得某些低能態(tài)粒子（如:振動激發(fā)態(tài)粒子）在存活時間內會產生大于德拜長度的位移。但是，針對特定的研究和應用需求，當這些粒子的作用不很重要時，可近似按LCE狀態(tài)對等離子體進行仿真分析[16]。

在LCE 狀態(tài)下，等離子體中的粒子組分及分布規(guī)律主要由非彈性碰撞過程決定。由于該狀態(tài)下非彈性碰撞是平衡的可逆過程，所以各種粒子的密度是自身勢能的函數(shù)，等離子體系統(tǒng)中所有粒子的吉布斯自由焓最小?；谶@樣的原理，在計算粒子密度與分布規(guī)律時，可以不考慮等離子體中復雜的微觀物理與化學過程，而是通過宏觀的熱力學統(tǒng)計規(guī)律來得到。在NLCE 狀態(tài)下，由于物理運動（漂移、對流與擴散）帶來的影響不能忽略，統(tǒng)計物理的方法無法準確描述等離子體中粒子的密度與分布規(guī)律，深入考察等離子體中微觀的物理與化學過程在所難免。具體說來，每一種粒子都需要一個質量守恒方程來描述，方程中既包括復雜的化學過程，又包括對流與擴散等物理過程。由于等離子體中粒子種類繁多，化學過程復雜，且這些過程的特征時間尺度差異巨大，所以相比而言，NLCE 狀態(tài)等離子體的建模與計算難度遠大于LCE 狀態(tài)等離子體[15]。

過去的數(shù)十年內，人們對非平衡態(tài)等離子體開展了大量仿真與實驗研究，但主要針對低氣壓等離子體[17]。在大氣壓及更高氣壓條件下，等離子體中的物理與化學過程極為復雜，研究工作相對滯后。本文主要針對高氣壓、非平衡態(tài)等離子體，對近十年來的仿真模型進行綜述，分別介紹LCE 和NLCE狀態(tài)下的建模與仿真方法，并重點介紹仿真研究的新進展。

2 局部化學平衡態(tài)等離子體的仿真模型

如引言中所述，在滿足LCE 條件的非平衡等離子體中，粒子密度與分布規(guī)律可以通過宏觀的熱力學統(tǒng)計規(guī)律計算得到，無需深入考察等離子體中微觀的物理和化學過程。研究表明，電子和重粒子的彈性碰撞使得電子的能量分布函數(shù)成平衡態(tài)的麥克斯韋分布，而非彈性碰撞和較高的電場則造成一個非平衡態(tài)的能量分布。一般認為，在等離子體的電離度小于0.1%后，電子和重粒子之間的彈性碰撞頻率較低，能量交換不足，電子溫度Te就會明顯高于重粒子的溫度Th。此時，等離子體可以用“雙溫度模型”描述，也就是仿真模型中有Te和Th兩個溫度變量。嚴格說來，等離子體中不同種類重粒子（離子、振動激發(fā)態(tài)粒子、基態(tài)粒子等）之間也有溫度差異，在Th<9 000K 時需采用“多溫度模型”[15]。但是溫度變量越多模型越復雜，使得目前多溫度模型應用很少，且通過雙溫度模型很容易推導得出多溫度模型的控制方程。因此，本文重點介紹雙溫度模型。

雙溫度模型是在考察Te和Th兩個溫度變量條件下，采用熱力學統(tǒng)計規(guī)律或化學動力學方法計算等離子體中任意位置的粒子組分，采用Navier-Stokes方程（簡稱N-S 方程）等流體方程計算等離子體中氣體密度、溫度與流速等物理特性，從而獲取NLTE、LCE 態(tài)等離子體中微觀粒子組分和宏觀物理特性的時空演變規(guī)律。

2.1 雙溫度條件下的粒子組分計算模型

在雙溫度條件下計算粒子組分主要基于質量作用定律、Gibbs 自由焓最小化原理或化學動力學方法。其中，基于質量作用定律和Gibbs 自由焓最小化的兩種方法是通過熱力學統(tǒng)計規(guī)律計算LCE 態(tài)等離子體中粒子組分[18-21]，而化學動力學方法是基于反應動力學的規(guī)律得到粒子組分。

上述 3 種方法都需要與元素化學計量平衡條件、道爾頓分壓定律、以及電荷準中性條件相結合，才能計算得出NLTE、LCE 態(tài)等離子體中任意位置的粒子組分[22]。

2.1.1 基于質量作用定律的計算模型

等離子體內部的電離和解離過程遵循質量作用定律原理，具體的表現(xiàn)形式為 Saha 電離方程和Guldberg-Waage 解離方程。應用質量作用定律計算粒子組分，不同學者推導出了各自不同的公式，形成了Potapov 方法[23-27]和Eindhoven 方法[23,28]等典型方法。在NLTE 條件下，這些方法都根據(jù)等離子體的雙溫度特點作了修正。下面將對較早出現(xiàn)的Potapov 方法作具體的闡述。

（1）Saha 電離方程。r+1 級電離方程式為:Ar?Ar+1+e-EI,r+1。非平衡態(tài)時，電子溫度Te和重粒子溫度Th不再相等，用雙溫度模型來考慮，兩者的比值Te/Th=θ。此時，配分函數(shù)Zr不再是電子和重粒子統(tǒng)一溫度的函數(shù)，而是電子溫度Te和重粒子溫度Th的函數(shù)。

式中，ne、nr、nr+1分別是電子、粒子Ar、粒子Ar+1的數(shù)密度；Zr、Zr+1分別是粒子Ar、Ar+1的配分函數(shù)，電子的配分函數(shù)值是2；k、h、Te、Tex、me分別是玻耳茲曼常數(shù)、普朗克常量、電子溫度、粒子激發(fā)溫度以及電子的質量；EI,r+1為上述r+1 級電離反應的電離能；ΔEI,r+1為等離子體內部因粒子間相互作用而導致的電離能跌落。

Tex為粒子激發(fā)溫度，其值與決定輸運過程的主要粒子碰撞方式有關。對于電離反應涉及到的原子和離子，粒子的激發(fā)對配分函數(shù)的計算貢獻最大，而粒子激發(fā)溫度在雙溫度條件下近似認為等于電子溫度。式（2）中，λd為等離子體的德拜長度，計算式為

式中，0ε 為真空介電常數(shù)；zt為第t 種粒子的帶電電荷；nt為第t 種離子的密度；υ 為等離子體中粒子種類數(shù)。

（2）Guldberg-Waage 解離方程。解離反應方程式為:AB?A+B-Ed。對于雙溫情況，考慮到解離反應主要是由于重粒子相互碰撞發(fā)生的，推出Zr(Te,Th)≈Zr(Th)(i=A,B,AB)。因此，對于非平衡態(tài)等離子體中解離反應涉及的各種粒子，其數(shù)密度關系用Guldberg-Waage 解離方程表示為

式中，ni、im、Zi（i=A, B,AB）分別代表粒子的數(shù)密度、質量和配分函數(shù)；Ed為該解離反應的解離能。

相比于Potapov 方法，Eindhoven 方法的Saha方程少了指數(shù)項1/θ，相應配分函數(shù)的指數(shù)項θ 也去掉了。研究發(fā)現(xiàn)，在熱力學平衡條件下，Potapov和Eindhoven 方法有著很好的一致性[29]。但在非熱力學平衡態(tài)時，Eindhoven 方法的結果更精確。對于Eindhoven 方法下的Saha 方程和Guldberg-Waage方程，在文獻[30]中有詳細表述。

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

系統(tǒng)在溫度和壓力不變情況下，對于各種可能的變動，平衡態(tài)的Gibbs 函數(shù)最小，這就是Gibbs自由焓最小化計算模型的理論依據(jù)[30-32]。

雙溫條件下的Gibbs 自由焓定義由André 等人給出，公式如下[31]

采用 Gibbs 自由焓最小法計算等離子體組分時，需預先求出每種粒子在壓力p0下的熵、焓和定壓熱容，計算比較繁瑣復雜，增加了工作量，因而這種方法不常被采用。

2.1.3 化學動力學模型動力學模型的主要理論依據(jù)是化學反應的Arrchenius 定理[31-38]，Arrchenius 定理反映了一定濃度下化學反應速率隨溫度的變化關系。

通常的化學反應表達式為

式中，Ai為參加化學反應的粒子i，υji為正向反應化學計量系數(shù)；υ 'ji為逆向反應化學計量系數(shù)；N 為粒子總數(shù)；L 為總的反應個數(shù)。當?shù)入x子體處于LCE狀態(tài)時，粒子物理運動的影響可以忽略，化學反應滿足可逆條件。如下，第j 個反應的速率方程可表示為

基于動力學的模型準確度最高，但該方法中正反應速率系數(shù)kf和逆反應速率系數(shù)kr很難獲得。目前，已有很多學者通過大量實驗獲得了反應速率系數(shù)并做了較為完善的歸納[39,40]，然而這些數(shù)據(jù)都還局限在一定的溫度范圍之內，高溫情況（如:Th>4000K）的數(shù)據(jù)還很缺乏。

2.2 雙溫度條件下的流體模型

基于粒子組分的計算模型，可以根據(jù)特定Te和Th組合求取等離子體的輸運參數(shù)，這些參數(shù)是建立流體模型的基礎。進一步，運用流體模型可以求取等離子體宏觀物理特性（包括Te和Th）的時空分布。因此，粒子組分計算模型與流體模型相結合，才能全面求取等離子體中宏觀的物理特性與微觀的粒子組分，以及它們的時空演變規(guī)律。

在NLTE、LCE 條件下，需要對LTE 條件適用的流體模型（一般為磁流體動力學模型）做相應的修改。簡單說來，質量守恒方程與動量守恒方程保持不變，但電子和重粒子分別采用各自的能量守恒方程。

在粒子組分已知的前提下，式（8）～式（12）的輸運參數(shù)（通常指的是擴散系數(shù)，粘滯系數(shù)，熱導率、電導率）可以由熱力學表達式計算得到。等離子體的輸運參數(shù)跟粒子質量、動量和能量的傳輸有關，而后者由分子的隨機運動和碰撞完成，可用Boltzmann 方程來描述。該方程是一個復雜的多重積分微分方程，直接求解非常困難，可以基于局部熱力學平衡態(tài)假設在流動平衡狀態(tài)附近，對其作Chapman-Enskog 展開近似得到輸運參數(shù)[41]。

雖然雙溫模型早在1981 年就由Chen 等人提出[34]，但是迄今的應用還偏少。一般認為電弧等離子體是LTE 狀態(tài)，但是在冷空氣注入[43]，等離子體邊界[44-49]以及弧后初期階段[50]，等離子體一般處于NLTE、LCE 狀態(tài)。采用雙溫模型仿真電弧等離子體非常少，目前僅限于穩(wěn)態(tài)情況下的氬氣以及氧氣電弧[51,52]。筆者采用雙溫模型仿真了SF6電弧等離子體，發(fā)現(xiàn)在強氣吹條件下，電弧邊緣和噴口下游有明顯偏離熱力學平衡的情況發(fā)生（如圖1 所示），說明有必要采用雙溫模型進行仿真分析。

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

3 非化學平衡態(tài)等離子體的仿真模型

目前，對NLCE 態(tài)等離子體的仿真主要有流體模型[53-56]，化學動力學模型[57-59]、PIC 模型[60,61]等。在大氣壓甚至更高氣壓條件下，PIC 模型的計算量往往難以承受，因此應用非常少。本文主要討論流體模型與動力學模型，以及這兩種模型相結合的仿真方法。

相比于LCE 狀態(tài)的等離子體，NLCE 狀態(tài)等離子體的仿真要復雜許多，這主要是因為如下幾個原因[15]:

（1）NLCE 狀態(tài)等離子體中，粒子組分不滿足Saha 方程和吉布斯最小自由焓原理，因此每一種粒子都必須采用一個守恒方程來描述。粒子守恒方程中不但包括等離子體中復雜的化學反應，而且包括漂移、擴散、對流等物理過程。等離子體中粒子種類往往很多，使得仿真模型需要對大量的守恒方程進行耦合求解，計算量非常大。

（2）與粒子生成/去除密切相關的化學反應時間尺度差異大，有的反應特征時間小于1ns，有的則在ms 量級甚至更長時間。為了能夠對每一種粒子進行精確計算，仿真模型的時間步長必須小于化學反應的最短特征時間；但為了求取粒子平衡的穩(wěn)態(tài)解，仿真模型的總體時間尺度必須大于化學反應的最長特征時間（或者物理運動過程的特征時間）。這進一步加大了仿真模型的計算量。

（3）模型中的輸入?yún)?shù)，包括遷移率、擴散率、化學反應率系數(shù)等，都是電子溫度或氣體溫度的函數(shù)，在仿真計算過程中不斷變化。尤其是化學反應率系數(shù)的數(shù)據(jù)資料還很不全面，大量的化學反應率系數(shù)只能通過估計和推算得來[62]。對于電子碰撞反應率，還常常假定電子能量符合Maxwell-Boltzman分布來獲取反應率系數(shù)[63,64]，這都增加了模型的復雜性和不確定性。

3.1 基于漂移-擴散方程的流體模型

在流體模型中，NLCE 態(tài)等離子體中的粒子守恒方程可以表示為

式中，Γi表示粒子i 的通量，Si表示因非彈性碰撞（化學反應）引起的粒子生成/去除源項。帶電粒子的通量主要包括漂移項與擴散項，Γi=μni/E-D ?ni；中性粒子的通量主要包括對流項和擴散項，表示為Γi=niv -D? ni。

這個方程一般需要聯(lián)立電子能量方程、氣體溫度方程和麥克斯韋方程進行求解。電子能量守恒方程表示如下[65]

式中，右邊第一項表示電子從電場中獲得的功率，第二項表示非彈性碰撞能量損失，第三項表示彈性碰撞（動量轉移）能量損失。Γε表示電子能量的通量，在假定電子能量滿足Boltzman 分布條件下[66]，

氣體溫度方程表示為[67]

式（15）中表明，除了對流與傳導（右邊第一項）帶來的溫度變化外，速度與壓力變化（右邊第二項），化學反應帶來的熱焓變化（右邊第三項）和離子焦耳熱（右邊第四項）對氣體溫度變化起著關鍵作用。

聯(lián)立式（13）～式（15）及麥克斯韋方程，就可以求解 NLCE 態(tài)等離子體中粒子密度及分布規(guī)律。因為每一種粒子就需要一個質量守恒方程，因此模型計算量非常大?，F(xiàn)有的模型大多數(shù)只針對原子氣體（如:He 和Ar），或者原子氣體中含極少量分子氣體（如:He+N2）的等離子體進行仿真，以避免復雜的粒子成分和化學過程帶來的計算量過大、程序收斂性差的問題[68-72]。同時，大多數(shù)模型僅限于1 維，只有少量2 維模型報道[73-75]。

上述模型對LCE 態(tài)等離子體同樣適用，但是由于LCE 態(tài)等離子體的氣體溫度比較高（接近或高于1eV），上述方程中的參數(shù)難以全面獲得，無法實現(xiàn)準確求解。比如:某些氣體分子的振動溫度與氣體溫度相當，等離子體中存在大量振動激發(fā)態(tài)粒子，而這些粒子的化學反應率系數(shù)在現(xiàn)有數(shù)據(jù)資料中嚴重不全。

3.2 化學動力學模型

相比于流體模型而言，化學動力學模型無需計算等離子體內部的物理場與粒子運動過程，因此計算量大大降低，程序易收斂。對比NLCE 和LCE 態(tài)等離子體，化學動力學模型在本質上是一致的，只是前者中化學反應不滿足可逆條件。也就是，式（6）也是 NLCE 態(tài)等離子體化學動力學模型的基本方程，但是式（7）不適用。雖然本質上幾乎沒有差別，但是化學動力學模型在NLCE 態(tài)等離子體仿真中應用廣泛，而在LCE 態(tài)等離子體仿真中應用較少。這主要是因為如下兩個原因:

（1）由于實驗水平的限制，化學反應率系數(shù)往往只能在4 000K 溫度以下測量得到[76,77]，這一溫度范圍內等離子體一般處于NLCE 態(tài)。對于更高溫度的LCE 態(tài)等離子體，往往需要通過數(shù)據(jù)擬合求取化學反應率系數(shù)。由于化學反應率系數(shù)往往是溫度的指數(shù)函數(shù)，擬合得到的數(shù)據(jù)誤差很大，最終使得仿真結果準確度不高[15]。

（2）對于流體模型而言，LCE 態(tài)等離子體比NLCE 態(tài)等離子體仿真難度小很多。前者可對成分復雜的等離子體進行3 維仿真，而后者一般對簡單成分等離子體進行1 維仿真。因此，對LCE 態(tài)等離子體，人們更趨向于采用流體模型；而對于成分復雜的NLCE 態(tài)等離子體，化學動力學模型是當前主要的選擇。

普通的化學動力學模型忽略空間變化，且不考慮等離子體邊界的影響，因此計算結果誤差較大，往往只能用于定性判斷和趨勢分析。對于成分復雜的NLCE 態(tài)等離子體，當前除了化學動力學模型外，其他模型的應用還有困難。因此，人們期望通過改進化學動力學模型，獲得更加準確的計算結果。加州大學伯克利分校Lieberman 教授率先開發(fā)了全局模型（global model）[78]。全局模型是一種優(yōu)化的化學動力學模型，它通過考察等離子體特性，預設粒子的時空分布狀態(tài)，并計算邊界條件的影響。

全局模型在低氣壓NLCE 態(tài)等離子體仿真中大量應用[79-81]。研究發(fā)現(xiàn)，低氣壓RF 放電等離子體中帶電粒子分布狀態(tài)與電負性（α）大小有關[82]:①α 很小時，一般中心區(qū)域是電負性的，呈拋物線形狀，邊界區(qū)域是電正性的；電子密度在中心區(qū)域均勻分布，在邊界區(qū)域逐步減小。②隨著α 值的增加，電正性邊界逐步消失。③如果α 繼續(xù)增加，中心區(qū)域離子分布逐步平滑。在電極邊界，一般采用Bohm 判據(jù)和熱通量分別估計帶電粒子和中性粒子的通量。通過設定帶電粒子分布狀態(tài)并引入邊界條件，可以更加準確的求取等離子體中粒子密度。

全局模型應用到高氣壓（如:大氣壓）等離子體仿真是近幾年才開始的[63,57-59]。相比于低氣壓等離子體，高氣壓等離子體有很大差別，具體說來主要包括如下幾點[83]:①粒子空間分布狀態(tài)不同；②電子平均自由程往往小于電極間距，因此在電負性α 值很大的時候，離子的焦耳熱比較重要；③等離子體內部碰撞頻繁，正離子在鞘層中的速率遠低于Bohm 速率，不能簡單使用Bohm 速率估計正離子的電極通量；④高氣壓下三體反應比較重要；⑤某些存活時間較長的粒子（如:O3）會累積到很高的密度，并通過等離子體側邊界對流或擴散出去，因而等離子體側邊界損失不可忽略?？紤]高氣壓條件下等離子體的特點，建立的全局模型具有較高的計算準確度。

對射流等離子體，近年來研究人員還開發(fā)了1維全局模型[84,85]，通過在射流橫截面上假定粒子的分布狀態(tài)，在射流方向上應用活塞流假設，求取射流方向上粒子密度分布規(guī)律，并揭示其內部機理。

3.3 全局模型與流體模型的聯(lián)合仿真方法

對比全局模型和流體模型，前者計算量小，但準確度低，且不能求取粒子的空間分布狀態(tài)；后者可用于求取粒子的時空演變規(guī)律，準確度高，但計算量大。全局模型往往用于仿真化學過程復雜的等離子體，而這類等離子體目前難以直接應用流體模型進行仿真。比如:He+O2+H2O 等離子體中常見粒子達數(shù)十種，常見化學反應有近千個[59]。復雜的粒子組分和化學反應相互耦合，且粒子生命周期與化學反應特征時間變化范圍很大（快到ns 以下，慢到ms 以上），會給流體模型帶來難以承受的計算量。但是，在眾多粒子和化學反應中，對等離子體特性及應用效果影響顯著的粒子和化學反應數(shù)量卻要少很多。如果只針對關鍵的粒子和化學反應進行分析，計算量將大大減小，流體模型就能夠得以應用，這將大幅度提升人們對復雜等離子體的認知水平。前提條件是，如何準確提取關鍵的粒子和化學反應。

一般說來，人們通過簡單比較化學反應率系數(shù)來提取關鍵的粒子和化學反應。這樣的方法對于化學過程簡單的等離子體有效，而復雜等離子體中大量粒子與化學反應相互耦合，難以實現(xiàn)準確提取。全局模型可以對等離子體特性進行整體化描述，且由于考察了等離子體空間分布狀態(tài)與邊界條件，其計算結果相比普通的化學動力學模型更加準確。因此，采用全局模型提取關鍵的粒子和化學反應，是一種相對準確、可行的方法。比如:應用全局模型對He+H2O（1～30×10-6）等離子體進行分析，從46 種粒子和577 個化學反應中提取出了21 種重要粒子和28 個主要的化學反應[57]。將這些提取出來的粒子和化學反應建立了簡化的全局模型，相比于原模型的計算誤差不超過60%，但計算速度提高了2 個數(shù)量級以上?？紤]到流體模型中粒子、化學反應、物理場與粒子運動的耦合更加復雜，提取關鍵粒子和化學反應，計算量減小幅度會更大。

全局模型和流體模型在特性上具有互補性，因此采用全局模型和流體模型的聯(lián)合仿真方法，對深入揭示復雜等離子體中的微觀機理非常有效[86]。具體說來，就是針對特定的研究和應用需求，應用全局模型全面考察等離子體中的理化過程，提取關鍵粒子和化學反應，進而建立流體模型對等離子體進行時空分辨的仿真研究。這種新的仿真方法已經(jīng)在He+O2、He+H2O 大氣壓冷等離子體仿真中得以較好應用[53,55-58]。圖2 給出了大氣壓He+O2等離子體的仿真結果，其中全局模型與流體模型得到的粒子密度都在同一數(shù)量級[86]。

圖2 大氣壓He+O2等離子體的全局模型和流體模型仿真結果Fig.2 Simulation results of global model and fluid model for atmospheric-pressure He+O2plasma

4 結論

非平衡態(tài)等離子體是當前國際高科技競爭的熱點領域，對其開展深入研究具有重要意義。由于等離子體自身的復雜性，目前的實驗技術需要與仿真研究相互配合，才能深入揭示此類等離子體的微觀機理，闡釋其宏觀特性變化規(guī)律，進而推動非平衡等離子體的產業(yè)化應用發(fā)展。

目前，大規(guī)模產業(yè)化應用的非平衡態(tài)等離子體主要是在低氣壓條件下產生的，因此仿真研究也主要集中于低氣壓等離子體。對大氣壓甚至更高氣壓的非平衡態(tài)等離子體，仿真研究的報道還偏少。為此，本文綜述了近年來高氣壓非平衡態(tài)等離子體的仿真模型，并重點強調了仿真方法的新進展，以饗讀者。

一般說來，熱力學平衡特征時間長于化學平衡特征時間，因而非平衡等離子體中往往包括局部化學平衡（LCE）和非化學平衡（NLCE）兩種狀態(tài)。在LCE 狀態(tài)下，等離子體中某一位置的粒子組分主要由該位置處的非彈性碰撞過程（化學反應）決定，輸運過程對粒子組分的影響可忽略。因此，可通過宏觀的熱力學統(tǒng)計原理（主要是質量作用定律和Gibbs 自由焓最小化原理）計算等離子體中任一位置的粒子組分，再結合N-S 方程等流體方程獲取等離子體特性的時空分布規(guī)律。通過化學動力學方法深入考察微觀粒子的理化過程，也可以求取等離子體中的粒子組分，且原理上比熱力學統(tǒng)計方法更加準確。但是，LCE 態(tài)等離子體的溫度比較高（一般不低于數(shù)千開），化學過程往往非常復雜，化學反應率系數(shù)非常缺乏，使得化學動力學方法的復雜度高、準確度低。目前，對LCE 態(tài)等離子體中粒子組分的求取主要是采用熱力學統(tǒng)計方法。

非平衡態(tài)等離子體的仿真模型必需考察粒子溫度的差異。針對滿足LCE 條件的非平衡態(tài)等離子體，可近似認為重粒子具有相似的溫度且遠低于電子溫度，因而將電子溫度Te和重粒子溫度Th分別作為模型的關鍵參數(shù)，由此建立的模型稱為“雙溫度模型”。雙溫度模型的相關理論體系已經(jīng)建立，貫穿粒子組分、物性參數(shù)和等離子體特性時空分布規(guī)律整個鏈條的仿真計算。盡管如此，由于粒子組分的計算方法還存有爭議，對電子和重粒子之間能量交換的機理還有待于進一步研究，目前該模型的應用還偏少，人們還是習慣假定局部熱力學平衡對非平衡態(tài)等離子體進行仿真研究。

相比于LCE 態(tài)等離子體，NLCE 態(tài)等離子體的仿真計算復雜得多。具體說來，每一種粒子需要一個質量守恒方程，方程中需要同時考慮化學反應和物理運動對粒子密度的影響。等離子體中粒子成分和化學過程比較復雜，尤其是多原子氣體或混合氣體形成的等離子體，常見粒子種類可達上百種，化學反應可達數(shù)千個。這些粒子的生命周期和化學反應的特征時間尺度可跨越6 個數(shù)量級以上。因此，現(xiàn)有的模型大多數(shù)只針對原子氣體（如:He 和Ar），或者原子氣體中含極少量分子氣體（如:He+N2）的NLCE 態(tài)等離子體。對于復雜的NLCE 態(tài)等離子體（如:空氣），應用流體模型進行仿真計算量過大，程序難以收斂。

近幾年來，全局模型與流體模型的聯(lián)合仿真方法得以提出，并針對NLCE 態(tài)等離子體成功實踐。全局模型是一種優(yōu)化的化學動力學模型，可以在假定空間分布的前提下對等離子體進行較為準確的定量分析，計算量小。應用全局模型系統(tǒng)的提取出關鍵的粒子和化學反應，并用之建立流體模型，可以在保證較高準確度的基礎上，大大降低流體模型的計算量、提高收斂性。全局模型與流體模型的聯(lián)合仿真方法，突破了化學過程復雜的NLCE 態(tài)等離子體仿真瓶頸，對于新興的等離子體納米和醫(yī)學應用（這些應用往往采用化學過程復雜的等離子體）具有重要價值。

[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]高尚惠,孫煜.熱力學與統(tǒng)計物理學[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]王偉宗,榮命哲,吳翊,等.平衡態(tài)與非平衡態(tài)等離子體物性參數(shù)計算模型研究[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]臧春艷,何俊佳,程禮椿.平衡態(tài)和非平衡態(tài)等離子體的微觀模型研究[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.