鐘武燁，呂 征，鄭劍平，楊啟法

(中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

熱離子反應(yīng)堆電源是采用熱離子能量轉(zhuǎn)換方式將反應(yīng)堆裂變釋放的熱能直接轉(zhuǎn)換成電能的能量轉(zhuǎn)換裝置，其具有比功率高、無機械轉(zhuǎn)動、所需比散熱面積小等優(yōu)點，是國際上最具優(yōu)勢的高效空間核電源之一[1]。目前，具有實際應(yīng)用價值的熱離子能量轉(zhuǎn)換器的結(jié)構(gòu)形式是充銫的工作在電弧工況的熱離子二極管[2]，其通常采用難熔金屬作電極對，發(fā)射極被加熱到高溫(～1 500 ℃)產(chǎn)生熱電子發(fā)射，充入電極間隙的銫蒸氣在其中電離形成等離子體，電子輸運通過等離子體區(qū)域后達(dá)到接收極，經(jīng)外電路的負(fù)載做功后重新回到發(fā)射極，形成回路。熱離子轉(zhuǎn)換器可看作以電子作為介質(zhì)的熱機[3]，電極間隙的勢能分布則決定該熱機的能量來源。因此，電子在電極間隙的勢能分布是研究熱電轉(zhuǎn)換特性的關(guān)鍵。

目前，對關(guān)于電極間隙的電子勢能分布已有較深入的了解，但一些基礎(chǔ)問題仍有待進(jìn)一步深入研究，如熱離子轉(zhuǎn)換器運行在電弧工況不同區(qū)域(包括飽和區(qū)、阻塞區(qū))時，對于電極鞘層的電子勢能分布的認(rèn)識仍不夠深入，忽略對阻塞區(qū)虛擬陰極的討論[4]；在討論阻塞區(qū)虛擬陰極時又無法給出數(shù)值上可解的方法[5]；缺乏對接收極電勢分布特征的研究等。導(dǎo)致該問題的原因主要是現(xiàn)有的模型并沒有很好地處理等離子體輸運及其邊界條件的匹配問題。本文基于動理學(xué)方法建立電弧工況下的等離子體輸運模型，對等離子體鞘層的電子勢能分布進(jìn)行分析，并將等離子體鞘層處理為輸運的邊界條件，自編程實現(xiàn)等離子體輸運方程與邊界條件的自適應(yīng)解耦求解。

1 理論建模

電子勢能分布的研究需建立等離子體的輸運模型。當(dāng)熱離子轉(zhuǎn)換器工作在電弧工況(圖1)時，在電極間隙中，等離子體與電極表面相接觸的區(qū)域形成電勢突然躍變的、非準(zhǔn)中性的等離子體鞘層，其尺度與德拜長度相當(dāng)，遠(yuǎn)小于電極間隙。因此，電極間隙的等離子體可分為電極附近的鞘層和中間的準(zhǔn)中性等離子體區(qū)[6]。在具體的建模過程中，在準(zhǔn)中性等離子體區(qū)建立帶電粒子的輸運模型，而等離子體鞘層則處理為輸運的邊界條件。

1.1 電弧工況等離子體輸運模型

熱離子能量轉(zhuǎn)換器電極間隙的等離子體屬于低溫弱電離的等離子體，此時必須考慮電子、離子以及中性原子3組分同時存在的多粒子碰撞輸運問題。等離子體中帶電粒子體系間存在復(fù)雜的相互作用，要精確描述等離子體的行為是極其困難的。目前，只能根據(jù)不同的條件和研究的問題，采用不同的近似方法[7]。文獻(xiàn)給出的建模方法可歸為唯象模型和基于介觀尺度的動理學(xué)模型兩大類[8]。唯象模型作為定性分析是有效的，作為定量分析計算仍有差距。本文的輸運方程采用動理學(xué)的模型。動理學(xué)模型是輸運過程嚴(yán)格的定量理論。動理學(xué)方法研究其輸運過程的基本思路是從粒子碰撞遷移所滿足的玻爾茲曼方程出發(fā)，然后根據(jù)轉(zhuǎn)換器等離子體輸運的具體特點對分布函數(shù)和碰撞輸運方程作合理簡化，建立輸運控制方程組，并給出輸運系數(shù)。

由于電極間隙遠(yuǎn)小于電極的軸向尺寸，故只需建立電極間隙準(zhǔn)中性等離子體輸運的一維模型，在建模過程中以發(fā)射極作為原點建立坐標(biāo)系。中性銫原子按理想氣體處理：靠近電極的中性原子溫度等于電極溫度，電極間隙各點的溫度隨電極溫度呈線性變化，不必考慮中性原子的輸運問題，但必須考慮由于中性原子的存在導(dǎo)致的散射效應(yīng)。等離子體輸運的描述的目的是獲得電極間隙等離子體濃度n、電子溫度Te、負(fù)電勢V、電子電流密度je及其能量流密度Se、離子電流密度ji及其能量流密度Si等分布狀態(tài)。

圖1 熱離子能量轉(zhuǎn)換器的工況分類及工作原理Fig.1 Operation mode and working principle of thermionic energy converter

首先，建立電子、離子、中性原子3組分共存條件下的玻爾茲曼方程：

式中：β=e,i,a(e、i、a分別表示電子、離子、中性原子)；α=e,i；v為速度；r為位置；e為電量；下標(biāo)in表示非彈性碰撞；Jαβ為碰撞積分項；fα為分布函數(shù)并假定其偏離平衡態(tài)不遠(yuǎn)。通過對玻爾茲曼方程的簡化，可建立等離子體參量的通量及連續(xù)性方程[9-10]。

1) 通量方程。帶電粒子的通量包括引起帶電粒子輸運的濃度梯度、溫度梯度以及弱電場力的貢獻(xiàn)，電子、離子的電流密度分別為：

