耦合電極直流等離子體物理場數(shù)值模擬研究

武劭恂, 司徒達(dá)志, 張子煒, 李 晴, 黃衛(wèi)星

(1. 四川大學(xué) 化學(xué)工程學(xué)院, 四川 成都 610065; 2. 中國廣核集團(tuán)中廣核研究院有限公司, 廣東 深圳 518124)

1 引 言

國內(nèi)近幾年對固體垃圾問題，包括高危固體廢棄物的重視，使探索出一種高效、清潔的處理方法變得日益重要。而等離子體技術(shù)具有能源利用率高、處理效果突出的特點(diǎn)，能夠滿足這一需求。目前在國外該技術(shù)已應(yīng)用到包括城市生活、醫(yī)院衛(wèi)生、工業(yè)廢置等諸多領(lǐng)域[1-2]。而該技術(shù)在應(yīng)用中卻存在一個問題，即由于等離子體溫度高于電極金屬的熔點(diǎn)，導(dǎo)致設(shè)備壽命減小，工作效率降低[3]。所以發(fā)展一種簡便可靠的耦合電極-電弧的研究方法則成為關(guān)鍵。

由于等離子體發(fā)生裝置尺寸小、結(jié)構(gòu)復(fù)雜，在實(shí)驗(yàn)或生產(chǎn)操作中難以對其電極狀態(tài)進(jìn)行觀測。除此之外，其具有的高溫-電磁相互耦合的特點(diǎn)更使得數(shù)據(jù)測量工作難以進(jìn)行。因此，數(shù)值模擬計算成為解決該問題的主要研究方法[4-6]。但同時，由于人們更多關(guān)注等離子體主體內(nèi)的物理場分布，而忽略了電極本身熱電子發(fā)射機(jī)制帶來的熱效應(yīng)，因此在傳統(tǒng)計算結(jié)果中對于電極溫度的計算需要改進(jìn)。有鑒于此，且考慮轉(zhuǎn)移弧與非轉(zhuǎn)移弧發(fā)生原理相似，因此本文針對自由燃燒電弧等離子體，將電極與電弧區(qū)耦合建立二維軸對稱的磁流體動力學(xué) (magnetohydrodynamic, MHD) 計算模型，并通過合理的設(shè)計方法，通過求解電子連續(xù)方程與附加面熱源將電弧對電極的熱效應(yīng)、流體域的電導(dǎo)率修正包含在內(nèi)，得到了與實(shí)驗(yàn)測量相近的溫度場分布。該模型的建立可以為工業(yè)用直流電弧等離子體炬的設(shè)計改良提供方法基礎(chǔ)。

2 數(shù)學(xué)模型

2.1 基本假設(shè)

研究宏觀條件下發(fā)生裝置內(nèi)部的溫度場分布，因此通過下面的合理假設(shè)對復(fù)雜的電弧模型進(jìn)行簡化，具體包括：

1) 在電弧區(qū)域采用局域熱平衡 (local thermodynamic equilibrium, LTE) 假設(shè)；

2) 等離子體為光學(xué)薄；

3) 等離子體的密度、比熱、速度、黏度、熱導(dǎo)率、電導(dǎo)率、雙極擴(kuò)散系數(shù)、凈輻射系數(shù)、氣體粒子數(shù)密度分布等熱力學(xué)參數(shù)和輸運(yùn)系數(shù)都僅為溫度的函數(shù)；

4) 忽略重力和黏性耗散；

5) 由于等離子體溫度高，對應(yīng)的聲速也高，使得當(dāng)?shù)亓鲃玉R赫數(shù)并不高，非特別說明時等離子體流動是遠(yuǎn)低于聲速的，視為不可壓縮流動，可忽略壓力功；

以上假設(shè)可以將耦合模型中控制方程的諸多項(xiàng)進(jìn)行簡化，并能夠合理地選擇需要的變物性函數(shù)關(guān)系。

2.2 控制方程

本文計算模型以自由燃燒電弧為基礎(chǔ)，使用定常、二維軸對稱條件，通過對納維-斯托克斯方程組、Maxwell 方程組各項(xiàng)進(jìn)行增減，并補(bǔ)充近電極區(qū)的電導(dǎo)率修正以及電極表面熱源，得到以下耦合電極的磁流體動力學(xué)方程組(MHD)。

(1)納維-斯托克斯方程組

式中zv 、rv 分別為速度的軸向與徑向分量；Bθ為磁感應(yīng)強(qiáng)度的周向分量；zj 、rj 分別為電流的軸向與徑向分量；T 為溫度；ek 為工質(zhì)氣體的導(dǎo)熱系數(shù)；pc 為定壓比熱；σ 為電導(dǎo)率；Nε 為凈輻射系數(shù)；kB為玻爾茲曼常數(shù)；e 為電子電荷。物理量單位均使用國際單位制。

