脉冲等离子体推力器放电电离特性的数值研究

杨　磊，刘向阳， 汤海滨，黄玉平，王宁飞

(1.北京精密机电控制设备研究所，北京　100076；2.北京理工大学 宇航学院，北京　100081；3.北京航空航天大学 宇航学院，北京　100191)



脉冲等离子体推力器放电电离特性的数值研究

杨磊1,3，刘向阳2， 汤海滨3，黄玉平1，王宁飞2

(1.北京精密机电控制设备研究所，北京100076；2.北京理工大学 宇航学院，北京100081；3.北京航空航天大学 宇航学院，北京100191)

推进剂烧蚀产物的电离过程直接制约着等离子体的产生及其特性，是决定脉冲等离子体推力器(PPT)工作效能的关键因素。为了准确描述PPT放电电离特性，提出了基于热化反应的烧蚀边界层和等离子层耦合模型，并以LES-6 PPT为例进行数值研究，揭示了脉冲电流激励下PTFE烧蚀与离子化过程所产生的多种组分在放电室内时空分布特性。结果表明，烧蚀并电离出的等离子体主要以一阶电离形式存在，其中C、F元素的离子化区域分别集中在上、下游放电通道内；高度电离且短期存在的C、F离子(二价和三价)主要出现在放电初期；等离子体各种组分浓度的峰值时刻与羽流实验数据基本吻合。该研究为探究PPT放电能量分配及损失机理提供了理论参考。

电推进；脉冲等离子体推力器；烧蚀；电离

0　引言

固体烧蚀型脉冲等离子体推力器(Pulsed Plasma Thruster，PPT) 由于具有比冲高、结构简单、平均功率小和控制灵活等特点，可应用于轨道转移、位置保持、阻力补偿、精确编队飞行及姿态控制等许多微型航天器的推进任务，成为当前电推进技术研究的热点和重要方向之一[1]。然而，长期以来PPT效率低下一直是限制其发展和空间应用的核心问题[2]。试验研究发现，其原因在于PPT工作时所消耗的推进剂只有小部分在微秒级的放电烧蚀过程中被电离和加速喷出产生推力[3]。因此，从推进剂烧蚀演化机制入手，研究PPT放电电离过程，定量揭示推进剂烧蚀组分尤其是离子组分的时空分布特性，这对探究放电能量分配及损失机理有着至关重要的作用，同时也是研制高性能PPT亟需解决的关键。

合理完善的烧蚀分析模型是准确描述PPT放电电离特性,揭示等离子体产生、流动、扩散等过程的基础。Mikellides通过壁面与等离子体层之间的热力学状态关系首次将烧蚀与等离子体流动结合起来描述PPT放电烧蚀过程[4]。此后，基于MACH2磁流体程序的一系列PPT数值研究[5-6]都是建立在Mikellides模型的基础上。然而该模型的边界层跳跃条件采用了简单的单值比例关系，不能体现放电过程中等离子体流动入口参数的动态变化[7]。针对这一缺点，Keidar等将PPT的放电室按空间位置划分为烧蚀边界层和等离子体层两个区域，通过克努森理论来刻画烧蚀和流动耦合机制，建立更为精准的K-B模型[8]。尹乐等在K-B模型的基础上开展了三维双温PPT工作过程的仿真研究[9]。然而由Anisimov一维平衡蒸发理论发展的K-B模型忽略了推进剂的热化反应对烧蚀边界层和等离子体层相互关联产生影响[10]。固体推进剂的热化反应最直观体现在电离组分的复杂多样性。受热化反应控制的这些组分状态变化很大，直接影响着烧蚀边界层与等离子体层间传热传质密切相关的关键参数，对PPT工作过程准确模拟的作用不容忽视。

为此，本文提出了基于热化反应的烧蚀边界层和等离子体层耦合模型，并以LES-6 PPT为例开展了PPT放电电离特性的数值研究，定量获取微秒脉冲电流激励下PTFE烧蚀和电离的各价C、F离子组分及其浓度存在时空分布规律及中性原子在脉冲放电末期及结束后的演化机制，初步揭示了推力器性能贡献来源及其效率影响机制。