為考慮離子電流密度的貢獻(xiàn)，總電流密度jtot包括電子電流密度和離子電流密度的部分。

jtot=je-ji

對于帶電粒子的能量通量，考慮了宏觀運動的動能、勢能以及熱擴(kuò)散的貢獻(xiàn)。

式中：k為玻爾茲曼常數(shù)；De,i、μe,i、ke,i、βe,i分別為濃度擴(kuò)散系數(shù)、遷移率、熱導(dǎo)率、熱擴(kuò)散比；Rie為離子受電子的作用力。

2) 連續(xù)性方程。帶電粒子通量的連續(xù)性是由體電離效應(yīng)決定的。

式中：Γ為電離函數(shù)；σ0(Te)為與電子溫度相關(guān)的等效電離截面；Na為原子濃度；n(Te)為熱平衡態(tài)下的等離子體濃度。

對于能量通量的連續(xù)性，重粒子在電離/復(fù)合過程中吸收/放出的能量由電子提供/帶走，忽略電子的輻射ΔSrad以及e-i、e-a的碰撞損失ΔSei、ΔSea，則有：

ΔSei-ΔSea=-Γ(Eion-eV)

式中，Eion為銫離子的第一電離能。

上述的通量及連續(xù)性方程組成了等離子體輸運的控制方程組，其中包含了je、ji、Se、Si、n、V、Te7個變量。對于其中的輸運系數(shù)，當(dāng)分布函數(shù)偏離平衡態(tài)不遠(yuǎn)時，基于動理學(xué)理論可對分布函數(shù)作Sonine多項式展開，從而將輸運系數(shù)分級展開，取其中的前3項即可使截斷誤差小于10%[11]。

1.2 電弧工況等離子體輸運的邊界條件

圖2 電弧工況飽和區(qū)可能的電子勢能分布Fig.2 Possible distribution of electron potential energy at saturation region of ignited mode

