有源等離子體在開放大氣環境下的反應擴散過程研究*

尹增謙 趙盼盼 董麗芳 房同珍

1)(華北電力大學數理系，保定 071003)

2)(河北大學物理科學與技術學院，保定 071002)

3)(中國科學院物理研究所，北京 100190)

(2010年1月29日收到;2010年4月26日收到修改稿)

有源等離子體在開放大氣環境下的反應擴散過程研究*

尹增謙1)?趙盼盼2)董麗芳2)房同珍3)

1)(華北電力大學數理系，保定 071003)

2)(河北大學物理科學與技術學院，保定 071002)

3)(中國科學院物理研究所，北京 100190)

(2010年1月29日收到;2010年4月26日收到修改稿)

利用一維模型用數值模擬的方法，研究了低氣壓開放環境下大氣等離子體在存在等離子體源的情況下的反應擴散過程.得到了考慮化學反應、擴散以及漂移共同作用下的大氣等離子體的主要成分隨等離子體注入流量的變化規律及一定流量情況下主要帶電成分隨時空的演化規律.將數值模擬結果與一個近似解析公式相銜接，估計達到穩態時維持一定電子密度所需要的等離子體流量，可據此進一步估計所需功率.

大氣等離子體，等離子體源，數值模擬，反應擴散過程

PACS:52.50.Dg，92.60.H － ，52.65.－y，92.70.Cp

1.引 言

由于大氣等離子體的產生不需要真空設備，成本低廉，而且具有極其廣泛的應用前景，例如廣泛應用于化工、材料、環境保護、醫療和等離子體隱身等領域［1—4］，近年來，大氣等離子體演化過程成為研究熱點.其中，在空間飛行器表面產生一層等離子體以改變其電磁或動力性質，是一項正在發展的新技術.在開放空間中，能否在所需要的時間段內保持一定空間分布的等離子體，使其達到所需要的顯著效應，是此項新技術成敗的關鍵.考慮到快的擴散、流動和復合過程，須有強的等離子體源補充這些損失，因此載體必須提供充分的功率.產生和維持一定濃度和厚度的等離子體所需要的等離子體流量和所需源的功率是必須研究的課題.

在大氣等離子體過程的早期研究中，Vidmar［5］編制了大氣化學程序用于計算大氣中等離子體的產生和維持過程，著重研究了等離子體壽命與高度及等離子體密度的關系，只考慮三體復合，在存在種子氣體的條件下，如果電子壽命為10－6s，為維持1012cm－3的電子密度，需要的功率為1 W/cm3.這一結果顯然偏小，因為他只考慮了反應過程而未考慮擴散.Ouyang等［6］也研究了只考慮化學反應的零維大氣等離子體的演化規律.Fang等［7］研究了零維情況下重復脈沖放電產生等離子體的演化規律，并進一步編制了一維程序，考慮到擴散過程以研究這種放電模式的等離子體衰減規律［8］.后一工作發展了偏微分方程的 QSSA(quasi-steady-state approximation)技術，所得主要結論是基于地面大氣環境，電子損失的主要途徑是原地復合;而在高空，擴散是主要的.這說明，零維計算是遠遠不夠的.

從實際應用出發，最重要的是回答維持一定密度和厚度的等離子體需要提供多大的等離子體流.知道了這個答案，才能解決在飛行器表面產生等離子體層的可能性問題.從物理上說，須在一定等離子體流量下，在一個半無限空間解一維擴散及不同組分間的反應速率方程.此處我們暫且不計相對于空間的流動，這樣得到的結果可給出需要等離子體流量的下限.

由于這一方程組是高度非線性的(一般考慮兩體及三體反應過程)，只能用數值解.然而，數值解難于在半無限空間進行，往往選擇在一定厚度內計算.這時，在選擇的邊界處使用人為的邊界條件使計算能進行下去.如果研究暫態過程［8］，這樣的邊界條件對結果影響不大.但是這樣的邊界處理會產生不穩定性，可能影響穩態的結果.