式(2)、(3)為動量方程，添加了洛侖茲力項(xiàng)，以實(shí)現(xiàn)磁場中帶電粒子受到洛侖茲力作用而產(chǎn)生的“箍縮”現(xiàn)象，這一現(xiàn)象發(fā)生在自由燃燒電弧中時，由于電弧受力向軸線方向壓縮，導(dǎo)致電流流通面積變小，電流密度增大，其產(chǎn)生的焦耳熱也會增大，因此不可忽略。而發(fā)生在非轉(zhuǎn)移弧時，由于存在固壁約束，電流的流通面積在氣流與洛侖茲力的共同作用下減少，此時洛侖茲力項(xiàng)可忽略。

(2) Maxwell 方程組

圖1 氬等離子子體凈輻射系數(shù) Fig.1 Profile of net emissionn coefficient of argon plasma

其中 σeff為有效電導(dǎo)率， n0為中性性粒子數(shù)密度， nT為總粒子數(shù)， μe為電子遷移率率，σ 為熱平衡狀態(tài)下的電導(dǎo)率。值得注意的是，除了電電子遷移率 μe以外，粒子數(shù)密度等都可以看作溫度度的函數(shù)，而電子遷移率，根據(jù)其定義

其中，em 為電子質(zhì)量，ev 為粒子碰撞頻率，thv 為電子熱運(yùn)動速度，eiQ 為電子-離子碰撞截面，eaQ 為電子-原子碰撞截面，這些參數(shù)并沒有溫度函數(shù)關(guān)系，意味電子遷移率eμ 需要用額外的方法去求解，本文使用了Bolsig 程序[12]求解玻爾茲曼方程得到了系列電子遷移率eμ ，用來修正電導(dǎo)率。

除了宏觀等離子體對電極的熱傳導(dǎo)、熱輻射外，電弧的行為也影響到電極的能量傳遞過程。陰極的過程具體表現(xiàn)為陰極金屬發(fā)射的電子具有一定的電子焓熱，使陰極表面溫度降低。同時陰極表面的電離層產(chǎn)生的離子，將會在較大的電場作用下向陰極表面運(yùn)動從而帶來一定的離子焓熱，使陰極表面溫度升高，除此之外還有本身的熱輻射，幾種因素共同作用，作為電極表面的額外能量輸入。陽極的過程表現(xiàn)為電子運(yùn)動到陽極表面釋放焓熱，存在熱輻射，而沒有離子焓。電極的面熱源形式如下[13]：

其中，ej 為近陰極區(qū)的電子電流，cφ 為陰極的逸出功，ij 為近陰極區(qū)的離子電流，iV 為工質(zhì)氣體一次電離電勢，ε 為表面發(fā)射系數(shù)，α 為斯特藩-玻爾茲曼常數(shù)。

在具體計算中，由于離子電流無法提前得知，需要根據(jù)電流連續(xù)性條件，將陰極金屬發(fā)射產(chǎn)生的熱電子電流與弧柱區(qū)求解Maxwell 方程組得到的電弧電流做比較，繼而通過離子電流補(bǔ)全守恒關(guān)系。熱電子發(fā)射方程(Richardson 方程)為

其中， A 為熱發(fā)射常數(shù)， kB為玻爾茲曼常數(shù)。離子電流 ji計算方法為，當(dāng)弧區(qū)電流大于發(fā)射電流時( jn＞ jR)， je= jR， ji= jn- jR；當(dāng)弧區(qū)電流小于發(fā)射電流時( jn＜ jR)， je= jn， ji= 0。

但由于在Fluent 軟件中，無法進(jìn)行面熱源的設(shè)置，式(9)、(10)無法直接使用，因此采用另外一種，類似激光加熱金屬表面的處理方法[14]對其進(jìn)行等效轉(zhuǎn)化，從而能夠正常加載。

2.3 模型歸納

2.2 節(jié)對該耦合模型的內(nèi)容以及設(shè)置方法進(jìn)行了描述，在此，為了方便與Fluent 軟件程序相接，將上述方程組歸納為統(tǒng)一的控制方程形式，即基本的NS 方程形式。

表1 項(xiàng)類歸納 Table 1 Components of equations

根據(jù)以上歸納，即可使用戶自定主標(biāo)量(user-defined scalar，UDS)定義新的控制方程，使用用戶自定義函數(shù)(user-defined function，UDF)寫入物性函數(shù)、邊界條件等。

3 計算模型

3.1 模型劃分及邊界條件

