背景压强对电推进羽流场影响的数值模拟

2022-11-09翁惠焱蔡国飙郑鸿儒刘立辉张百一贺碧蛟

北京航空航天大学学报 2022年10期

翁惠焱蔡国飙郑鸿儒刘立辉张百一贺碧蛟

(1. 北京航空航天大学宇航学院, 北京 100083;2. 北京航空航天大学航天器设计优化与动态模拟技术教育部重点实验室, 北京 100083;3. 长光卫星技术股份有限公司, 长春 130000)

在电推力器地面长时间寿命试验过程中,真空舱内溅射产物返流沉积和舱内背景压强对推力器工作性能评估及羽流场参数诊断影响较大。由于背景压强存在,推力器放电室压强比在轨工作时高,推力器工作状态发生改变;舱内中性气体密度较真空环境高十几个量级,使电荷交换(charge exchange, CEX)碰撞发生频率提高、发生位置变化,探针测量的流场参数出现较大偏差;此外,背景压强还会影响电推力器稳定性[1]、栅极表面溅射[2]等参数的评估和测试,导致对推力器在轨工作状态的分析和设计不准确。目前,对于背景压强的影响研究主要有试验、粒子仿真和一维模型3 种方法。

试验研究集中在背景压强对推力器地面工作参数的影响[3],尤其是对于霍尔推力器,由于自身结构设计导致其受背景压力影响更大,试验研究也较多[4]。 Walker 和Gallimore[5]先后完成了推力器冷流与热试情况下的实验测量工作,并对冷流下的舱内压强分布进行了直接模拟蒙特卡罗(direct simulation Monte Carlo, DSMC)仿真分析,一维模型将真空舱内的计算区域虚拟化为若干区域,分别用公式代表真空泵吸附作用、返流和粒子入口流量,结合壁面条件,求解出整个流场中的压强分布。 Cai[6]给出了不同真空泵组合形式下的计算模型,并对地面试验真空舱中的中性气体压强分布进行了计算。 Frieman[7]进一步增加了底部圆顶泵并进行了分析。这些研究中一维模型和冷流仿真都造成了羽流流场的失真,无法对电推力器工作时的真实羽流场进行影响分析。

目前,在离子发动机羽流的数值模拟中维持背景压强的方法有以下2 种:①采用虚拟粒子维持背景压强。在每次迭代计算时生成临时的虚拟粒子,维持背景压强[8]。临时产生的虚拟粒子参与碰撞,改变计算粒子的相关参数。而虚拟粒子的参数变化不会被记录下来,在下一次需要虚拟粒子时重新生成虚拟粒子,其参数为最初设定的参数,改变状态的只有计算粒子[9]。 ②使用计算粒子维持背景压强。背景压强粒子按照独立的粒子类型参与计算[10]。认为背景压力是具有一定温度和压强的中性推进剂原子。方法1 由于不用记录背景气体粒子的相关参数和运动轨迹,能够大幅降低计算资源的占用,提高计算速度。方法2 对计算资源要求高,计算效率低。以往研究中采用计算粒子维持背景压强的研究较少,认为虚拟粒子能够满足精度要求。随着计算机技术的不断发展,目前的计算资源能够满足方法2 的实现,因此有必要对计算粒子模拟压强的情况进行深入分析。

李小平等[11]采用二维粒子网格/蒙特卡罗碰撞(particle in cell/Monte Carlo collision,PIC/MCC)方法对离子栅极系统的影响进行了仿真分析,其中推力器喷出的中性粒子采用计算粒子,背景中性原子密度根据克拉伯龙方程均匀布置。 Jian等[12]采用基于浸入式有限元方法的PIC(IFEPIC)方法和虚拟粒子的方式仿真分析了背景压强对羽流场及羽流对栅极溅射率的影响。 Korkut等[13]采用推力器喷出的粒子建立压强的方式进行了PIC-DSMC 方法的粒子仿真分析。王军伟等[14]采用计算粒子计算分析了舱内冷板布置方式对羽流场的影响。综合来看,以往羽流场仿真研究主要采用虚拟粒子建立压强,导致背景中性粒子分布失真。而对背景压强的仿真研究中为了加速粒子平衡,往往忽视了羽流场参数的计算精度,两方面研究中尚未开展虚拟粒子和计算粒子2 种建立背景压强方式下流场参数的分析对比。本文针对离子推力器羽流分别开展计算粒子、虚拟粒子和绝对真空情况下的仿真对比分析,得到背景压强及其不同计算方法对羽流场计算结果的影响。

