有源等離子體在開放大氣環(huán)境下的反應(yīng)擴散過程研究*

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

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

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

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

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

有源等離子體在開放大氣環(huán)境下的反應(yīng)擴散過程研究*

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

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

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

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

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

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

大氣等離子體，等離子體源，數(shù)值模擬，反應(yīng)擴散過程

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

1.引 言

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

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

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

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

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

2.反應(yīng)擴散模型和數(shù)值計算方法

2.1.反應(yīng)擴散模型

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

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

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

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

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

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

其中kB為波爾茲曼常數(shù)，得

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

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

表1 主要反應(yīng)物列表

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

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

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

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

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

2.4.QSSA方法的應(yīng)用

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

3.計算結(jié)果分析

3.1.模擬條件

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

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

3.2.模擬結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

4.半無限空間的反應(yīng)擴散過程解析分析

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

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

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

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

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

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

4.2.半無限空間的反應(yīng)擴散過程

當?shù)入x子體中接近穩(wěn)態(tài)分布時，特別在計算邊界處，反應(yīng)消減項將占主體地位，可以由如下的定性模型刻畫反應(yīng)擴散過程:

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

此定解問題的解為［12］

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

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

4.3.從數(shù)值模擬結(jié)果和近似解析表達式估計穩(wěn)態(tài)解

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

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

這一結(jié)果和圖3顯示的數(shù)值模擬結(jié)果接近，說明我們的模擬已經(jīng)接近穩(wěn)態(tài)分布.峰值數(shù)值較模擬結(jié)果略高，可能因為近似解析式在此處不甚適用.

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

5.結(jié) 論

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

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

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

［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，F(xiàn)ang 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計算流體力學的理論、方法及應(yīng)用 (北京:科學出版社)第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數(shù)學物理方法(第二版)(北京:高等教育出版社)第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).