本文采用二維軸對稱模型，圖2 為模型結(jié)構(gòu)，為了顯示陰極平臺，圖示已作適當(dāng)放縮，取半徑35 mm，長40 mm 的計算域，其中包含陰極固體域、流體域、陽極固體域3 個部分。陰極固體域主體半徑1.6 mm，尖端平臺為30 度半錐角，平臺半徑為0.1 mm，材料為釷鎢合金；流體域使用氬氣工質(zhì)氣體；陽極固體域?yàn)榇笃桨?，材料為銅。由于本文只探究電極溫度分布，對于溫度超過金屬熔點(diǎn)的情況不添加熔化-凝固模型，一律默認(rèn)水冷材料，同時水冷材料不會影響電弧特性[9]。

網(wǎng)格建立過程中，采取先疏后密，并不斷調(diào)節(jié)亞松弛因子，每階段計算收斂后，使用自適應(yīng)網(wǎng)格，將溫度梯度較大區(qū)域以及電弧主要分布區(qū)域網(wǎng)格加密，不斷調(diào)整計算最終達(dá)到要求。

圖2 耦合電極自由燃燒電弧計算模型 Fig.2 Calculation model for coupling electrode free combustion arc

本文采用200 A 弧電流，電流連續(xù)方程與磁矢量勢方程在全計算域計算，電子連續(xù)方程只在流體域計算，使用的工質(zhì)氣熱力學(xué)參數(shù)與輸運(yùn)參數(shù)分別來自MURPHY[15]與張曉寧[16]。

所用程序使用Visual Studio2017 編寫，F(xiàn)luent18.1 求解，除了直接得到上述標(biāo)量值，對計算域中總(分)電流值、焦耳熱值、熱輻射值、電場強(qiáng)度等黑箱值均單獨(dú)分配存儲空間UDM 以便分析。

表2 邊界條件 Table 2 Boundary conditions

3.2 計算過程

該模型需要將電磁場與一般流場涉及到的物理模型耦合計算，需要使用Fluent 軟件的UDF 與UDS功能，用以引入某些物性參數(shù)隨溫度變化的函數(shù)，并添加自定義標(biāo)量電勢φ 、磁矢量勢的軸向與徑向分量zA 、rA，以及電子數(shù)密度en 對應(yīng)的標(biāo)量方程。

通過上文對控制方程組的描述可以得知，除了電子連續(xù)方程較為特殊之外，其余方程求解中需要使用到的物性參數(shù)均與溫度分布密切相關(guān)，因此從溫度來設(shè)想一次迭代過程如下：由初始流場溫度計算出局域熱平衡電導(dǎo)率σ ，同時得到當(dāng)前溫度下的雙極擴(kuò)散系數(shù) AD 、局域熱平衡電子數(shù)密度eqn 用以求解式(8)電子連續(xù)方程，得到當(dāng)前溫度下的電子數(shù)密度en ，該溫度下還可求解出1σ 參與的電流連續(xù)方程得到電勢，即得到電場強(qiáng)度分布，以上計算結(jié)果可以得到有效電導(dǎo)率effσ ，于是其他過程均以有效電導(dǎo)率effσ 為基礎(chǔ)繼續(xù)計算。可以發(fā)現(xiàn)，雖然電磁場與流場耦合計算，但是電子方程和以電勢為基礎(chǔ)的方程組的求解相對獨(dú)立，溫度作為影響因素僅對電導(dǎo)率的計算有影響，而電子方程的求解反過來又會影響電勢方程的求解，所以在大的耦合模型下，這兩個方程又作為另一對耦合形式出現(xiàn)。此外，各方程添加的源項(xiàng)，極易引起迭代計算的震蕩或者發(fā)散，因此為了得到穩(wěn)定且較為準(zhǔn)確的物理場分布，一方面需要保證高質(zhì)量的網(wǎng)格以及合理的亞松弛因子，另一方面各項(xiàng)參數(shù)采集以及形式變化需要經(jīng)過仔細(xì)斟酌。

4 計算結(jié)果分析

4.1 速度-壓力分布

氣體進(jìn)口處速度較低，但由于陰極端部電流密度很大，工質(zhì)氣被迅速加熱電離后，等離子體在洛侖茲力的作用下壓縮，熱量集中，從而使等離子體加速形成射流。如圖3、4，在陰極附近，即溫度最高的區(qū)域形成了一個由于熱膨脹造成的高壓區(qū)與陽極軸線附近形成了動能-勢能轉(zhuǎn)化高壓區(qū)。

圖3 速度分布 Fig.3 Map of speed distribution

圖4 壓力分布 Fig.4 Map of pressure distribution

4.2 溫度分布