為獲得等離子體輸運微分方程的邊界條件，需建立帶電粒子通量、能量通量的邊界值與等離子體邊界參量的關(guān)系，其中的關(guān)鍵點是等離子體鞘層的電子勢能分布。在電弧工況的不同區(qū)域，等離子體鞘層電子勢能的分布狀態(tài)是有差異的，所有可能的勢能分布狀態(tài)如圖2、3所示，其中VE,C分別為發(fā)射極、接收極鞘層電勢的躍變，φE,C分別為發(fā)射極、接收極功函數(shù)，Ud、U分別為電弧壓降及輸出電壓。對于發(fā)射極鞘層，其分布在飽和區(qū)是單調(diào)下降的(圖2)；而在阻塞區(qū)則是先增后降(圖3)，形成了電子發(fā)射的勢壘ΔVE。另一方面，對于接收極鞘層，其電勢分布并不依賴于具體的工作區(qū)域，其正負(fù)性取決于邊界條件本身的匹配。以下依據(jù)電極鞘層內(nèi)電子勢能的分類對邊界條件進(jìn)行討論，并將接收極表面的電勢定義為零電勢參考點。

圖3 電弧工況阻塞區(qū)可能的電子勢能分布Fig.3 Possible distribution of electron potential energy at obstructed region of ignited mode

1) 發(fā)射極鞘層的邊界條件

電弧工況發(fā)射極鞘層的邊界條件需區(qū)分為飽和區(qū)與阻塞區(qū)兩種情況進(jìn)行討論。

在飽和區(qū)，發(fā)射極附近的電子勢能是單調(diào)下降且躍變較大，由此產(chǎn)生顯著的附加電場促進(jìn)發(fā)射，即產(chǎn)生肖脫基效應(yīng)，從而獲得較大的輸出電流密度[12]。邊界通量需考慮帶電粒子的熱發(fā)射以及來自于等離子體的反射。

jeE=(1-r1)jSE(EE)-

式中：jeE、jiE、SeE為發(fā)射極鞘層與等離子體交界面上的通量；TeE為交界面上的電子溫度；nE為交界面上的等離子體濃度；jSE(EE)為熱電子的肖脫基發(fā)射項；r1、r2為考慮了發(fā)射的動理學(xué)反射特性而引入的參數(shù)[13]。

在阻塞區(qū)，發(fā)射極鞘層的電子勢能是先增后降的，即熱電子發(fā)射的勢壘增大了ΔVE，形成虛擬陰極，由于未促進(jìn)發(fā)射附加電場，肖脫基效應(yīng)也就不存在[14]，此時的發(fā)射電流密度是飽和發(fā)射電流密度jS0。

2kTeE+je(Ud+ΔVE)

2) 接收極鞘層的邊界條件

當(dāng)接收極的電勢躍變VC>0時(對應(yīng)于圖2a與圖3a)，有：

式中：jeC、jiC、SeC為接收極鞘層與等離子體交界面上的通量；f0與f1為描述接收極鞘層內(nèi)電子運動的各向異性而引入的函數(shù)[15]。

當(dāng)接收極的電勢躍變VC<0時(對應(yīng)于圖2b與圖3b)，有：

2 模型求解與算法流程

根據(jù)上述建立的輸運方程和采用的邊界條件，輸運系數(shù)依賴于等離子體參量本身，因此輸運方程本身是非線性的，故本文建立的模型只能數(shù)值求解。另一方面，輸運的邊界條件依賴于工作參數(shù)所處的工況區(qū)域(飽和區(qū)或阻塞區(qū))，但實際計算中并不能根據(jù)輸入條件(工作參數(shù))預(yù)知其工況區(qū)域。因此，在數(shù)值求解的過程中，應(yīng)使得邊界條件與輸運方程自適應(yīng)匹配求解。本文采用的算法流程如圖4所示。

1) 輸入計算所必要的參數(shù)：基本物理參量、熱離子轉(zhuǎn)換器的運行參數(shù)，如電極溫度(TE、TC)、銫溫TCs、電極的吸銫功函數(shù)(φE、φC)、電極間隙d、總電流密度(輸出電流密度)jtot等。

圖4 熱離子轉(zhuǎn)換器電弧工況等離子體輸運算法流程Fig.4 Algorithm flow of plasma transportation of TEC at ignited mode

2) 給定接收極鞘層與準(zhǔn)中性等離子體區(qū)交界處的電子溫度TeC和接收極鞘層的電勢躍變VC作為迭代初始值，根據(jù)VC的正負(fù)性選擇相應(yīng)的接收極邊界條件模型，計算此邊界的jeC、jiC、nC及SeC、SiC。