本工作首先仍在一維有限空間數值模擬有等離子體源的反應擴散方程，使用人為的邊界條件.計算進行到充分長的時間后，利用計算結果取得的數據代入一個近似解析解以求得達到穩態時的等離子體密度和厚度，以此反推達到一定密度和厚度所需的等離子體流.采用這一方案還可節省計算時間.

2.反應擴散模型和數值計算方法

2.1.反應擴散模型

對于發生器發出的實際等離子體，粒子的空間演化過程除了純擴散以外，還有遷移和反應過程的影響.為研究純擴散、遷移、反應共同作用下的粒子數密度的時空演化規律，首先分析純擴散和遷移的影響.在一維情況下，粒子的通量(定義為單位時間內通過單位面積的粒子數)為

式中Γk為第k種粒子的通量，μk，Dk分別為第 k種粒子的遷移率和擴散系數，nk為第k種粒子的粒子數密度，E為電場強度.對于右邊第一項前面的±號，正離子取“+”，負離子取“－”，對于中性粒子，由于沒有電場作用下的遷移，μk=0.

為了便于數值計算，把(1)式改寫為等效的擴散形式.對大氣等離子體做簡化:1)電中性假設，即n+=n－+ne，其中 ne為電子密度;2)通量相等假設，即Γ+=Γ－+Γe;3)正離子、負離子、中性粒子的擴散系數相等，正離子、負離子的遷移率相等，即μi=μ－=μ+，Di=D－=D+.上述各式中的 n+，Γ+分別為總的正離子密度和通量，n－，Γ－分別為總的負離子密度和通量，μ+，D+分別為正離子的遷移率和擴散系數，μ－，D－分別為負離子的遷移率和擴散系數，Di為中性粒子的擴散系數.

根據以上三個假設，得到電場強度分布以及電子、負離子、各種正離子、中性粒子的通量分別為

式中 Γe，Γ－m，Γ+m，Γm分別為電子、負離子、正離子和中性粒子的通量.再考慮到反應項的貢獻，得到普遍的一維大氣等離子體反應擴散模型的基本公式

其中，Pk為單位時間、單位體積內第 k種成分的反應生成量，Lknk為單位時間、單位體積內第k種成分的反應消耗量.顯然，反應項 Pk－Lknk由與該成分有反應關系的其他成分的濃度和相應的反應速率系數決定，使得等離子體中各種成分時空演化的計算變得復雜，必須利用數值模擬的方法來計算.另一方面，如(3)式所示，利用愛因斯坦關系

其中kB為波爾茲曼常數，得

這樣只需要溫度和擴散系數就可以了.

本文的模型中考慮了34種反應物(見表1)和145個反應式，其反應速率系數由文獻［9］中查得.

表1 主要反應物列表

2.2.泛定方程的離散化方法

普遍的一維大氣等離子體反應擴散模型的基本公式(4)可以化簡為如下的標準泛定方程:

其中u為粒子數密度，D為等效的擴散系數，f相應于反應項Pk－Lknk.采用時間、空間的中心差分格式［10］，并做代換得到恒穩的實用差分格式采用這種三層的菱形格式，在數值計算時必須處理好初始條件和邊界條件.

2.3.初始條件和邊界條件的處理

由于菱形的差分格式不是自啟動的，所以在數值計算時，假設初始時刻與前一時間步長的時刻各種成分的密度空間分布對應相等.而對于邊界條件，假設左邊界點上的各成分密度與左邊空間格點上的密度相同.注入項的處理是必須要注意的一個問題，在本工作中，在每一個時間步長中使電子，在最左邊的第一個格點增加一個定值Q，電子的相對比例為100%，21%，79%，以保證電中性，即的比例按照 O2，N2在大氣當中的比例產生.當時間步長 t0改變時，最左邊的第一個格點增加的定值Q也隨之改變，以保證乘積Qt0(t0為時間步長)為定值.盡管右邊界為無限空間邊界，但是在數值模擬時只能計算有限的空間范圍，所以對于右邊界上的格點采用插值外推的方法得到各成分密度，即

2.4.QSSA方法的應用

鑒于不同化學反應速率系數有數量級的差異，造成方程組的剛性，為節省計算時間，在本工作的數值計算中，仍然采用QSSA方法［6］，即根據時間步長Δt與成分τk的壽命之比，采用不同的計算方法，以提高計算效率.