1　计算模型及方法

针对传统模型的缺陷，本文以PPT常用的固体推进剂PTFE为例，首先建立其热化反应模型来获取PTFE烧蚀与离子化过程所产生的组分信息；其次构建能够表征热化反应效应的烧蚀边界层和等离子体层的耦合模型，用以描述定常条件下的烧蚀与流动机制；该模型考虑了传统模型所忽略的PTFE热化反应所造成的烧蚀多组分及其质量分布特性对边界层热力学状态尤其是外边界速度的影响机理。最后在此基础上结合放电模型及磁流体动力学方程定量化地揭示微秒非稳态工作过程中PPT放电室空间的电离特性。

1.1热化反应模型

本文重点研究等离子体层PTFE高温等离子体系组分及其热物性，主要包括一价至三价C、F离子，C、F原子及电子等。

在局部热力学平衡的假设基础上，PTFE高温等离子体系组分可满足如下的Saha方程[11]：

(1)

式中ne为电子数密度；n(r)为r次电离粒子的数密度；me为电子质量；εi(r)为r次到r+1次电离所需电离能；Z(r)为r次电离的配分函数(r=0，1，2)。

为了计算上述Saha方程,需要进一步补充以下假设和方程：

(1)等离子体宏观电中性假设

(2)

式中ni为第i种等离子组分的数密度；r为电离阶次。

(2)数密度守恒

根据固体推进剂PTFE的组成，即(C2F4)n，可得

(3)

(3)气体分压定律

(4)

式中pir=nirkT表示第i种成分的分压。

由Saha方程及式(2)、式(3)、式(4)即可确定PTFE烧蚀与离子化产生的多种离子组分、原子浓度、电离程度等组分信息与局部热力学状态(温度和压强等)间的关系。根据确立的组分及浓度信息，进一步获取建模中所需的物性参数变化规律，如相对分子质量、电导率、比焓等。

1.2烧蚀边界层和等离子体层耦合模型

PPT的放电室按空间位置所划分的烧蚀边界层和等离子体层如图1所示。其中，烧蚀边界层(也称为克努森层、蒸气层或动力层)用于描述PTFE壁面(0处)与流体入口(1处)间烧蚀过程的过渡层，并且烧蚀边界层的外边界参数可以表征固体推进剂的烧蚀特性；等离子体层则描述了等离子体的产生及其流动过程，反映了PPT放电电离特性。

本文基于克努森理论改进了K-B模型，建立了适用于PPT烧蚀边界层特点的分析模型:

(5)

式中di=mi/2kTi(i=0,1)；m0和m1分别为壁面和流体入口粒子的平均分子质量；α=V1/(2kT1/m)0.5；k为Boltzmann常数；erf为误差函数；erfc为余误差函数，且有erf(α)=1-erfc(α)；β为Anisimov比例系数；γ=(j+5)/(j+3)为绝热指数；j为内部自由度，体现PTFE热化反应所产生的组分结构的影响。

图1　PPT烧蚀与等离子体层结构示意图

式(5)描述了边界0和1处热力学状态参数的本质关系，即若已知PTFE壁面温度T0、外边界速度V1即可由式(5)解算出等离子体层入口的热力学状态参数如温度T1、密度n等。其中，T0可由式(6)所示的传热模型，即半无限大固体材料非稳态传热方程、PTFE汽化分解条件及传递到推进剂表面净热流Qs所确立的边界条件获取[12]。

(6)

式中Q=I2Rp/(h·w)为单位面积的等离子体焦耳热；I为放电电流；Rp为等离子体电阻，可通过热化反应模型所提供的电导率来获取；h和w分别为烧蚀壁面的高和宽；Δhj为等离子体在计算时刻一定温度、压强下比焓值与在初温T0时比焓值之差，而比焓值可相应地由PTFE热化反应所产生的多种组分基于平动、转动、振动及电子模式所混合的比内能来表征[13]。

进一步，n0由PTFE的平衡蒸气压公式计算：

(7)

式中pc=1.847×1015Pa；Tc=20 815 K。