3) 以流程2計算得到的參量為邊界值，采用龍格-庫塔方法對輸運微分方程進(jìn)行積分，得到準(zhǔn)中性等離子體區(qū)Te、n、je、ji、V、Se、Si的分布，特別地，記發(fā)射極邊界上的電子電流密度積分為je_inte(TeC,VC)、離子電流密度積分為ji_inte(TeC,VC)。

4) 先假定工作點處于阻塞區(qū)，根據(jù)阻塞區(qū)發(fā)射極邊界條件計算非單調(diào)的電勢躍變情況VE及ΔVE。若解存在，則認(rèn)定虛擬陰極的存在(即處于阻塞區(qū))；若解不存在，則認(rèn)定處于飽和區(qū)(即虛擬陰極不存在，發(fā)射極電勢躍變是單調(diào)的)，并據(jù)此計算相應(yīng)的電勢躍變VE。

5) 根據(jù)邊界條件的自適應(yīng)判斷，計算相應(yīng)工況下發(fā)射極邊界上的電子電流密度je_cal(TeC,VC)、離子電流密度ji_cal(TeC,VC)，并與積分所得邊界值聯(lián)立成方程組je_cal(TeC,VC)-je_inte(TeC,VC)=0與ji_cal(TeC,VC)-ji_inte(TeC,VC)=0。采用牛頓迭代法求解方程組，直至收斂到設(shè)定的精度。若發(fā)散，則重新選擇初始條件進(jìn)行計算。

6) 輸出計算結(jié)果。

3 結(jié)果討論

由于程序的輸入包含電極溫度、電極功函數(shù)、銫溫、電極間隙、總電流密度(輸出電流密度)等7個參數(shù)，在結(jié)果討論中不適宜采用控制變量的方式。Bullis和Rasor[16]在其綜述中指出，盡管熱離子轉(zhuǎn)換器的工作參數(shù)有差異，其輸出特征可采用離子裕度β及參量pCsd來區(qū)分。其中，離子裕度β由發(fā)射極的發(fā)射特性決定，用來表征熱發(fā)射離子中和熱發(fā)射電子的程度，其臨界值為1；參量pCsd由輸運過程決定，用來表征電子輸運過程的散射強度，其臨界值為20 mile·Torr(1 Torr≈133.322 Pa，1 mil=25.4 μm。

為在較寬的參數(shù)范圍內(nèi)驗證計算模型的準(zhǔn)確性，受限于電極材料，本文研究的熱離子轉(zhuǎn)換器是電極對在給定的電極溫度下運行，僅在可變的銫溫下進(jìn)行輸出，此時離子裕度β被限定而參量pCsd可在較寬的范圍取值。為便于與文獻(xiàn)的數(shù)據(jù)進(jìn)行比較[17]，程序計算中輸入的工作參數(shù)列于表1，選定的電極對材料為Re-Mo，限定電極的工作溫度TE=1 700 K、TC=925 K，電極間隙d=0.25 mm。根據(jù)表1，由于電極特性的制約，離子裕度β限制在小于1的區(qū)域，而參量pCsd則覆蓋了較寬的范圍。將上述工作參數(shù)作為程序的輸入，獲得本文所有的計算結(jié)果。

表1 程序計算中輸入的工作參數(shù)Table 1 Input working parameter for procedure

1) 計算結(jié)果的有效性判定

首先，為判定計算結(jié)果的合理性，需將等離子體的德拜長度rD(表2)與特征尺寸(電極間隙)進(jìn)行比較，德拜長度均遠(yuǎn)小于電極間隙，說明電極間隙形成了等離子體，所計算的工況都處于電弧工況；并且隨著總電流密度的增大，工況狀態(tài)逐漸從電弧工況的阻塞區(qū)進(jìn)入飽和區(qū)，在飽和區(qū)發(fā)射極鞘層電場值滿足EE遠(yuǎn)大于0。因此，計算結(jié)果驗證了邊界條件分區(qū)設(shè)置的合理性。