1 仿真方法

本文涉及的仿真分析均使用电推进羽流效应仿真软件EX-PWS[15-16]实现。 EX-PWS 采用混合粒子网格(hybrid PIC)方法计算带电粒子在电场中的运动过程,使用DSMC 方法[17]处理粒子间的弹性碰撞和电荷交换碰撞。

1.1 控制方程

电子动量方程是EX-PWS 仿真软件计算电推进羽流效应时的控制方程,如下:

式中:me为电子质量;ne为电子密度;ve为电子速度;e为电子电荷量;E为电场强度;B为磁场强度;p为压强;vei为电子和离子之间的碰撞频率;vi为离子速度。

在计算本文涉及的离子推力器羽流场时,考虑到等离子体产生的自洽磁场及推力器自身磁场对于羽流流场的影响较小,因此在计算中忽略磁场的作用。此外,在本文主要考虑的羽流范畴,由于空间尺度较大,等离子体的德拜长度远小于分子自由程,等离子体总体上可视为准电中性,即可以认为电子密度等于离子密度。同时,混合PIC方法将电子视为流体,采用电子运动方程计算,而将中性粒子和离子视为计算粒子,参与流场中的迭代计算。

在描述等离子体中粒子间的碰撞时,通常用碰撞频率与等离子体频率的比值vei/ω≈lnND/ND来判断,其中ND为德拜球中的带电粒子数。当对于电子温度取值为1. 0 eV 时,其比值约为1.7 ×10-3,远远小于1,说明离子推力器的等离子体羽流流场可以视为无碰撞。此时,对于无碰撞同时忽略磁场影响的电子的动量方程可以演化为

结合理想气体状态方程(式(3)),则流场内等离子体的电势和电子密度之间的关系可以由玻尔兹曼关系式(式(4))给出[9]:

式中:ϕ为电势;nref为0 处的参考电子密度,本文取为1 ×1016m-3;Te为电子温度;k为玻尔兹曼常数。

从式(2) ～式(4)可以看出,为了确定电势,需要得到电子温度。对于氙气工质离子推力器的羽流流场中电子温度的实验测量结果,根据测量位置不同,从0.5 eV 到3.0 eV 不等[18]。

1.2 粒子运动

在获得各网格内粒子数密度后,需要根据粒子所处的位置给予相应的权重,将电荷分配至网格节点上。 EX-PWS 采用非结构化四面体网格,对于网格节点上的电荷储存位置采用面元中心,在分配电荷时采用体积分配法,如图1 所示。取粒子所在位置点与网格的4 个面组成的小四面体的体积和总四面体网格的体积之比作为各表面中心点处获得的粒子电荷分配权重。

图1 网格中电荷分配示意图Fig.1 Schematic diagram of charge distribution in the grid

EX-PWS 软件中推动粒子运动采用蛙跳算法,如图2 所示。蛙跳算法计算原理简单,且具有二阶精度,适合编程实现。采用蛙跳算法后,对于单个粒子的运动微分方程可以离散为如下有限差分方程:

图2 蛙跳算法示意图Fig.2 Schematic diagram of leapfrog algorithm

式中:m为质量;v为速度;F为粒子受到的力;x为粒子运动所在位置。

1.3 粒子碰撞

粒子间的碰撞是影响流场分布的一个主要因素,在电推进羽流流场仿真中,通常包括弹性碰撞及电荷交换碰撞。在EX-PWS 软件中,粒子间碰撞使用DSMC 方法进行处理。采用非时间计数器,分子作用势模型选用可变硬球模型。本文涉及的碰撞种类如式(7) ～式(12)所示。其中,MEX代表动量交换碰撞,CEX 代表电荷交换碰撞。