此外，外边界速度V1对式(5)的解算尤为关键。传统的烧蚀边界层模型为计算方便，通常假定V1为当地声速或接近声速。而本文则采用边界1和2处质量和动量守恒方程来确立V1，见式(8)。

(8)

式中n2、T2和m2分别为边界2处的密度、温度和粒子的平均分子质量。

相比传统的计算方法，式(8)考虑了等离子体层热化反应所造成的边界处组分平均分子质量的差异，即m1/m2的变化对V1的影响[14]。

综上所述， 式(5)～式(8)将烧蚀边界层与等离子体层的热力学状态参数紧密结合起来，形成可表征PTFE热化反应的耦合模型。

1.3数值计算方法

为了揭示PPT微秒量级工作过程中放电电离特性，需将上节定常模型与PPT放电过程与等离子体流动动态结合，构成一系统的数值计算体系，见图2。

图2　数值计算模型及流程图

热化反应模型为其余各模型提供物性参数。其中基于基尔霍夫定律的放电模型如下：

(9)

本文采用ode45算法来解算放电模型。放电电流所表征的焦耳热通过传热激励PTFE壁面，造成边界0处热力学状态参数(如T0)的动态变化，继而由烧蚀边界层和等离子体层耦合模型获取等离子体流动入口，即边界1处相应的热力学状态参数。

进一步，采用磁流体动力学方程(MHD)来描述层边界1和边界2间等离子体流动：

(10)

(11)

(12)

(13)

式中ρ、u和p分别为等离子体的密度、速度和压强；B为磁感应强度；j为电流密度；x为从推进剂烧蚀表面指向放电室出口的方向。

为了提高计算效率同时保证稳定性，本文采用显式二步MacCormack格式对上述等离子体层流动模型进行求解[15]。

以电极板间刚开始放电为初始条件，即在0～6 mm长度的放电室中：ρ=ρ0，u=0，T=T0，B=B0。其中，初始密度ρ0依据通道内初始参数设定；初始温度T0=300 K；初始磁感强度B0=0。

对于边界条件，采用图1中烧蚀边界层外边界1处的热力学状态参数作为流动入口边界参数， 即

(14)

考虑到PPT放电工作过程中，其出口流动可达到超声速条件，即出流边界条件(边界2)可采用简单的一阶外推获得。

2　计算结果与讨论

本文以常用的LES-6 PPT电参数(C=2 μF，L0=34 nH，R0=30 mΩ)和几何构型(h=30 mm，w=10 mm)等参数[15]为仿真对象，据图2所示的数值模型及计算流程，开展PPT放电工作过程的数值模拟研究，揭示PPT放电室内多种等离子体组分的时空演化机制。

图3是LES-6 PPT 脉冲放电工作过程中C元素等离子体数密度(分别是一价至三价的C离子，即C+、C2+和C3+)在放电室空间内随时间变化的等值线分布图。图3的纵向坐标表示放电时间，横向坐标表示放电室的长度(从推进剂烧蚀表面到放电室出口)。通过对比图3(a)、(b)和(c)可看出，C2+和C3+数密度的量级是1021，比C+所对应的1023小了2个数量级，这表明LES-6 PPT工作过程中烧蚀并电离出的C离子主要是以一阶电离形式存在。进一步由C+、C2+和C3+的时空分布可揭示C离子的产生和演化规律：在脉冲放电初期(0～1 μs)，放电电流最大，在放电室上游区域烧蚀并电离出大量C+离子，同时相对应于大脉冲电流所产生的强自感磁场-电磁加速效应非常快，刚刚电离的一价离子迅速被加速至放电室下游通道；C离子在下游加速同时进一步电离为高价离子，尤其是在0～0.5 μs时间范围内，出现了较高浓度的C2+和C3+；在脉冲放电中后期(1～3 μs)，随着放电半周期结束，放电能量开始减弱，电离出三价离子的概率降低，使得下游通道会呈现C2+，即C2+在下游分布存在的时间略晚于C3+；同样这也类似于C+。随着放电能量进一步降低，C2+的密度大幅下降，此时放电室下游通道及出口主要是C+的分布。