3.計算結果分析

3.1.模擬條件

模型中假設等離子體由氣體放電產生，帶電粒子主要是電子(e)和兩種正離子，在放電過程中電子以一定的速率產生.為保持電中性條件，正負電荷相同，兩種正離子按照大氣中的比例產生.

選取60 km高空作為模擬環境，其壓強為20 Pa，大氣溫度為253 K;電子溫度為2 eV，約合23200 K，離子溫度為500 K，電子擴散系數為1.76×107cm2/s，等離子體中各種離子擴散系數相差不大，為300 cm2/s［11］.假設初始時刻模擬空間中不存在等離子體，各粒子的初始粒子數密度與本底大氣相同［11］.空間步長選為 0.1 cm，時間步長為 10－12s，選取51個格點進行模擬.所以模擬區域的邊界為x=l=5 cm.

3.2.模擬結果與分析

圖1為源點處電子數密度在不同流量(實線為流量 1018cm－2·s－1，短劃線為流量 1016cm－2·s－1，虛線為流量 1014cm－2·s－1)下隨時間的演化，為了便于觀察，圖中采用雙對數坐標.如圖1所示，模擬初期，電子數密度都迅速增長.隨著時間的推移，電子數密度先后出現了增長減緩的趨勢，流量越高減緩趨勢出現得越早.尤其是流量為 1018cm－2·s－1時，電子數密度在10－5s左右達到一個極大值6×1012cm－3，隨后衰減，在模擬時間內已經趨于穩定.這種現象是反應過程與空間過程共同作用的結果.流量越大，趨于穩定的時間越短，是由于在60 km高空，粒子數密度越大，電子的壽命越低［5］，其衰減就越快與注入量達到平衡.在等離子體低目標技術中，為了滿足低目標技術要求，等離子體的臨界密度至少要在很短的時間內達到1010cm－3以上.下面著重對等離子體流量為1018cm－2·s－1的情況進行研究.

圖1 不同流量下源點處電子數密度隨時間的演化 實線、短劃線和虛線分別代表流量為 1018，1016和 1014cm －2·s－1的情況

圖2為零維和一維等離子體源點處電子數密度隨時間的演化圖.一維等離子體過程中電子在源點的流量設為 1018cm－2·s－1，為了便于比較，我們將零維過程中電子的注入量設為 1019cm－3·s－1，即兩個過程在單位時間、單位體積內總的注入量相同.如圖2所示，零維過程中電子數密度隨著等離子體注入量的增加迅速升高，在10 μs左右達到峰值，隨后開始迅速下降，到150 μs左右趨于穩定.而一維等離子體過程中并沒有出現零維過程中那樣明顯的電子數密度的峰值，而是在10 μs左右增長趨勢減緩，到50 μs左右開始出現衰減，并趨于穩定.零維等離子體過程隨著注入量的增加電子數密度所達到的峰值明顯比一維等離子體過程的最大值要高，且零維過程也比一維過程衰減得要快.這種現象說明，空間過程對低氣壓大氣等離子體的演化起著很重要的作用.

圖3為流量1018cm－2·s－1的情況下電子數密度在不同時刻的空間分布.如圖所示，電子在空間過程的作用下隨著時間的推移不斷由等離子體產生源處向外擴散.在60 km高空由于空氣含量稀薄和氣壓較小的影響，電子向外擴散得很快，200 μs時在模擬的邊界(5 cm)處電子數密度已經達到了1011cm－3.并且在等離子體源附近，等離子體分布隨著時間的推移逐漸趨于穩定，越靠近源點趨于穩定的時間越短.

圖 4 為流量 1018cm－2·s－1的情況下數密度隨時間的演化圖，圖中虛線為數密度，實線為數密度.如圖所示，這兩種正離子數密度的時間演化過程與電子數密度具有相似的趨勢，都在10 μs增長趨勢減緩，在50 μs左右達到極大值.和數密度最大值的比值約為28∶37，這與注入的比例21∶79產生了很大的差別.