氙原子与一价氙离子和二价氙离子间弹性碰撞截面由Dalgarno 和Williams[19]给出:

式中:针对氙原子与一价氙离子的电荷交换碰撞k1= -0.882 1,k2=15.126 2。对于氙原子与二价氙离子电荷交换碰撞时的截面有k1= -2.703 8,k2=35.006。

2 计算条件

本文使用的推力器模型为兰州空间技术物理研究所研制的LIPS-200 型离子推力器,性能指标为:推力(40 ±4) mN,比冲(3 000 ±300) s,推进剂利用率大于80%,功耗1 300 W[21]。表1 给出了与仿真相关的推力器的工作参数,表2 给出了推力器的仿真参数设置[22]。文献[15]给出了采用本文计算模型得到的羽流场仿真结果与实验测量数据的对比分析,证明了本文计算方法的准确性。

表1 LIPS-200 型离子推力器工作参数[21]Table 1 Parameters of LIPS-200 ion thruster[21]

表2 LIPS-200 型离子推力器仿真基本参数[22]Table 2 Basic simulation parameters of LIPS-200 ion thruster[22]

图3 给出了背景压强计算时采用的计算域示意图。整个计算域为长度1.5 m、直径2 m 的圆柱体,内部设置推力器模型,为了表征返流区粒子分布情况,计算域包括了推力器前方1 m 和后方0.5 m 范围。推力器前方壁面设置外径1.0 m、内径0.8 m 的环形边界作为模拟真空泵表面,认为当粒子撞击该表面时则被真空泵吸附,在计算域中删除该粒子,即自由边界条件。本文计算中共涉及的3 个仿真算例初始条件设置如表3 所示。

表3 背景压强处理方式对比的算例安排Table 3 Case arrangement for comparison of background pressure simulation methods

图3 背景压强算例计算域示意图Fig.3 Simulation domain for background pressure cases

在计算粒子算例中,计算域中的初始压强为0,所有边界面设置为漫反射,仅在推力器前方环形泵的位置设置为自由边界。当推力器开始工作后,喷射出的粒子撞击计算域边界后被中和为慢速运动的原子,成为背景压强,除少数粒子被真空泵吸收外,大部分粒子充满整个计算域空间。参与计算的所有粒子均进行位置跟踪和信息更新。

在虚拟粒子算例中,将0.25 mPa 的背景气体视为虚拟粒子,仅参与碰撞计算。当计算粒子运动到某一网格中将要与背景原子发生碰撞时,根据克拉伯龙方程给出粒子数密度,按照环境温度根据Maxwell 分布给出热运动速度。碰撞后背景粒子状态不变,计算粒子运动速度及方向按照碰撞模型发生变化。计算粒子与背景粒子间的碰撞类型包括弹性碰撞及电荷交换碰撞。由于使用虚拟粒子建立压强,所有边界面设置为自由边界,认为背景压强恒定不变。

在绝对真空算例中,将背景压强设置为0,同时所有边界面设置为自由边界,认为气体分子撞击在表面后随即被吸收,从计算域中删除。绝对真空算例主要用来对比分析虚拟粒子和计算粒子2 种背景压强处理下的流场结果与绝对真空情况下的差别。

3 计算结果

本节分别针对3 种仿真算例中的中性粒子和电荷交换离子分布进行对比分析。

3.1 中性粒子分布

图4 给出了计算粒子建立背景压强的算例中中性粒子数密度分布。可以看出,中性粒子的空间分布可以分为5 个区域。区域①为返流区,以推力器出口为界,在推力器后方区域,这部分区域中的粒子主要是后方和侧方舱壁的反射,以及束流区中部分返流粒子。区域②为电荷交换离子所在的“翼”形区域,这一区域的粒子能量、数密度及运动方向等信息是研究人员普遍关心的。区域③为地面真空舱内试验时涉及的真空泵影响区,由于真空舱的真空泵布置不同,真空泵影响区的位置也不同,同时根据真空泵的抽气能力不同,仿真中设置的吸附系数也应该不同,需要进行折算。区域④为推力器前方舱壁返流影响区,这部分反射粒子由于运动方向与束流正好相反,朝向推力器运动,并与束流直接混合,会对羽流束流区的测量结果造成直接影响。区域⑤为束流区域,在推力器前方0.4 m 内,同时又根据束流发散角的限制,使得区域⑤成为一个类似圆台的结构。