为分析烧蚀产生的中性气体的分布特性，图4给出了相应C原子数密度在放电室空间内随时间变化的等值线分布图。

(a)C+　　　　　　(b)C2+　　　　(c)C3+

图4　C原子数密度的时空分布图

由图4可知，C原子最大密度存在于脉冲放电初期放电室上游区域，并随着下游区域逐渐下降，与C离子的变化趋势一致。很显然这是由于C离子基本是C原子的分布进一步电离产生的。从下游区域来看，C原子主要集中于脉冲放电中后期，此时放电能量减弱，烧蚀组分的电离程度较低。C原子与其离子分布的最大区别是在脉冲放电末期及放电结束后(>3 μs),此时由于PPT所特有的延时汽化效应，放电室内依然充满一定浓度的中性C原子。

图5相应给出了 F元素等离子体(分别是一价至三价的F离子，即F+、F2+和F3+)数密度在放电室空间内随时间变化的等值线分布图。比较图5各价F离子峰值数密度可知， LES-6 PPT 烧蚀并电离出的F离子主要是以一阶电离形式存在。对比图3，F离子的峰值数密度比C离子小了一个数量级，这表明PPT放电电离过程中所产生等离子体的主要组分是以C+为主。此外，F离子在空间分布上与C离子产生较大的差异。不同于C离子主要集中在放电室上游区域，脉冲放电初期F离子主要在下游区域内被电离产生，这进一步说明，F元素的电离程度要远低于C元素。造成的原因主要是F原子的外层电子所激发跃迁的能量相对大于C原子，尤其体现是在更价的离子上。对比F+、F2+和F3+的时空分布可看出，与C离子类似，越高价的F离子在下游区域内电离产生的时刻越早，尤其是在0～0.5 ms内，大量的F3+产生并加速排出放电室，继而随着放电能量的降低，下游区域内分别呈现F2+和F+的分布。

(a)F+　　　　　(b)F2+　　　　　(c)F3+

图6描述了F原子在放电过程及放电结束后在推力室空间内随时间变化的等值线分布图。

图6　F原子数密度的时空分布图

对比图4和图6可知，中性F原子的峰值密度要比C原子高一个数量级，这主要是由于PTFE推进剂分子式C2F4所决定的。由图6可知，在脉冲放电初期，F原子集中于放电室上游区域并随着下游放电通道逐渐下降，与C原子的变化趋势一致。在整个脉冲放电中后期，F原子浓度持续从推进剂向下游区域内扩散，此时放电能量减弱，相应组分的电离程度较低。同样，由于PPT所特有的延时汽化效应，脉冲放电末期及放电结束后(>3 μs)放电室内存在着以F原子为主(密度大于C原子)的中性烧蚀气体。

数值模拟的结果需要实验数据的支持。Thomassen等人在LES-6 PPT 羽流排出产物诊断中发现存在C、F原子及其一价、二价、三价电离形式的等离子体[16]。为了与实验结果进行对比，本文分别给出了放电室出口处C、F原子及其离子数密度变化的数值解，见图7。图7表明在放电初期PPT烧蚀并电离产生了二价和三价的C、F离子；电离越充分的组分持续的时间越短。

进一步，表1对比了C、F原子及其各价离子数密度峰值时刻的数值解与Thomassen等测试的实验值。从表1可看出，数值结果基本符合实验值的范围。

(a)C原子

(b)F原子

粒子CC+C2+C3+FF+F2+F3+实验值1.51.20.80.31.30.80.50.3数值解1.91.280.80.581.91.00.760.6

PPT电离特性也可由电子数密度的时空分布来体现，见图8。从图8中可看出脉冲放电初期电离程度较高，大量C、F原子的外层电子激发跃迁最终形成大量自由电子，同时也相应产生高价的C、F离子；在脉冲放电中后期，放电室空间电子数逐渐下降，这主要是由于一方面随着放电能量的减少降低工质气体的电离度及相应产生的电子数量，另一方面可能是由于复合作用所导致。当C、F离子与自由电子碰撞时可能会俘获电子从而降低所带电荷数或形成中性原子，尤其是C、F离子在电磁作用下朝放电室下游区域加速过程与电子复合效应最明显(该区域电子数浓度最低)，从而导致高价C、F离子的损失。