表2 計算結(jié)果有效性判定所需的數(shù)據(jù)Table 2 Datum needed to judge validity of computing result

2) 伏安特性

圖5示出本文計算得到的伏安特性及其與文獻(xiàn)值的對比，計算值與實驗值符合較好。圖6比較了TCs=551 K與591 K時的伏安特性，計算表明，當(dāng)參量pCsd<20 mile·Torr時，飽和區(qū)的輸出電流密度表現(xiàn)出較明顯的飽和特征；而當(dāng)參量pCsd>20 mile·Torr時，飽和電流密度的特征不明顯。

圖5 本文計算的伏安特性曲線與文獻(xiàn)值的對比Fig.5 Comparison of U-I characteristic between this paper and lecture results

圖6 本文計算的不同銫溫時伏安特性的比較Fig.6 Computed U-I characteristic in this paper at different temperatures of cesium

3) 電子勢能分布

當(dāng)熱離子轉(zhuǎn)換器的銫溫為551 K時，電極間隙電子勢能的變化如圖7所示，涵蓋了電子從發(fā)射極費米面熱發(fā)射輸運到接收極費米面的全過程。整體而言，電子勢能在發(fā)射極鞘層的躍變大于接收極鞘層的躍變，而在準(zhǔn)中性等離子體區(qū)則相對平緩(其中躍變發(fā)生的尺度僅有德拜長度的量級，為便于展示，圖片對躍變的尺度進(jìn)行放大)。

對于特定的工作點(例如jtot=6 A/cm2)，電子始末位置的能級差，即兩電極費米面的能級差，構(gòu)成了對外輸出電壓。兩電極表面勢能的差用于維持電弧狀態(tài)，屬于固有的內(nèi)部壓降。

圖7 TCs=551 K時電子勢能隨電極間位置的變化Fig.7 Dependence of electron potential energy on location of interelectrode at TCs=551 K

對于發(fā)射極鞘層的電勢躍變，當(dāng)電流密度較小時(jtot=1 A/cm2)，發(fā)射極鞘層的電子勢能是先增后減的非單調(diào)躍變，在發(fā)射極的鞘層內(nèi)形成了附加的發(fā)射勢壘(即所謂的虛擬陰極)，相當(dāng)于增加了發(fā)射極功函數(shù)，從而抑制了電子發(fā)射，但發(fā)射極功函數(shù)的提升有利于獲得更大的電壓輸出，此時轉(zhuǎn)換器工作在阻塞區(qū)。隨著電流密度的增大，發(fā)射極鞘層的電子勢能變?yōu)閱握{(diào)向下躍變，其量值隨電流密度而增大，此勢能躍變使電子熱發(fā)射的過程中得到電場的激勵(即肖脫基效應(yīng))，同時該電場較好地阻礙等離子體中的電子返回發(fā)射極；為維持大電流密度，就必然需要更高的電弧壓降(內(nèi)壓增大)，從而輸出電壓也就降低了，此時轉(zhuǎn)換器工作在飽和區(qū)。

圖8 TCs=551 K時發(fā)射極鞘層電子電流密度的構(gòu)成Fig.8 Consistency of electron current density in sheath of emitter at TCs=551 K

整體而言，在輸出電流密度遞增的過程中，發(fā)射極鞘層的電子勢能分布經(jīng)歷了由非單調(diào)變成單調(diào)的過程。對發(fā)射極鞘層電子電流密度的構(gòu)成分析可進(jìn)一步說明產(chǎn)生此分布特征的機理，如圖8所示，在小電流密度的阻塞區(qū)，從等離子體返回的電子電流密度超過凈電子電流密度，相當(dāng)于發(fā)射極鞘層“累積”了電子，這將增加電子發(fā)射的勢壘，形成了虛擬陰極；在大電流密度的飽和區(qū)，熱發(fā)射的電子電流密度遠(yuǎn)大于從等離子體返回的電子電流密度，此時在發(fā)射極附近產(chǎn)生強烈的電離(圖9a)，發(fā)射極附近產(chǎn)生大量離子，形成了促進(jìn)電子發(fā)射的單調(diào)向下的電子勢能分布。

對于接收極鞘層的電勢躍變，在輸出電流密度遞增的過程中，其電子勢能分布經(jīng)歷了從向下躍變(VC<0)到向上躍變(VC>0)的過程。產(chǎn)生這種分布的原因是在輸出電流密度從0開始增大的過程中，實驗觀察到最先發(fā)生電離的區(qū)域是接收極的鞘層[18]，由此從擴(kuò)散工況進(jìn)入飽和工況的阻塞區(qū)，因此需要向下躍變的電子勢能分布促使接收極附近發(fā)生電離(圖9b)。而隨著輸出電流密度的增大，電離的峰值逐漸向發(fā)射極移動，接收極一側(cè)趨向于電離平衡(圖9a)，此時接收極鞘層表現(xiàn)為準(zhǔn)中性等離子體鞘層的性質(zhì)，即熱速率較大的電子運動到接收極表面逸走而留下富余的正離子，從而使電子勢能向上躍變。

以上關(guān)于電極間隙電子勢能的分布特性，在轉(zhuǎn)換器的銫溫為591 K時也同樣存在(圖10)。

圖9 TCs=551 K時電離速率隨電極間位置的變化Fig.9 Dependence of rate of ionization on location of interelectrode at TCs=551 K

圖10 TCs=591 K時電弧工況不同工作點下對應(yīng)的電子勢能分布Fig.10 Dependence of electron potential energy at different working points of ignited mode at TCs=591 K

4) 根據(jù)輸出功率特性對電子勢能分布特性進(jìn)行確認(rèn)

