热离子能量转换器铯电弧工况下电子势能的分布特性

钟武烨，吕 征，郑剑平，杨启法

(中国原子能科学研究院 反应堆工程技术研究部，北京 102413)

热离子反应堆电源是采用热离子能量转换方式将反应堆裂变释放的热能直接转换成电能的能量转换装置，其具有比功率高、无机械转动、所需比散热面积小等优点，是国际上最具优势的高效空间核电源之一[1]。目前，具有实际应用价值的热离子能量转换器的结构形式是充铯的工作在电弧工况的热离子二极管[2]，其通常采用难熔金属作电极对，发射极被加热到高温(～1 500 ℃)产生热电子发射，充入电极间隙的铯蒸气在其中电离形成等离子体，电子输运通过等离子体区域后达到接收极，经外电路的负载做功后重新回到发射极，形成回路。热离子转换器可看作以电子作为介质的热机[3]，电极间隙的势能分布则决定该热机的能量来源。因此，电子在电极间隙的势能分布是研究热电转换特性的关键。

目前，对关于电极间隙的电子势能分布已有较深入的了解，但一些基础问题仍有待进一步深入研究，如热离子转换器运行在电弧工况不同区域(包括饱和区、阻塞区)时，对于电极鞘层的电子势能分布的认识仍不够深入，忽略对阻塞区虚拟阴极的讨论[4]；在讨论阻塞区虚拟阴极时又无法给出数值上可解的方法[5]；缺乏对接收极电势分布特征的研究等。导致该问题的原因主要是现有的模型并没有很好地处理等离子体输运及其边界条件的匹配问题。本文基于动理学方法建立电弧工况下的等离子体输运模型，对等离子体鞘层的电子势能分布进行分析，并将等离子体鞘层处理为输运的边界条件，自编程实现等离子体输运方程与边界条件的自适应解耦求解。

1 理论建模

电子势能分布的研究需建立等离子体的输运模型。当热离子转换器工作在电弧工况(图1)时，在电极间隙中，等离子体与电极表面相接触的区域形成电势突然跃变的、非准中性的等离子体鞘层，其尺度与德拜长度相当，远小于电极间隙。因此，电极间隙的等离子体可分为电极附近的鞘层和中间的准中性等离子体区[6]。在具体的建模过程中，在准中性等离子体区建立带电粒子的输运模型，而等离子体鞘层则处理为输运的边界条件。

1.1 电弧工况等离子体输运模型

热离子能量转换器电极间隙的等离子体属于低温弱电离的等离子体，此时必须考虑电子、离子以及中性原子3组分同时存在的多粒子碰撞输运问题。等离子体中带电粒子体系间存在复杂的相互作用，要精确描述等离子体的行为是极其困难的。目前，只能根据不同的条件和研究的问题，采用不同的近似方法[7]。文献给出的建模方法可归为唯象模型和基于介观尺度的动理学模型两大类[8]。唯象模型作为定性分析是有效的，作为定量分析计算仍有差距。本文的输运方程采用动理学的模型。动理学模型是输运过程严格的定量理论。动理学方法研究其输运过程的基本思路是从粒子碰撞迁移所满足的玻尔兹曼方程出发，然后根据转换器等离子体输运的具体特点对分布函数和碰撞输运方程作合理简化，建立输运控制方程组，并给出输运系数。

由于电极间隙远小于电极的轴向尺寸，故只需建立电极间隙准中性等离子体输运的一维模型，在建模过程中以发射极作为原点建立坐标系。中性铯原子按理想气体处理：靠近电极的中性原子温度等于电极温度，电极间隙各点的温度随电极温度呈线性变化，不必考虑中性原子的输运问题，但必须考虑由于中性原子的存在导致的散射效应。等离子体输运的描述的目的是获得电极间隙等离子体浓度n、电子温度Te、负电势V、电子电流密度je及其能量流密度Se、离子电流密度ji及其能量流密度Si等分布状态。

图1 热离子能量转换器的工况分类及工作原理Fig.1 Operation mode and working principle of thermionic energy converter

首先，建立电子、离子、中性原子3组分共存条件下的玻尔兹曼方程：

式中：β=e,i,a(e、i、a分别表示电子、离子、中性原子)；α=e,i；v为速度；r为位置；e为电量；下标in表示非弹性碰撞；Jαβ为碰撞积分项；fα为分布函数并假定其偏离平衡态不远。通过对玻尔兹曼方程的简化，可建立等离子体参量的通量及连续性方程[9-10]。

1) 通量方程。带电粒子的通量包括引起带电粒子输运的浓度梯度、温度梯度以及弱电场力的贡献，电子、离子的电流密度分别为：

