本文“热流模型”部分的内容基本源于Warner N Z博士的博士论文Theoretical and experimental investigation of Hall thruster miniaturization[1]. 本文是对其的一份中文翻译与再整理，原作者保留所有知识权利。而后的有限元仿真部分则是笔者自行运算结果。
摘要
霍尔推进器是一种高比冲、高推力密度的小型电推进器装置。工作时由于粒子碰撞和粒子-壁面碰撞而产热。实际工况中，为防止铁磁材料（如磁芯）温度超过居里点（Curie Temperature）而发生退磁，导致磁场位形畸变甚至失效，需要对其产热与传热机制进行详细建模与仿真分析。使用理论公式计算热通量，并作为条件与约束输入有限元仿真程序，可得到高可信的温度分布结果，可为实际设计与验证提供有效理论参考。
能量守恒关系
稳定状态下的功率平衡关系可以表述为以下关系[1]：
$$
P_d=P_{exhaust}+P_{radiation}+P_{conduction}
$$
霍尔推进器主要能量输入为阳极与阴极之间的放电电压及线圈绕组电流，在此统一合并为$P_d$。能量损失的途径主要包括热损耗、热辐射与热传导。热损耗主要发生在下游羽流区域，是指电子、离子与中性粒子合为等离子体相态向外喷出；热辐射主要发生在外壁面向宇宙背景空间之间，是真空中霍尔推进器主要散热途径之一，但同时，碰撞区的高温粒子也会以辐射方式向通道壁面加热，造成温升；热传导发生在装配元件之间，实际中可通过散热/隔热板等方式进行合理调控。
能量损失可划分为以下几种途径，分别为：Ions Ejected in Plume, Ions Lost to Walls, Electrons Ejected in Plume, Electrons Ejected to Walls, Electrons Lost to Anode. 此外，辐射散热和接触热阻也被考虑在本模型中。接下来，将按照分类分别简要讨论这几种途径的热流模型公式。
热流模型
Ions Ejected in Plume
离子主要出现在碰撞区域与羽流区域；其能量可通过机械运动动能、热力学动能及电离能描述：
$$
P_{ip} = [\frac{1}{2}m_ic^2+\frac{5}{2}kT_{ip}+eV_i]\frac{I_b}{e}
$$
其中，离子速率c通过目标比冲与理想缩放模型计算得到，或者通过加速效率和工作电压计算得到。羽流中的离子温度$T_{ip}$通常被估计为1eV（电子伏特，1 eV ≈ 1.602×10⁻¹⁹ J，等离子体物理中常用的能量单位）或者通过其他实验数据得到。如果推进工质是Xe，则其电离能$V_i=12.1eV$。羽流中的束流电流通常被估计为放电电流的70%，即$I_b \approx 0.7I_d$。有文献提供了束流电流的理论公式：
$$
I_b=\frac{\dot{m}_a e}{m_i}\eta_u
$$
其中，$\dot{m}_a$为质量流量，更多详细信息可参考文献[2]。
Ions Lost to Walls
离子在壁面的复合现象是通道壁面温升的最直接因素之一。离子被鞘层加速冲向壁面，最终将能量传递到壁面造成加热。能量总和包含：形成离子所需的电离能、离子自身热能、预鞘层中所获动能、鞘层电场传递给离子的势能：
$$
P_{iw}=[\frac{1}{2}m_iv_B^2+\frac{5}{2}kT_i+eV_i+e(\phi_s-\phi_w)]\frac{I_w}{e}
$$
文献[3]研究了径向（radial）、周向（azimuthal）方向的离子温度，大约在个位数电子伏特数量级。离子定向运动的能量转化关系通过鞘层效应得以体现，因此径向和周向离子温度是公式中$T_i$的主要输入项。
鞘层及其内部的作用机制应当被更详细讨论。等离子体中的电子由于质量极小，会更容易在壁面发生沉积，形成负的电场，该电场则会在接下来排斥后续电子而反过来吸引正离子。这个排斥-吸引作用在壁面附近形成了一个电势急剧下降的区域，称为鞘层。最终，鞘层内的电场强度达到平衡，壁面电位稳定，鞘层厚度不再变化。
鞘层向离子传递的势能可通过鞘层电位和壁面电位之间的差值进行估算。对于具有薄平面鞘层并且离子温度远小于电子温度的等离子体系统，其鞘层电位可表示为：
$$
\phi_s=-\frac{kT_e}{2e}
$$
经过近似假设，壁面电位表示为：
$$
\phi_w=\frac{kT_e}{2e} \left [ \ln(2\pi\frac{m_e}{m_i})-1\right ]
$$
文献[3]指出，对于一个300V放电电压的霍尔推进器，放电区域的电子温度范围在25-30eV。因此，对于一个25eV的电子温度输入，传递给离子的势能\phi_s-$\phi_w$可计算得出为132eV。
离子向壁面的电流$I_w$可表示为束流电流$I_b$的某个比例。当离子与壁面接触时，会与电子发生复合反应，随后以中性原子（如Xe原子）的形式反弹。中性原子可能会被再次电力，也可能从推进器中逸出。对于从通道中加速出来的束流离子，在再次电离并进入尾流之前，已经在壁面经历了一次符合过程。在该过程中，
$$
I_w=\alpha_wI_b
$$
对于$\alpha_w$的计算，可参考文献[3]。在本文后续仿真所依赖的输入数据中，该参数被估算为50%，即束流中每两个离子中就有一个经历了复合功率损失。这意味着有50%的束流离子最终与通道壁发生相互作用并沉积能量。
Electrons Ejected in Plume
用于中和羽流的电子来自阴极，通过阴极与临近等离子体中的微小电位差而被推入束流中。该电位差表示为$\phi_{nc}$，典型值在15-20V，赋予了电子初始动能。电子同样携带热能，可通过羽流中的电子温度$T_{ep}$计算，该温度可在稳定推进工况中通过实验手段得到，典型值在3eV或更低[4]。从而，电子损失到羽流中的总能量$P_{ep}$可表述为：
$$
P_{ep}=[\frac{5}{2}kT_{ep}+e\phi_{nc}]\frac{I_b}{e}
$$
Electrons Ejected to Walls
与壁面上发展完成的鞘层接触的平衡等离子体会接收$2kT_e$的平均能量。这个结果来自对球坐标系下修正麦克斯韦分布函数的积分：
$$
E_{ew}=\int\limits_0^\infty\int\limits_0^{\frac{\pi}{2}}\int\limits_0^{2\pi}\tilde{f}(w)(w\cos\theta)(\frac{1}{2}m_ew^2)w^2\sin\theta \mathrm{d}\chi\mathrm{d}\theta\mathrm{d}w
$$
被排斥的电子具有运动速率$w$，$f(w)$描述了其在三维球坐标系下的速率分布。电子速率的修正分布假设在远离壁面的区域，电子服从标准麦克斯韦分布，该修正分布考虑了减速势$\varphi$对电子密度的影响：
$$
\tilde{f}(w)=n_\infty \left ( \frac{m_e}{2\pi k T_e} \right ) ^\frac{3}{2} \exp(-\frac{1}{kT_e}\left(\frac{m_ew^2}{2}+e\phi\right))
$$
上式$E_{ew}$的积分结果为：
$$
E_{ew}=\frac{1}{4}n_\infty\bar{c}_e e^{-\frac{e\phi}{kT_e}}(2kT_e)
$$
括号前的系数用于确定到达壁面的电子通量。括号内的$2kT_e$项表示每个到达电子在壁面处获得的平均能量。该数值可用于估算推进器壁面因电子碰撞而接收的总功率$P_{ew}$，其中，电子温度与壁电流的估算方法如前所述：
$$
P_{ew}=[2kT_e]\frac{I_w}{e}
$$
Electrons Lost to Anode
阳极表面同样形成了鞘层。扩散至阳极的电子携带少量热能，与阳极鞘层发生排斥作用。文献[5]表明近阳极鞘层电势$\phi_{na}$约为4V。该过程与Electrons Ejected to Walls作用机制类似，故阳极吸收的能量功率可表示为：
$$
P_{ea}=[2kT_{ea}]\frac{I_d}{e}
$$
阳极附近的电子温度（$T_{ea}$）通常低于电离区域的温度。在电离过程中，电子释放能量后会冷却，并通过与中性粒子碰撞向阳极方向散射，具体数值取决于电离区域的结构特征，详细讨论可见文献[3]。
Excitation and Radiation
中性粒子可通过碰撞过程被激发到高能级状态，当激发态粒子返回基态时，会释放多于能量。推进器表面辐射和吸收的热量将导致推进器发射，这种辐射功率为：
$$
P_{rad}=[\psi e V_i]\frac{I_b+I_w}{e}
$$
系数psi是推进器内每个粒子产生的激发能与电离能之比。当氙气电子温度高于其电离能（12.1eV）时，该系数约为2.2，且与温度参量没有显著相关性。
考察激发与辐射需要深入探究粒子能级跃迁、粒子碰撞的详细作用机制。本文暂只对辐射分配情况作简要估算：若激发主要发生在推进器出口附近的高密度区域，辐射能量将半数逸出推进器，另一半被壁面吸收；若激发发生在环形空间更深层区域，则可采用壁面面积与推进器出口面积的比值进行估算。综合这两种方法，能量占比通常介于0.5-0.8之间。为简化并保守估计，直接假设50%辐射功率作用于推进器内外壁面，25%功率被阳极表面吸收。
Thermal Contact Resistance
霍尔推进器装配过程中，两种接触材料之间并非完全理想接触。表面光洁度与清洁度、材料特性、接触压力、间隙流体（背景气体）及其压力、接触表面温度等，将会对材料间的热量传递产生影响。即使在真空环境中，由于装配压力，我们认为这些微小空隙内的背景气体将不逸出，因而遵循对流传热规律。
对接触热阻的简化建模可以在界面处引入一个虚拟薄层，通过定义热阻/（热导）的方式引发温度梯度差异。对于总面积为$A$且温度差异为$\Delta T$的接触区域，通过界面层的热传递表示为：
$$
Q = A h_c \Delta T = \frac{\Delta T}{R_t}
$$
式中，$R_t$定义为接触热阻，具有单位[K/W]. 在COMSOL软件中，接触热阻可直接通过固体传热接口的“热接触”边界条件——“等效薄热阻层”的方式定义，参数为层热导流率或者层热阻，辐射热导率等。因此，本文不采用原文中赋予该层实际厚度和材料属性的方法，而是直接使用COMSOL的“等效薄热阻层”功能进行定义。
在实际仿真分析中，需对不同材料间的层流热导率作仔细考察，或者实验验证，以准确模拟接触热阻效果。
有限元仿真分析
程序构建
按照公开资料，可在COMSOL软件中直接构建如上图所示的霍尔推进器的几何模型。包括以下主要部件：阳极、通道壁、磁屏蔽、线圈、磁芯、钢支撑结构。在实际研究中，可根据CAD模型或其他3D模型直接导入，以实现对真实工件的高精度仿真。
对材料的定义至关重要。为了简化，本仿真尽量采用了COMSOL内置材料，包括线圈的铜材料（Copper）、铁芯的软铁材料（Soft Iron (Without Losses)）、阳极石墨材料（Graphite）、钢支撑材料（Structural steel）. 本文独立定义了通道壁氮化硼材料的材料属性，相关性能参数可通过公开商业资料查询。
物理场（边界条件、初始条件）的定义是仿真的核心。本文使用固体传热（ht）接口，通过“热接触”边界条件模拟壁面间的接触热阻，通过“热源”条件模拟线圈欧姆发热，通过“热通量”边界条件模拟离子、电子分别向壁面、阳极的热沉积，以及粒子辐射对壁面和阳极的加热。使用“表面对环境辐射”模拟霍尔推进器对环境散热，这是在真空环境中霍尔推进器的唯一散热方式。
使用“物理场控制网格”对模型的计算网格进行构建，单元大小为“细化”。实际研究中，一方面可以先选择较粗网格以快速得到结果并建议物理正确性，另一方面可根据需求与模型复杂度精细化定义网格分布。后者则是研究工作中的常见手段。
结果与讨论
对物理场进行计算，可得温度分布结果如下：
由温度分布图可以看出，推进器的主要发热位置出现在线圈及通道壁面。本文参考了文献[1]中的热沉积定义——均匀施加在通道壁面，因而温度梯度并未明显体现在通道壁面的轴向上。实际工况下，位于下游通道口附近的电子约束区（碰撞区）的粒子密度与发生碰撞的频度远远高于上游漂移区域，因而壁面轴向应会在下游通道口出现最大温度，并向上游方向呈现明显递减的温度梯度。后续可考虑函数定义等方式改善热沉积分布。
此外，外部钢支撑的温度明显低于内部组件的发热区域。两个部分的主要接触都发生在阳极出口附近的连接上，而对接触热阻的定义也全部发生在这里。这表明接触热阻的引入有效对热量传递产生影响，得到了与其他文献相一致的温度分布结果[1,4]。
可以对接触热阻带来的温度梯度做进一步考察：定义一维截线，考察由钢支撑到通道壁的温度变化情况，如下图右侧红色点划线所示。
绘制这个截线上的温度数值。由图可见，在连接处，温度从 566.682°C 跃升至 646.869°C，并且并不连续。这其实与先前的接触热阻定义方式直接关联。当规定了“等效薄层热阻”的层热导流率，软件在解算傅里叶导热定律$\vec{Q} = -k \nabla T$时，在稳态下热流密度Q需连续，而k在界面突变，必然导致温度梯度调整，最终表现为温度T在界面处不连续。
参考文献
[1] Warner N Z. Theoretical and experimental investigation of Hall thruster miniaturization[D]. Massachusetts Institute of Technology, 2007.
[2] Goebel D M, Katz I, Mikellides I G. Fundamentals of electric propulsion[M]. John Wiley & Sons, 2023.
[3] Szabo, Joseph J .Fully kinetic numerical modeling of a plasma thruster[J].ph.d.thesis massachusetts inst.of technology, 2001.
[4] Wang Y M, Xu S, Levchenko I, et al. Approach to a simplified numerical optimization of low-power Hall thrusters[J]. Vacuum, 2018, 152: 173-183.