為確認(rèn)上述的熱離子轉(zhuǎn)換器電極間隙電子勢能分布特性，進(jìn)行最大輸出功率的討論。一般而言，當(dāng)電池的內(nèi)外壓降相等時有最大的電功率輸出，這是電池的一基本屬性。根據(jù)電子勢能分布，可獲得輸出電壓、電弧壓降與電流密度的關(guān)系，如圖11所示，電流密度越大，輸出電壓越小，而電弧壓降則越大，兩者的交點為內(nèi)外壓相等的工作點。同時，根據(jù)計算的伏安特性，可獲得輸出電功率密度與電流密度的關(guān)系(圖12)，當(dāng)轉(zhuǎn)換器銫溫為551、591 K時，在對應(yīng)內(nèi)外壓相等的工作點jtot=5、14 A/cm2下獲得最大的輸出電功率。這說明了本文關(guān)于熱離子轉(zhuǎn)換器的計算結(jié)果是自洽的，并與一般電池的功率特性相一致。

圖11 輸出電壓、電弧壓降與輸出電流密度的關(guān)系Fig.11 Dependence of output voltage and plasma voltage drop on output current density

圖12 輸出電功率密度與電流密度的關(guān)系Fig.12 Dependence of output power density on output current density

4 結(jié)論

本文基于動理學(xué)方法建立了描述熱離子轉(zhuǎn)換器電極間隙電弧工況等離子體輸運的方程及相應(yīng)的邊界條件，并設(shè)計了一套基于牛頓迭代法的計算程序?qū)崿F(xiàn)其自適應(yīng)解耦求解，詳細(xì)討論了等離子體的電子勢能分布特性。

1) 電子勢能在發(fā)射極鞘層的躍變大于接收極鞘層的躍變，而在準(zhǔn)中性等離子體區(qū)則相對平緩。隨著輸出電流密度的增大，發(fā)射極鞘層的電子勢能躍變特征將從先增后減的非單調(diào)變?yōu)閱握{(diào)下降，這是由發(fā)射極鞘層內(nèi)電子電流密度構(gòu)成的差異導(dǎo)致的。

2) 同時接收極鞘層的電子勢能則從向下躍變變?yōu)橄蛏宪S變，這是由接收極鞘層附近的電離狀態(tài)決定的。

3) 基于等離子體鞘層電子勢能分布的分析，改進(jìn)了等離子體輸運的邊界條件模型，所計算的伏安特性曲線與文獻(xiàn)值較為吻合，說明本文的改進(jìn)模型是合理的。