圖3 流量為1018cm－2·s－1時電子數密度的空間分布

圖4 流量為1018cm－2·s－1時數密度隨時間的演化實線為數密度，虛線為數密度

圖5 流量為1018cm－2·s－1時 N+2數密度的空間分布

圖6 流量為1018cm－2·s－1時數密度的空間分布

鑒于邊界條件選取的任意性和無法判斷是否達到穩態，用以上數值結果判斷達到穩態所需等離子體流的強度是不充分的，需要做進一步討論.由于邊界處等離子體已很稀薄，可以采用近似的模型.

4.半無限空間的反應擴散過程解析分析

我們所做的主要近似就是將復雜的化學反應簡單化，而且只考慮最具實用價值的電子組分.我們先考慮純擴散，再加入由于化學反應(主要是復合)引起的密度減少.

4.1.半無限空間的純擴散過程

對于等離子體發生器，若只考慮所產生粒子的純擴散過程而忽略化學反應和漂移過程，當有持續的粒子流注入時，其密度分布的時空演化可由如下的模型來描述:在一維半無限空間x∈(0，∞)中，在左邊界面x=0上有持續的粒子流注入，單位時間內單位面積上注入的粒子數為DQ(D為擴散系數).若初始的粒子數密度為均勻分布n0，則粒子數密度n(x，t)的定解問題為

可以證明此定解問題的解為［12］

我們來分析空間中某一個定點的粒子數密度隨時間的變化情況.

說明隨著時間的推移，半無限空間上任意一點的密度都在單調無限增大，即當t→∞時，n(x，t)→∞.而對有限的時間 t，當 x→∞ 時，n(x，t)→n0，即粒子注入的影響需要一定的時間才能到達確定的地點.

4.2.半無限空間的反應擴散過程

當等離子體中接近穩態分布時，特別在計算邊界處，反應消減項將占主體地位，可以由如下的定性模型刻畫反應擴散過程:

和較完備方程組比較，在方程中用一簡單的衰減項代替反應項.k=1/τ，τ是一等效壽命.

此定解問題的解為［12］

同樣，我們來分析空間中某一個定點的粒子數密度隨時間的變化情況.當 t→∞時，n(x，t)→Q，即隨著時間 t的增大，空間某一點的粒子密度趨向于一個穩定的空間分布.其峰值為而產生的等離子體總量可表示為二者乘積 Γτ，與擴散系數無關.所以，如果只考慮產生的等離子體總量，零維計算就可滿足需要，而為研究空間分布和密度，則需一維計算.這在具體應用中還是很重要的.

從這一模型我們可以做一些簡化的估計.首先是半寬度，它不受人為控制，完全決定于壽命和等效擴散系數.從 Vidmar［5］的數據看，在很寬的空間范圍，如果達到1012cm－3的電子密度，電子壽命是10－5s.而等效的擴散系數我們卻不知道.如果按電子分量為 1.76×107cm2/s算，這一半寬度為 13 cm，如果按離子分量為300 cm2/s計算，為0.5 mm.實際情況處于二者之間.為達到一定的電子密度n，所需電子流量為

4.3.從數值模擬結果和近似解析表達式估計穩態解

在(12)式中，我們知道電子壽命 τ=1/k可以設為10－5s.所以可以從計算所得到的邊界處的電子密度值n(l，t)求有效擴散系數D.

在t=2 ×10－4s，x=5 cm 處，n(l，t)=1 ×1010cm－3，計算得到D=4.8 ×104cm2/s，計算得到穩態時的峰值為1.4×1013cm－3，分布的半寬度為0.7 cm.

這一結果和圖3顯示的數值模擬結果接近，說明我們的模擬已經接近穩態分布.峰值數值較模擬結果略高，可能因為近似解析式在此處不甚適用.

所以，在我們所設的高度，要得到密度為1012cm－3，厚度1 cm左右的等離子體層，所需等離子體流為1018cm－2·s－1.至于由這一結果估計所需源的功率，可以參考文獻［13］.所得數值模擬結果的合理性是因為我們的模擬區域寬度(5 cm)遠大于最后得到的等離子體分布的半寬度(0.7 cm).