为考虑离子电流密度的贡献，总电流密度jtot包括电子电流密度和离子电流密度的部分。

jtot=je-ji

对于带电粒子的能量通量，考虑了宏观运动的动能、势能以及热扩散的贡献。

式中：k为玻尔兹曼常数；De,i、μe,i、ke,i、βe,i分别为浓度扩散系数、迁移率、热导率、热扩散比；Rie为离子受电子的作用力。

2) 连续性方程。带电粒子通量的连续性是由体电离效应决定的。

式中：Γ为电离函数；σ0(Te)为与电子温度相关的等效电离截面；Na为原子浓度；n(Te)为热平衡态下的等离子体浓度。

对于能量通量的连续性，重粒子在电离/复合过程中吸收/放出的能量由电子提供/带走，忽略电子的辐射ΔSrad以及e-i、e-a的碰撞损失ΔSei、ΔSea，则有：

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

式中，Eion为铯离子的第一电离能。

上述的通量及连续性方程组成了等离子体输运的控制方程组，其中包含了je、ji、Se、Si、n、V、Te7个变量。对于其中的输运系数，当分布函数偏离平衡态不远时，基于动理学理论可对分布函数作Sonine多项式展开，从而将输运系数分级展开，取其中的前3项即可使截断误差小于10%[11]。

1.2 电弧工况等离子体输运的边界条件

图2 电弧工况饱和区可能的电子势能分布Fig.2 Possible distribution of electron potential energy at saturation region of ignited mode

为获得等离子体输运微分方程的边界条件，需建立带电粒子通量、能量通量的边界值与等离子体边界参量的关系，其中的关键点是等离子体鞘层的电子势能分布。在电弧工况的不同区域，等离子体鞘层电子势能的分布状态是有差异的，所有可能的势能分布状态如图2、3所示，其中VE,C分别为发射极、接收极鞘层电势的跃变，φE,C分别为发射极、接收极功函数，Ud、U分别为电弧压降及输出电压。对于发射极鞘层，其分布在饱和区是单调下降的(图2)；而在阻塞区则是先增后降(图3)，形成了电子发射的势垒ΔVE。另一方面，对于接收极鞘层，其电势分布并不依赖于具体的工作区域，其正负性取决于边界条件本身的匹配。以下依据电极鞘层内电子势能的分类对边界条件进行讨论，并将接收极表面的电势定义为零电势参考点。

图3 电弧工况阻塞区可能的电子势能分布Fig.3 Possible distribution of electron potential energy at obstructed region of ignited mode

1) 发射极鞘层的边界条件

电弧工况发射极鞘层的边界条件需区分为饱和区与阻塞区两种情况进行讨论。

在饱和区，发射极附近的电子势能是单调下降且跃变较大，由此产生显著的附加电场促进发射，即产生肖脱基效应，从而获得较大的输出电流密度[12]。边界通量需考虑带电粒子的热发射以及来自于等离子体的反射。

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

式中：jeE、jiE、SeE为发射极鞘层与等离子体交界面上的通量；TeE为交界面上的电子温度；nE为交界面上的等离子体浓度；jSE(EE)为热电子的肖脱基发射项；r1、r2为考虑了发射的动理学反射特性而引入的参数[13]。

在阻塞区，发射极鞘层的电子势能是先增后降的，即热电子发射的势垒增大了ΔVE，形成虚拟阴极，由于未促进发射附加电场，肖脱基效应也就不存在[14]，此时的发射电流密度是饱和发射电流密度jS0。

2kTeE+je(Ud+ΔVE)

2) 接收极鞘层的边界条件

当接收极的电势跃变VC>0时(对应于图2a与图3a)，有：

式中：jeC、jiC、SeC为接收极鞘层与等离子体交界面上的通量；f0与f1为描述接收极鞘层内电子运动的各向异性而引入的函数[15]。

当接收极的电势跃变VC<0时(对应于图2b与图3b)，有：

2 模型求解与算法流程

根据上述建立的输运方程和采用的边界条件，输运系数依赖于等离子体参量本身，因此输运方程本身是非线性的，故本文建立的模型只能数值求解。另一方面，输运的边界条件依赖于工作参数所处的工况区域(饱和区或阻塞区)，但实际计算中并不能根据输入条件(工作参数)预知其工况区域。因此，在数值求解的过程中，应使得边界条件与输运方程自适应匹配求解。本文采用的算法流程如图4所示。

1) 输入计算所必要的参数：基本物理参量、热离子转换器的运行参数，如电极温度(TE、TC)、铯温TCs、电极的吸铯功函数(φE、φC)、电极间隙d、总电流密度(输出电流密度)jtot等。