图4 计算粒子建立的中性气体分布Fig.4 Neutral gas distribution created by computed particles

从图4 可以看出,对于采用计算粒子建立背景压强的算例,返流区中性粒子数密度约为7 ×1016m-3,分布较为均匀。电荷交换离子分布区域和返流区的情况相近。在真空泵影响区可以看出,真空泵周围粒子数密度迅速降低,但由于真空泵吸附面积有限,压强下降区域范围有限。而在束流前方舱壁的影响区可以看出,返流的粒子数密度基本与束流区持平,影响范围更广。

图5 显示了计算粒子算例中背景中性粒子的分布情况。计算时设定所有推力器喷出的离子和原子撞击壁面后成为背景中性粒子,此处不包含推力器喷出的中性粒子。从图中可以看出,在计算粒子算例中,影响中性粒子分布的主要是背景粒子。对比图4 可知,在达到平衡后,背景中性粒子数密度远高于推进剂中的中性粒子。

图5 计算粒子建立的背景压强中性气体分布Fig.5 Neutral gas distribution of background pressure created by computed particles

图6 给出了使用虚拟粒子建立背景压强时的中性粒子空间分布情况。为了方便与计算粒子的算例进行对比,在图6 中已经叠加了计算时采用的虚拟背景压强0.25 mPa。根据文献[23]中给出的压强与吸附面积的关系式计算得到,当真空泵设定为吸附系数为1. 0 的外径1. 0 m、内径0.8 m的圆环时,处于热平衡状态下的压强约为0.32 mPa,与本文选取的0.25 mPa 较为接近。存在一定差距的原因是:热平衡状态下,认为粒子宏观运动速度为0,而推力器喷出的粒子速度远远大于0。虚拟粒子中背景压强0.25 mPa 取值于计算粒子算例中推力器出口径向位置的平均压强。

图6 虚拟粒子建立的背景压强中性气体分布Fig.6 Neutral gas distribution of background pressure created by virtual particles

可以看出,与图4 相比,除束流区域⑤的云图比较相近外,其他区域都有所不同。这说明虚拟粒子算例无法表达粒子与边界面的相互作用,对返流区等粒子数密度较小的区域参数分布影响较大。图7 给出了绝对真空条件下的羽流场。可以看出,中性粒子数密度明显降低,分布范围也仅在推力器出口附近。

图7 绝对真空环境下的中性气体分布Fig.7 Neutral gas distribution in absolute vacuum case

图8 显示了中心轴线上中性粒子数密度的对比。图中黑色虚线代表计算粒子算例中的背景中性粒子[BG];红色实线代表计算粒子算例中流场的中性粒子总密度,包括未电离推进剂工质和背景气体;绿色虚线代表虚拟粒子算例中的中性粒子总密度;蓝色点划线代表真空算例中的流场中性粒子总密度,无背景气体。

图8 推力器出口轴线上不同模拟情况下中性粒子数密度对比Fig.8 Comparison of neutral particle number density at thruster outlet axis in different simulation cases

显然,在不受真空舱壁面影响的真空情况下,中性粒子数密度远低于在真空舱内存在背景压强的情况。计算粒子算例中给出的羽流中的中性粒子数密度与虚拟粒子算例中给出的结果在推力器附近区域相近,在算例中距离出口0.4 m 之内,表明在该区域推进剂粒子占主导地位。在距离推力器出口0.4 m 以外区域,从真空舱壁返回的中性背景粒子占主导地位,表明真空舱壁对羽流场具有显著影响,真空舱尺寸越小,影响越大。同时,从分析中也可以看到,当真空舱尺寸足够大时,真空舱壁对羽流场的影响有限,在推力器近场区域以推进剂粒子为主,在该区域计算粒子与虚拟粒子仿真结果差异并不明显,而虚拟粒子可以大幅提高计算速度,因此,虚拟粒子计算方法更适合近场羽流效应计算。