图8　电子数密度的时空分布图

综上所述，PPT性能贡献来源可简单分为等离子体电磁加速及中性气体的气动效应，且有必要通过下式对这两者定量分析，揭示其效率影响机制。

(15)

(16)

式中mi为带电粒子质量；IEM为电磁加速冲量；MC、niC、MF、niF分别为放电室出口处C或F离子的质量和密度；u为速度；A=h×w。

等离子体质量变化规律及电磁加速对推力器性能贡献的计算结果见图9。可看出，等离子体的质量(即C、F各价离子质量)约为1.2 mg，而相应电磁冲量贡献约为19.4 μN·s。相对于文献[15]所给出的LES-6 PPT仿真结果(采用传统Mikellides模型，等离子体质量约2.6 μg)，本文采用改进的基于热化反应的烧蚀边界层与等离子体耦合模型所获得的等离子体的质量与Thomassen等测试结果更为吻合(等离子体的质量约占总脉冲烧蚀质量10 μg的10%，即1 μg)。

图10为LES-6 PPT推进剂壁面温度的数值解。TG-DSC测试结果[17]表明PTFE的分解温度(质量烧蚀损失起始温度>740 K)远大于其熔融点600 K。对比传统模型，改进模型所计算的壁面温度在主放电结束后依然保持在700～800 K，这是能够造成延时汽化现象的必要条件，更加符合PPT实际脉冲放电末期工况。进一步对比LES-6 PPT元冲量可知：在一次脉冲放电过程中，仅有10%的烧蚀质量被电磁加速，这些基于电磁加速的带电粒子对元冲量的贡献约占总量的70%；烧蚀质量绝大部分是脉冲中后期以及主放电结束后由延时气化效应产生的中性气体。这些基于气动效应加速的中性烧蚀质量，对元冲量的贡献仅占总量的30%，是造成PPT低效率的一个重要因素。

图9　LES-6 PPT 等离子体质量和电磁加速冲量

图10　不同模型下LES-6 PPT 壁面温度的数值解

近年来，基于放电能量多点多时段分布式匹配释放方式来提高PPT效率的途径受到广泛关注。这种方法是在PPT主放电结束后，期望通过第二次(或多次)脉冲放电来电离和加速延时汽化产生的中性气体来提高PPT性能。因此该方法能够有效应用的关键在于揭示放电室的电离情况，尤其是未电离的中性气体组分产生的时间及在电极间空间分布情况；而本文的研究可为摸清二次放电能量释放的最佳匹配时间和位置，进一步探索高性能PPT的放电模式提供定量化和系统化的理论研究手段。

3　结论

(1)基于热化反应所建立的烧蚀边界层和等离子体层耦合模型能够定量化地描述PPT放电烧蚀与离子化过程所产生的多种组分在放电室内时空分布特性，其性能预测结果符合PPT实际工况。

(2)LES-6 PPT 的等离子体主要是以一阶电离形式存在。C、F元素的电离区域分别集中在上、下游放电通道内，并且C元素的电离程度远高于F元素。

(3)脉冲放电初期通道出口处呈现大量短期存在高度电离的二价和三价的C、F离子，并且电离越充分的组分持续的时间越短。

(4)增强脉冲放电后期烧蚀气体的电离程度，减少延时汽化效应，对提高效率是十分必要的。

[1]毛根旺，唐金兰,等. 航天器推进系统及其应用[M].西安：西北工业大学出版社，2009.

[2]Burton R L, Turchi P J. Pulsed plasma thruster [J]. Journal of Propulsion and Power, 1998, 14 (5):716-735.

[3]Hiroyuki K, Ryosuke N, Kimiya K, et al. Plasma acceleration processes in an ablative pulsed plasma thruster[J]. Physics of Plasma, 2007 (14):033506.

[4]Mikellides P G. Theoretical modeling and optimization of ablation-fed pulsed plasma thruster [D]. USA: The Ohio State University (Doctor), 1999.