图4 热离子转换器电弧工况等离子体输运算法流程Fig.4 Algorithm flow of plasma transportation of TEC at ignited mode

2) 给定接收极鞘层与准中性等离子体区交界处的电子温度TeC和接收极鞘层的电势跃变VC作为迭代初始值，根据VC的正负性选择相应的接收极边界条件模型，计算此边界的jeC、jiC、nC及SeC、SiC。

3) 以流程2计算得到的参量为边界值，采用龙格-库塔方法对输运微分方程进行积分，得到准中性等离子体区Te、n、je、ji、V、Se、Si的分布，特别地，记发射极边界上的电子电流密度积分为je_inte(TeC,VC)、离子电流密度积分为ji_inte(TeC,VC)。

4) 先假定工作点处于阻塞区，根据阻塞区发射极边界条件计算非单调的电势跃变情况VE及ΔVE。若解存在，则认定虚拟阴极的存在(即处于阻塞区)；若解不存在，则认定处于饱和区(即虚拟阴极不存在，发射极电势跃变是单调的)，并据此计算相应的电势跃变VE。

5) 根据边界条件的自适应判断，计算相应工况下发射极边界上的电子电流密度je_cal(TeC,VC)、离子电流密度ji_cal(TeC,VC)，并与积分所得边界值联立成方程组je_cal(TeC,VC)-je_inte(TeC,VC)=0与ji_cal(TeC,VC)-ji_inte(TeC,VC)=0。采用牛顿迭代法求解方程组，直至收敛到设定的精度。若发散，则重新选择初始条件进行计算。

6) 输出计算结果。

3 结果讨论

由于程序的输入包含电极温度、电极功函数、铯温、电极间隙、总电流密度(输出电流密度)等7个参数，在结果讨论中不适宜采用控制变量的方式。Bullis和Rasor[16]在其综述中指出，尽管热离子转换器的工作参数有差异，其输出特征可采用离子裕度β及参量pCsd来区分。其中，离子裕度β由发射极的发射特性决定，用来表征热发射离子中和热发射电子的程度，其临界值为1；参量pCsd由输运过程决定，用来表征电子输运过程的散射强度，其临界值为20 mile·Torr(1 Torr≈133.322 Pa，1 mil=25.4 μm。

为在较宽的参数范围内验证计算模型的准确性，受限于电极材料，本文研究的热离子转换器是电极对在给定的电极温度下运行，仅在可变的铯温下进行输出，此时离子裕度β被限定而参量pCsd可在较宽的范围取值。为便于与文献的数据进行比较[17]，程序计算中输入的工作参数列于表1，选定的电极对材料为Re-Mo，限定电极的工作温度TE=1 700 K、TC=925 K，电极间隙d=0.25 mm。根据表1，由于电极特性的制约，离子裕度β限制在小于1的区域，而参量pCsd则覆盖了较宽的范围。将上述工作参数作为程序的输入，获得本文所有的计算结果。

表1 程序计算中输入的工作参数Table 1 Input working parameter for procedure

1) 计算结果的有效性判定

首先，为判定计算结果的合理性，需将等离子体的德拜长度rD(表2)与特征尺寸(电极间隙)进行比较，德拜长度均远小于电极间隙，说明电极间隙形成了等离子体，所计算的工况都处于电弧工况；并且随着总电流密度的增大，工况状态逐渐从电弧工况的阻塞区进入饱和区，在饱和区发射极鞘层电场值满足EE远大于0。因此，计算结果验证了边界条件分区设置的合理性。

表2 计算结果有效性判定所需的数据Table 2 Datum needed to judge validity of computing result

2) 伏安特性

图5示出本文计算得到的伏安特性及其与文献值的对比，计算值与实验值符合较好。图6比较了TCs=551 K与591 K时的伏安特性，计算表明，当参量pCsd<20 mile·Torr时，饱和区的输出电流密度表现出较明显的饱和特征；而当参量pCsd>20 mile·Torr时，饱和电流密度的特征不明显。

图5 本文计算的伏安特性曲线与文献值的对比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) 电子势能分布

当热离子转换器的铯温为551 K时，电极间隙电子势能的变化如图7所示，涵盖了电子从发射极费米面热发射输运到接收极费米面的全过程。整体而言，电子势能在发射极鞘层的跃变大于接收极鞘层的跃变，而在准中性等离子体区则相对平缓(其中跃变发生的尺度仅有德拜长度的量级，为便于展示，图片对跃变的尺度进行放大)。

对于特定的工作点(例如jtot=6 A/cm2)，电子始末位置的能级差，即两电极费米面的能级差，构成了对外输出电压。两电极表面势能的差用于维持电弧状态，属于固有的内部压降。