電極耦合模型的固體域溫度分布是預(yù)估陰極壽命的重要依據(jù)，因此在自由燃燒電弧模型的基礎(chǔ)上需要補(bǔ)充陰極金屬發(fā)射電子對流體域電弧計算的影響，以及陰極弧根的形位對陰極表面的附近熱源的影響。圖5 中該模型計算得到流體域最高溫度位于陰極尖端附近，約為23 150 K，陰極最高溫度位于陰極尖端，約為4 000 K，陽極最高溫度位于軸線附近，且呈近二次項(xiàng)分布，約為850 K。溫度分布范圍與圖6 中前人研究結(jié)果基本一致?？梢园l(fā)現(xiàn)，在不采用電極水冷或添加劑熔融擴(kuò)散的前提下[13]，采用200 A 電流強(qiáng)度勢必會造成陰極的燒蝕。由于在計算過程中，氣體逐漸電離，電極溫度逐漸升高，可以設(shè)想，在此計算模型的基礎(chǔ)上，如果采用非穩(wěn)態(tài)模型計算，當(dāng)陰極溫度升高至熔點(diǎn)時所用計算時間即為電極的工作壽命，繼而可以加入水冷與熔化-凝固模型對實(shí)際設(shè)備進(jìn)行較為準(zhǔn)確的設(shè)計優(yōu)化，這也是本文的出發(fā)點(diǎn)。

圖5 計算溫度分布 Fig.5 Calculated temperature distribution

圖6 文獻(xiàn)中溫度分布[11] Fig.6 Temperature distribution results from literature[11]

4.3 電子數(shù)密度分布

此處的電子數(shù)密度為求解電子連續(xù)性方程得到，由于該標(biāo)量方程涉及到的物性參數(shù)較多，數(shù)量級普遍較大，采集時存在一定誤差，此處不再展示其余粒子數(shù)密度等分布圖。

圖7 電子數(shù)密度分布Fig.7 Map of electronic number density distribution

由圖7 中計算結(jié)果可知，電子數(shù)密度分布規(guī)律與溫度分布規(guī)律近似，并且16 000～23 000 K 區(qū)域?qū)?yīng)電子數(shù)密度最大的區(qū)域，這是因?yàn)橐淮坞婋x的氬離子在14 000 K 左右會發(fā)生二次電離，但在25 000 K 以下由二次電離產(chǎn)生的電子數(shù)與一次電離提供的電子存在數(shù)量級的差別，直到高于25 000 K 時，一次電離的氬原子數(shù)才會呈現(xiàn)明顯下降。所以，雖然本文中計算溫度最高達(dá)到了23 000 K，但是使用的粒子數(shù)相關(guān)物性仍采用一次電離數(shù)據(jù)，使用的氬氣電離能也采用了一次電離能，計算電子數(shù)密度分布規(guī)律基本合理。

4.4 電勢-電流密度分布

與電勢、電流密切相關(guān)的參數(shù)是電導(dǎo)率，由以上計算得到的電子數(shù)密度對局域熱平衡電導(dǎo)率進(jìn)行修正，最終得到電勢為13.6 V(圖8)，電流密度最大為7×108A?m-2位于陰極尖端附近(圖9)。由于本文關(guān)于近陰極區(qū)的計算未包含鞘層的電位降計算，所以采用電勢補(bǔ)償?shù)姆椒╗17]，由ZHOU 等[18]的研究，當(dāng)電流為200 A 大小時，鞘層內(nèi)電勢差約為11.7 V，則補(bǔ)償后電勢差為25.3 V，與實(shí)驗(yàn)測量值26 V 較為接近。相較于一般電導(dǎo)率σ 為溫度的單調(diào)函數(shù)，有效電導(dǎo)率effσ 在修正過程中需要使用電子遷移率eμ ，這一參數(shù)在很大程度上與電場強(qiáng)度E 大小相關(guān)。由電弧的結(jié)構(gòu)可知，在近電極區(qū)的較小軸向距離內(nèi)存在很高的電勢差，這意味著該區(qū)域的場強(qiáng)較大，因此，場強(qiáng)分布非線性，導(dǎo)致除了氣體完全未電離和完全電離兩種情況外，修正后的電導(dǎo)率分布關(guān)于溫度并不是單調(diào)關(guān)系，由此導(dǎo)致計算結(jié)果里電勢的數(shù)值會低于一般自由燃燒電弧，而電弧主體溫度值接近，其中差別最明顯的就是在近陰極區(qū)的溫度分布。

圖9 電流密度密度分布Fig.9 Map of density distribution of current density

5 結(jié) 論

針對耦合電極的直流等離子體進(jìn)行了二維軸對稱數(shù)值模擬。主要結(jié)論如下：

1) 速度-壓力變化關(guān)系與一般自由燃燒電弧相似，速度高值位于電弧中心，壓力高值位于電弧兩端；

2) 流體域的溫度-電子數(shù)密度分布規(guī)律基本一致，符合氬氣正常電離過程的數(shù)量變化規(guī)律；固體域溫度結(jié)果合理，可以近似預(yù)測電極壽命；

3) 提出了針對Fluent 軟件功能的計算模型，以及在實(shí)驗(yàn)數(shù)據(jù)缺乏下的近似處理方法。