[5]Cassibry J T, Francis Thio Y C, et al. Numerical modeling of a pulsed electromagnetic plasma thruster experiment [J]. Journal of Propulsion and Power, 2006, 22(3): 628-635.

[6]Henrikson E M, Mikellides P G. Modeling of ablation-fed pulsed plasma thruster operation using a new approach to the ablation process [R]. AIAA 2008-4645.

[7]杨磊，刘向阳，武志文，等.基于克努森理论的脉冲等离子体推力器烧蚀边界层数值研究[J]. 固体火箭技术, 2013,36(3): 346-351.

[8]Keidar M, Boyd I D. Ionization and ablation phenomena in an ablative plasma accelerator [J]. Journal of Applied Physics, 2004, 96(10):5420-5428.

[9]尹乐, 周进, 缪万波，等. 脉冲等离子体推力器推进剂烧蚀传热计算[J]. 推进技术, 2010, 31(5) : 623-628.

[10]Zaghloul M R. On the vaporization of PTFE and heated compound-materials in ablation-controlled arcs [J].Journal of Applied Physics: 2004, 95 (7): 3339-3343.

[11]Mitchner M, Charles H, Kruger J. Partially ionized gases [M]. John Wiley, 1994.

[12]赵镇南.传热学[M].北京：高等教育出版社，2002.

[13]Schmahl S C. Thermochemical and transport processes in pulsed plasma microthrusters: a two-temperature analysis [D]. USA: The Ohio State University, 2002.

[14]Yang Lei, Liu Xiang-yang, Wu Zhi-wen, et al. Study on the wall ablation of heated compound-materials into discharge plasmas based on a modified model [J]. Applied Physics Letters, 2014,104:084102.

[15]肖利杰, 刘向阳, 杨磊，等. 脉冲等离子体推力器烧蚀建模与仿真 [J]. 推进技术，2011,32(6): 788-793.

[16]Thomassen K I, Vondra R J. Exhaust velocity studies of a solid PTFE pulsed plasma thruster [J]. Journal of Spacecraft, 1972, 9(1): 61-64.

[17]Yang Lei, Liu Xiang-yang, Wang Si-yu, et al. Different ablation models for the wall-plasma interaction process in pulsed plasma thruster [J].High Voltage Engineering, 2013, 39(9): 30798-30805.

(编辑：吕耀辉)

Numerical analysis on ionization characteristics in pulsed plasma thruster discharge

YANG Lei1,3, LIU Xiang-yang2, TANG Hai-bin3, HUANG Yu-ping1, WANG Ning-fei2

(1.Beijing Research Institute of Precise Mechatronic Controls, Beijing100076, China； 2.School of Aerospace Engineering, Beijing Institute of Technology, Beijing100081, China ；3.School of Astronautics, Beihang University, Beijing100191, China)

Ionization process of propellant ablation products plays a vital role in pulsed plasma thruster (PPT). Its effects on plasma generation and behaviors directly determine the PPT performances. To accurately describe the ionization characteristics in PPT discharge, a dynamic coupling model of ablating boundary layer and plasma layer was presented on the basis of thermochemical reaction. Numerical solutions of LES-6 PPT using this model were studied, and the ablation process of PTFE excited by current pulses and the temporal and spatial variation of the different plasma species in the discharge channel were revealed. It is found that plasma is mainly composed of monovalent ions, of which carbon ions and fluorine ions are focused on the upstream and downstream discharge chamber, respectively. A large number of short-lived and highly ionized (divalent and trivalent ions) carbon and fluorine species are created in initial stage, and results of peak time of the various plasma species show a good agreement with plume experimental data. It is valuable for theoretically studying the discharge energy distribution and loss mechanism in PPT.

electric propulsion；pulsed plasma thruster；ablation；ionization

2015-03-27；

2015-04-15。

中国博士后科学基金(2015M570918)。

杨磊(1982—)，男，博士后，研究方向为空间电推进技术。E-mail:ylforeverbit@126.com

V439

A

1006-2793(2016)02-0194-07

10.7673/j.issn.1006-2793.2016.02.008