图7 TCs=551 K时电子势能随电极间位置的变化Fig.7 Dependence of electron potential energy on location of interelectrode at TCs=551 K

对于发射极鞘层的电势跃变，当电流密度较小时(jtot=1 A/cm2)，发射极鞘层的电子势能是先增后减的非单调跃变，在发射极的鞘层内形成了附加的发射势垒(即所谓的虚拟阴极)，相当于增加了发射极功函数，从而抑制了电子发射，但发射极功函数的提升有利于获得更大的电压输出，此时转换器工作在阻塞区。随着电流密度的增大，发射极鞘层的电子势能变为单调向下跃变，其量值随电流密度而增大，此势能跃变使电子热发射的过程中得到电场的激励(即肖脱基效应)，同时该电场较好地阻碍等离子体中的电子返回发射极；为维持大电流密度，就必然需要更高的电弧压降(内压增大)，从而输出电压也就降低了，此时转换器工作在饱和区。

图8 TCs=551 K时发射极鞘层电子电流密度的构成Fig.8 Consistency of electron current density in sheath of emitter at TCs=551 K

整体而言，在输出电流密度递增的过程中，发射极鞘层的电子势能分布经历了由非单调变成单调的过程。对发射极鞘层电子电流密度的构成分析可进一步说明产生此分布特征的机理，如图8所示，在小电流密度的阻塞区，从等离子体返回的电子电流密度超过净电子电流密度，相当于发射极鞘层“累积”了电子，这将增加电子发射的势垒，形成了虚拟阴极；在大电流密度的饱和区，热发射的电子电流密度远大于从等离子体返回的电子电流密度，此时在发射极附近产生强烈的电离(图9a)，发射极附近产生大量离子，形成了促进电子发射的单调向下的电子势能分布。

对于接收极鞘层的电势跃变，在输出电流密度递增的过程中，其电子势能分布经历了从向下跃变(VC<0)到向上跃变(VC>0)的过程。产生这种分布的原因是在输出电流密度从0开始增大的过程中，实验观察到最先发生电离的区域是接收极的鞘层[18]，由此从扩散工况进入饱和工况的阻塞区，因此需要向下跃变的电子势能分布促使接收极附近发生电离(图9b)。而随着输出电流密度的增大，电离的峰值逐渐向发射极移动，接收极一侧趋向于电离平衡(图9a)，此时接收极鞘层表现为准中性等离子体鞘层的性质，即热速率较大的电子运动到接收极表面逸走而留下富余的正离子，从而使电子势能向上跃变。

以上关于电极间隙电子势能的分布特性，在转换器的铯温为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时电弧工况不同工作点下对应的电子势能分布Fig.10 Dependence of electron potential energy at different working points of ignited mode at TCs=591 K

4) 根据输出功率特性对电子势能分布特性进行确认

为确认上述的热离子转换器电极间隙电子势能分布特性，进行最大输出功率的讨论。一般而言，当电池的内外压降相等时有最大的电功率输出，这是电池的一基本属性。根据电子势能分布，可获得输出电压、电弧压降与电流密度的关系，如图11所示，电流密度越大，输出电压越小，而电弧压降则越大，两者的交点为内外压相等的工作点。同时，根据计算的伏安特性，可获得输出电功率密度与电流密度的关系(图12)，当转换器铯温为551、591 K时，在对应内外压相等的工作点jtot=5、14 A/cm2下获得最大的输出电功率。这说明了本文关于热离子转换器的计算结果是自洽的，并与一般电池的功率特性相一致。

图11 输出电压、电弧压降与输出电流密度的关系Fig.11 Dependence of output voltage and plasma voltage drop on output current density

图12 输出电功率密度与电流密度的关系Fig.12 Dependence of output power density on output current density

4 结论

本文基于动理学方法建立了描述热离子转换器电极间隙电弧工况等离子体输运的方程及相应的边界条件，并设计了一套基于牛顿迭代法的计算程序实现其自适应解耦求解，详细讨论了等离子体的电子势能分布特性。

1) 电子势能在发射极鞘层的跃变大于接收极鞘层的跃变，而在准中性等离子体区则相对平缓。随着输出电流密度的增大，发射极鞘层的电子势能跃变特征将从先增后减的非单调变为单调下降，这是由发射极鞘层内电子电流密度构成的差异导致的。

2) 同时接收极鞘层的电子势能则从向下跃变变为向上跃变，这是由接收极鞘层附近的电离状态决定的。

3) 基于等离子体鞘层电子势能分布的分析，改进了等离子体输运的边界条件模型，所计算的伏安特性曲线与文献值较为吻合，说明本文的改进模型是合理的。