5.結 論

利用一維模型，用數值模擬的方法研究了有源大氣等離子體在低氣壓開放環境下的反應擴散過程，得到了考慮化學反應、擴散以及漂移共同作用下的大氣等離子體的主要成分隨等離子體注入流量的變化規律.使用一個近似的解析解與數值模擬結果銜接，估計達到穩態時的電子分布.

本文提供了一種研究開放環境下有源等離子體的反應擴散問題的方法，可以用以估計產生一定密度厚度等離子體所需要的流量，并進一步估計所需要的功率.

作者衷心感謝中國科學院物理研究所王龍先生對本工作的指導.

［1］Liang Z W，Sun H L，Wang Z J，Xu J，Xu Y M 2008Acta Phys.Sin.57 4292(in Chinese)［梁志偉、孫海龍、王之江、徐 杰、徐躍民2008物理學報57 4292］

［2］Alexeff I，Anderson T 2006IEEE Trans.Plasma Sci.34 166

［3］Roth J R 2005Phys.Plasmas12 057103

［4］Jiang Z H，Hu X W，Liu M H，Gu C L，Pan Y 2003Chin.Phys.Lett.20 885

［5］Vidmar R J 1990IEEE Trans.Plasma Sci.18 733

［6］Ouyang J M，Guo W，Wang L，Shao F Q 2004Chin.Phys.13 2174

［7］Fang T Z，Ouyang J M，Wang L 2005Chin.Phys.Lett.22 2888

［8］Ouyang J M，Shao F Q，Wang L，Fang T Z，Liu J Q 2006Acta Phys.Sin.55 4974(in Chinese)［歐陽建明、邵福球、王龍、房同珍、劉建全2006物理學報55 4974］

［9］Fang T Z 2006Ph.D.Dissertation(Beijing:Chinese Academy of Sciences)(in Chinese)［房同珍 2006博士學位論文(北京:中國科學院)］

［10］Wu J H，Han Q S 1988The Theory，Method and Application of Computational Fluid Dynamics(Beijing:Science Press)p103(in Chinese)［吳江航、韓慶書1988計算流體力學的理論、方法及應用 (北京:科學出版社)第103頁］

［11］Heicklen J 1976Atmospheric Chemistry(New York:Academic Press)p3

［12］Gu C H，Li D Q，Chen S X，Zheng S M，Tan Y J 2007Method of Mathematical Physics(secondary edition)(Beijing:Higher Education Press)p49(in Chinese)［谷超豪、李大潛、陳恕行、鄭宋穆、譚永基2007數學物理方法(第二版)(北京:高等教育出版社)第49頁］

［13］Wang L，Zhang H F，Shao F Q 2004Wuli(Physics)33 18(in Chinese)［王 龍、張海峰、邵福球2004物理33 18］

PACS:52.50.Dg，92.60.H － ，52.65.－y，92.70.Cp

Numerical simulation of reaction-diffusion process of air plasma with a plasma source in open atmospheric environment*

Yin Zeng-Qian1)?Zhao Pan-Pan2)Dong Li-Fang2)Fang Tong-Zhen3)
1)(Department of Mathematics and Physics，North China University of Electric Power，Baoding 071003，China)
2)(College of Physics Science and Technology，Hebei University，Baoding 071002，China)
3)(Institute of Physics，Chinese Academy of Sciences，Beijing 100190，China)
(Received 29 January 2010;revised manuscript received 26 April 2010)

A one-dimensional model is used to investigate the reaction-diffusion process of air plasma with a plasma source by numerical simulation in low pressure open atmospheric environment.The spatial evolutions of main charged species’densities in the chemical reaction，drift and diffusion processes are obtained for a given plasma flow.The results of numerical simulation are combined with an approximate analytical formula to estimate the necessary plasma flow for electron density when a stable state is formed，with which the necessary power of the plasma source also can be estimated.

air plasma，plasma source，numerical simulation，reaction-diffusion process

*國家自然科學基金(批準號:10705049，10975043)資助的課題.

*Project supported by the National Natural Science Foundation of China(Grant Nos.10705049，10975043).