图9 给出了距离推力器出口0.5 m 半径范围内中性粒子数密度随角度变化,0°代表轴线x=0.5 m 位置,90°代表与出口平面共面,垂直于轴线0.5 m 位置。可以看出,背景中性粒子数密度在计算粒子算例中占主导,与虚拟粒子数密度相近,两者沿径向分布较为均匀。而在真空算例中,中性粒子数密度则随角度增大而降低,整体较有背景压强的算例低1 ～2 个数量级。这反映出2 个现象:①计算粒子算例在达到平衡后,流场中心区域分布较为均匀。 ②在地面试验过程中,由于舱内背景压强的限制,必然导致流场失真,无论采用计算粒子还是虚拟粒子进行仿真都无法避免。在本文算例中,虚拟粒子压强按照计算粒子算例中的推力器径向位置压强设置,在实际工程中往往需要根据真空规获得舱内压强,真空规的布置位置对测量结果较大。

图9 推力器出口0.5 m 径向中性粒子数密度对比Fig.9 Comparison of radial neutral particle number density at thruster outlet of 0.5 m

3.2 电荷交换离子分布

电荷交换离子由于易受到电场影响轰击航天器表面,其分布情况是工程部门关心的重点之一。图10 和图11 给出了3 种工况下轴线和径向0.5 m半径上的电荷交换离子数密度分布对比。可以看出,在计算粒子和虚拟粒子情况下,电荷交换离子数密度分布相似,但计算粒子的算例中数密度要比虚拟粒子的情况稍高。在真空算例中,在推力器出口处,电荷交换离子数密度与有背景压强的情况相近,但随着距离出口的距离越大,数密度逐渐降低,在距离推力器出口1 m 处的数密度较出口处下降了2 个数量级以上。而有背景压强的算例则在轴线上变化不大。此外,可以看出轴线对比中数值波动较大,可能与轴线网格较为稀疏及电荷交换碰撞发生的随机性有关。

图10 推力器出口轴线上不同模拟情况下电荷交换离子数密度对比Fig.10 Comparison of CEX ion number density at thruster outlet axis in different simulation cases

图11 给出的径向分布中电荷交换离子数密度先下降,在40°左右位置又升高,这与电荷交换离子的分布有关,电荷交换离子由于速度较低,容易在电场驱动下沿径向运动,在推力器出口附近产生“翼”形区域。因此,在半径0.5 m 径向取值时会先下降后升高。

图11 推力器出口0.5 m 径向电荷交换离子数密度对比Fig.11 Comparison of radial CEX ion number density at thruster outlet of 0.5 m

仿真中,离子的运动由网格中各节点间的电势差驱动。为了进一步分析电荷交换粒子分布情况,图12 和图13 给出了在计算粒子和真空环境条件下的电势分布云图。可以看出,电势均为负值,这是由设定推力器出口平面为电势零点造成的,推动粒子运动的是各节点间电势差。通过对比分析可以看出,计算粒子与绝对真空环境算例中电势在轴线上变化较为接近,但在径向上,计算粒子电场作用范围更广,但电势梯度相对较低。

图12 计算粒子算例中的电势分布Fig.12 Potential distribution in computed particle case

图13 绝对真空环境下的电势分布Fig.13 Potential distribution in absolute vacuum case

图14 和图15 分别给出了计算粒子和真空环境条件下的电荷交换离子数密度分布云图。根据前面分析可知,计算粒子和虚拟粒子工况在轴线和径向上的电荷交换离子数密度均较为接近,都明显高于绝对真空环境。通过对比可知,在计算粒子算例中,轴线束流区是电荷交换离子分布的主要区域,而在绝对真空情况下,电荷交换离子主要集中在羽流核心区域和推力器出口附近的“翼”形区域。说明舱内残留的中性气体的存在明显提高了电荷交换离子的产生,导致电荷交换离子分布与在轨绝对真空环境相差较大。因此,如何在地面真空舱内测量得到电荷交换离子的参数性质,对仿真结果进行验证分析是电推进羽流场测量的一个难点。