外磁场作用下的磁等离子体动力学过程仿真

罗杨 陈茂林 苏冬冬 许诺 王忠晶 韩志聪 赵豪

1) (北京机械设备研究所,北京 100854)

2) (西北工业大学,燃烧、热结构与内流场重点实验室,西安 710072)

磁等离子体动力学推力器是空间高功率电推进装置的典型代表,磁等离子体动力学过程是其核心工作机制.为深入理解外磁场对其工作特性的影响,本文采用粒子云(particle in cell,PIC)方法结合基于自相似准则的缩比模型,进行外加磁场作用下磁等离子体动力学推力器工作过程的建模仿真,通过与实验结果对比验证模型和方法的可靠性,并重点分析推力器点火启动过程的等离子特性参数分布,以及外磁场和阴极电流对推力器工作性能的影响.研究结果表明:阴阳极放电电弧构建是推力器启动和高效工作的关键步骤;外磁场强度较低工况不利于构建稳定放电电弧,等离子体束流集中于轴线附近,推力主要产生机制是自身场加速;外磁场强度较高时,阴阳极放电电弧稳定,推力产生主要机制是涡旋加速,推力、比冲随外磁场强度线性增大;推力器效率随阴极电流和外磁场强度增大而增大;放电电压随阴极电流增大而增大,但随外磁场强度的增大表现出先减小后增大的趋势.

1 引言

磁等离子体动力学推力器(magnetoplasmadynamic thruster,MPDT)是通过洛伦兹力产生推力的航天推进装置,也称洛伦兹力加速器,可实现兆瓦级功率和百牛级推力性能,是高功率、大推力电推力器的重要发展方向.MPDT 在美国、前苏联/俄罗斯、日本、德国、意大利、中国均得到了广泛关注和重点研究[1,2],实验室内实现了数千瓦至数兆瓦功率的MPDT 测试,最大测试推力达100 N 以上.

研究人员对MPDT 进行了大量的理论分析、仿真模拟和实验测试等工作.在仿真模拟研究方面早期主要采用流体模型:Brushlinskii 和Morozov[3,4]最早于20 世纪60 年代采用流体模型对大功率MPDT 开展数值模拟,获得了MPDT 轴向密度分布与电流的关系;Kimura 等[5]提出了单温度的二维计算模型,并使用该模型对MPDT 的工作过程进行模拟;Niewood[6]和Sheppard 等[7]提出了双温度轴对称数值模型来研究计算推力器内部的流场分布,考虑了非平衡电离过程和中性气体的影响;Caldo 等[8]采用双温度模型来研究了反常运输对MPDT 流动的影响.Rudolph[9]在MPDT 流体模拟过程中采用多重网格和多重时间步长法,解决了时间依赖性问题.Heiermann[10]、赵博强等[11]、Sary等[12]、Tahsini[13]采用流体模型对MPDT 内部的电离过程和多种输运过程进行了详细建模仿真,实现了对MPDT 推力器的整体性能计算.

流体模型适用于MPDT 的高密度等离子体特性,具有计算适用性强、计算约束条件少、计算速度快的特点.但是流体作为基于宏观的平均化方法,流体模型无法捕捉到等离子体的一些瞬时现象.粒子云(particle in cell,PIC)[14]等基于动力学过程的等离子体粒子模拟方法是捕捉MPDT工作过程中等离子体演化过程和瞬时现象的优选方法.由于MPDT 高等离子体密度特性对PIC 仿真时间步长和空间步长的限制,使得其对计算资源和计算机时需求巨大,故早期粒子模拟方法在MPDT仿真中应用较少.Tang 等[15-17]开展了PIC方法在MPDT 模拟中的应用探索,采用改变重粒子质量的加速计算方法减少了计算机时,实现了低功率下MPDT 工作过程模拟;Li 等[18]采用PIC方法开展了无外加磁场作用下的MPDT 工作过程模拟.本文采用PIC 方法结合基于自相似准则的缩比模型进行外加磁场作用下的MPDT 工作过程模拟,解决MPDT 极高等离子体密度对模拟时间步长、空间步长的限制,在减小计算量的同时,进行真实工况下MPDT 工作过程模拟和性能影响因素分析.

2 模型与方法

2.1 MPDT 结构与工作原理

根据磁场产生方式的不同,MPDT 分为两大类:自身场磁等离子体动力学推力器(self-field MPDT,SF-MPDT)和附加场磁等离子体动力学推力器(applied-field MPDT,AF-MPDT).本文主要开展AF-MPDT 工作过程的仿真研究,其主要结构包括阴极、阳极和电磁线圈,如图1 所示.

图1 AF-MPDT 结构和工作原理示意图Fig.1.Schematic of AF-MPDT structure and operating principle.

AF-MPDT 主要有3 种加速机制产生推力:气动加速、自身场加速和附加场加速,附加场加速又包括霍尔加速和涡旋加速.

1)气动加速.气动加速分量是通过电弧欧姆加热将电能转换为热能来提高推进剂的总焓,进一步喷管将热能转换为动能而产生推力.当 MPDT在高功率工况工作时,电磁加速占主导,气动加速分量可忽略.

2)自身场加速.图2(a)为自身场加速原理.MPDT 工作时,电流轴向分量会产生一个周向自感磁场Bθ,该磁场又与电流相互作用产生洛伦兹力:其中径向电流jr与周向自感磁场Bθ,产生一个轴向的力 (jr×Bθ),该作用力直接对推力产生贡献;而轴向的电流jz与自感磁场Bθ相互作用产生径向的洛伦兹力 (jz ×Bθ),该作用力会导致轴线附近的等离子体密度增大,并产生一个不平衡的压力,间接对推力产生贡献.由于感应磁场Bθ与电流线性相关,自身场加速产生的推力与电流表现为平方关系.自身场加速机制在SF-MPDT 中贡献较大,而在AF-MPDT 中贡献较小.

图2 MPDT 加速机理 (a)自身场加速;(b)霍尔加速;(c)涡旋加速Fig.2.MPDT acceleration mechanism:(a) Self-field acceleration;(b) hall acceleration;(c) vortex acceleration.

3)霍尔加速.图2(b)为霍尔加速原理.在AFMPDT 中,外部电磁线圈会在阴阳极间产生轴向磁场Bz和径向磁场Br.根据欧姆定律,轴向磁场Bz和径向电场作用下,等离子体会感生出周向电流jθ,与自身场机制相似,jθ会与附加磁场产生箍缩力 (jθ×Bz)和轴向力 (jθ×Br).

4)涡旋加速.图2(c)为涡旋加速原理.在AFMPDT 中,阴阳极间的轴向和径向电流分量,也会与外部电磁线圈产生的径向磁场Br和轴向磁场Bz相互作用,产生周向电磁力 (jr×Bz) 和 (jz ×Br),使等离子体产生角向涡旋,宏观上表现为自旋的等离子体径向膨胀,在物理/磁喷管的作用下再转化为轴向作用力,产生推力.涡旋加速产生的推力与电流和外加磁场的乘积IB成正比,Fradkin[19]的实验也验证了该推力关系.

2.2 物理模型

为了便于仿真结果与实验结果对比分析,采用Myers 团队设计的AF-MPDT 推力器[20]作为计算模型,其几何结构如图3 所示,图中标注单位为mm,虚线框内为仿真模拟区域.

图3 MPDT 几何结构与计算区域Fig.3.MPDT geometry and simulation area.

2.3 仿真方法

2.3.1 PIC 方法

PIC 方法是等离子体粒子模拟的常用方法,其核心过程包括电磁场求解、带电粒子在电磁场作用下运动、以及粒子碰撞三部分[21].

MPDT 工作过程中的电磁场分布及演化可以通过麦克斯韦方程描述,由于MPDT 工作过程中电磁场的时变量较小,可以忽略,故可将麦克斯韦方程简化为

式中,B=∇×A为磁感应强度,E=-∇φ为电场强度,μ0为真空磁导率,J为电流,ρ为电荷密度,ε0为真空介电常数,A为矢量磁位,φ为电势.进一步将(1)式改写为泊松方程形式,便于数值求解:

带电粒子的运动通过牛顿-洛伦兹方程进行描述,如下:

式中,m为粒子质量,v为粒子速度,e为单位电荷,s为粒子位置,t为时间.

PIC 方法处理粒子间的碰撞一般使用蒙特卡罗碰撞(Monte Carlo collision,MCC)方法[22].在MPDT 模拟过程中,考虑如下碰撞类型:电子与中性分子间的弹性碰撞、激发碰撞、电离碰撞,离子与中性分子间的弹性碰撞和电荷交换(charge exchange,CEX)碰撞,电子-电子、电子-离子、离子-离子间的库仑碰撞,以及电子的Bohm 碰撞.

2.3.2 边界条件设置

仿真中,主要包含阴极、阳极、绝缘壁面、对称轴和自由边界5 种边界类型.阴极边界设置为Neumann 边界,固定法向场强E⊥=σ/ε0,切向场强为E//=0,其中σ为阴极表面电荷密度;固定阴极的轴向电流密度为,径向电流密度和周向电流密度均为0,其中,I为阴极电流,Rcathode为阴极半径;阴极作为粒子吸收边界,粒子碰撞阴极后注销.考虑离子碰撞阴极壁面的二次电子发射过程,对于钨阴极,当氩离子的入射能量介于10—100 eV 之间时二次电子发射系数为0.095 左右.

阳极设置为零电势的Drichlet 边界和粒子吸收边界,粒子碰撞阳极后注销.阴极阳极间的绝缘壁面为Neumann 边界,可认为壁面表面的场强是由累积电荷得到,且垂直于壁面,即法向场强E⊥=σ/ε0,切向场强为E//=0.对于离子,绝缘壁面为吸收边界,离子碰撞绝缘壁面后注销;对于电子,分别考虑在壁面沉积、反射、激发一个二次电子、以及激发两个二次电子的情况,4 种情况的概率分别为

式中,ξ为电子能量,Pdep(ξ)为电子沉积壁面的概率,Pref(ξ)为电子发射反射的概率,P1se(ξ)为电子撞击出一个二次电子的概率,P2se(ξ)为电子撞击出两个二次电子的概率.

自由边界主要是粒子的出口,粒子运动越过自由边界后注销;计算模型是二维轴对称模型,中心轴线为对称边界,粒子运动出对称边界时进行镜面反射处理;电磁场方程均使用Neumann 边界为自由边界,可认为电场、磁场在自由边界处的法相分量为零.

2.4 基于自相似准则的缩比模型

PIC 仿真过程中,空间步长需要小于德拜长度,时间步长需要小于等离子体频率,并同时满足CFL 条件.对于MPDT,电子数密度一般大于1020m—3,对应德拜长度一般小于10—6m,单个自由度的仿真网格数一般需要数万至数十万,等离子体频率约1012Hz,对应时间步长约10—12s.要同时满足空间步长和时间步长需求,也会给计算量带来极大挑战,因此需要一些加速方法在保证仿真准确性的同时,减少计算量.

目前对全粒子PIC 模拟有两大类加速方法,一种是直接改变介电常数以及重粒子质量来加速计算的方法[23],另一种是Khayms[24]和Taccogna等[25]在霍尔推力器仿真中使用基于自相似准则的缩比模型方法,其基本特征是保证推力器性能参数成比例缩放.

鉴于MPDT 和霍尔推力器工作过程和物理特征相似,本文采用基于自相似准则的缩比模型方法,结合PIC 粒子模拟方法,进行MPDT 仿真模型构建和加速计算.为了获得准确的推力器性能参数,缩比模型中,相似准则选取如下:1)推力器的比冲Isp不变;2)磁场对电子的约束效果不变;3)电子平均自由程与推力器特征尺寸比值不变.

缩比模型中,将推力器特征尺寸缩小至1/τ,则物理量电流I、质量流量、数密度ne,ni、电场E、磁场B、推力F、比冲Isp和效率η的缩比系数分别为1/τ,1/τ,1/τ,τ,τ,τ,1,1 和1/τ.其中自变量为计算输入时需要调整的量,因变量表示根据调整后的输入量而发生改变的物理量.缩比模型将推力器尺寸缩小为原来的1/τ后,电子数密度同样也扩大为原来的τ倍,使得德拜长度缩小1/τ1/2,可以认为缩比系数τ对单方向的网格数量的缩小量为1/τ1/2,对于轴对称模型,网格数量缩小1/τ.本文仿真中τ取为8 × 106.

3 仿真结果与分析

3.1 中性气体流场分布

中性气体流场分布一般采用直接蒙特卡罗模拟(direct simulation Monte Carlo,DSMC)方法求解,将等离子体放电和输运过程模拟的PIC 方法与中性气体流场模拟的DSMC 方法耦合,会大幅增大计算需求.本文将中性气体流场与等离子体流场解耦处理,即先单独采用DSMC 方法模拟MPDT 中的中性气体流场,然后采用电推进仿真中常用的PIC/MCC 方法处理带电粒子与中性气体原子间的相互作用,即将具备一定分布的气体原子作为背景,参与带电粒子的碰撞处理,但其分布不随模拟过程变化,这样可在保证较好的计算精度的同时,大幅缩减计算机时.

针对Ar 工质,流量为=0.1 g/s 的情况下的气体初始流场,Ar 原子数密度分布和流线如图4所示,流场内中性粒子数密度量级约在1020—1021m—3量级.

图4 Ar 原子数密度分布和流线Fig.4.The number density and streamlines of Ar.

3.2 等离子体密度分布演化

通过推力器内部等离子体特性参数的演化,可以分析推力器内部放电和带电粒子的输运过程.针对外磁场B=0.03 T,阴极电流I=1000 A 工作参数下的MPDT 的工作过程进行仿真,在缩比模型中,计算时间步长为1 × 10—18s,对应实际情况的时间步长8 × 10—12s.

3.2.1 电子数密度

MPDT 的放电由阴极发射的电子电流主导,故首先分析放电过程中的电子数密度分布随时间的演化过程,不同时间步的电子数密度分布仿真结果如图5 所示.

图5 电子数密度分布Fig.5.Distribution of electron number density.

MPDT 点火启动过程中,电子数密度分布随时间的变化可初步分为3 个阶段:放电初期,电子从阴极发射后,受外加磁场约束,沿磁力线螺旋运动,呈扩张型分布;放电中期,电子在阴极尖端下游堆积,逐渐形成一个高密度分布区域,并向阳极方向延伸,构成一条连通阴阳极的高密度电子分布带,形成放电电弧;放电末期,放电电弧逐渐稳定,等离子体区域进一步扩张,并在阴阳极下游形成完整的高密度的稳态分布,推力器进入稳定工作状态.推力器在达稳态工作状态后,在推力器外部,电子束流分为两个部分,一个集中在阴极尖端的轴线处向下游扩散,另一部分电子随磁力线向下游扩散,这和实验观测到的现象一致[1].稳态时,MPDT阳极尖端的电子数密度达1021m—3,阴阳极下游的电子数密度达1020m—3.

3.2.2 离子数密度

阴极发射电子在电磁场共同作用下运动的同时,会碰撞中性气体分子发生电离,产生新的电子和离子.由于离子完全由电离过程产生,离子数密度分布随时间的演化过程更真实地反映MPDT 的电离过程,不同时间步的离子数密度分布仿真结果如图6 所示.

由图6 可知,在放电初期,离子数密度同样由小增大,并随着放电进行,其分布也逐渐稳定.在放电初期,离子数密度的分布与电子分布相似,其原因是离子由电子电离碰撞产生的,在高电子数密度处离子数密度高;随着放电进行,离子数密度的分布逐渐偏离电子分布,呈现出高密度区域逐渐由阴极尖端附近向整个下游区域扩张并向轴线附近聚焦的趋势,离子分布偏离电子分布的原因是离子质量大,受磁场的约束远小于电子,故运动轨迹和最终分布也与电子不同.稳态工作时,离子数密度最高达到1021m—3,主要仿真区域离子数密度约1020m—3,与电子数密度相当.

图6 离子数密度分布Fig.6.Distribution of ion number density.

3.3 外加磁场对MPD 性能影响及其与实验数据的对比分析

外加磁场是影响MPDT 推力器性能的重要参数,在阴极电流I=1000 A,流量为=0.1 g/s工况下,开展外磁场B=0.008,0.017,0.025,0.034,0.052 和0.069 T 情况下的MPDT 工作过程仿真,分析不同磁场强度下的电子数密度分布和推力器性能参数,并进行仿真结果与实验结果对比,验证仿真模型和方法的有效性.

3.3.1 不同外磁场的电子数密度分布

不同外磁场作用下,MPDT 稳态工作时的电子数密度分布如图7 所示.

图7 不同磁场强度下的电子数密度分布Fig.7.Electron number density distribution with different applied magnetic field.

从图7 可知,在B=0.008 T 情况下,电子仅在阴极下游的轴线附近形成高密度区域,未能形成高效的阴阳极电弧放电;随着磁场强度的增大,中轴线高密度区域逐渐向径向移动,在外磁场强度B=0.017—0.340 T 之间时,电子数密度出现两个高密度区域,分别位于轴线附近和推力器出口下游区域,构成稳定的放电电弧;外磁场强度继续增大,轴线附近的高密度区域逐渐缩小,当磁场B＞0.052 T 时,中心轴线的高密度区域基本不存在.分析认为:在外磁场强度较小时,推力器的自感磁场比重较大,自感磁场主要为周向磁场,约束电子在轴线附近聚集;随着外磁场的增加,自感磁场的比重降低,外磁场对电子的约束成为主要部分,外磁场将使得电子沿着磁力线向下游运动.

3.3.2 外磁场对推力器性能的影响

根据不同外磁场作用下的仿真结果,计算了MPDT 稳态工作时的放电电压U、推力F和效率η,如图8 所示.

由图8 可知,当B＜ 0.053 T 时,放电电压随外磁场增大而降低,B≥ 0.052 T 时放电电压逐渐变大.其原因是阴极表面电荷密度决定了阴极电势即放电电压的大小,而阴极表面电子的运动受磁场影响严重:当外磁场较小时,自感磁场对电子的约束占据主导,电子向阴极聚集,增大了阴极表面电荷密度,放电电压升高;随着外加磁场增大,电子逐渐跳出自感磁场约束,远离阴极,使得放电电压降低;当外磁场进一步增大时,外加磁场对电子约束增强,电子难以跨越阴极磁力线沿径向扩散,使得阴极表面电荷密度增大,放电电压再次升高.

图8 MPDT 性能随外磁场的变化Fig.8.MPDT performance changes with applied magnetic field.

不同磁感应强度条件下的推力结果表现为随磁感应强度B近似线性增大,与Myers[20]的实验结论一致.在B=0.069 T 情况下的推力值为0.98 N,略小于Myers[20]的实验结果1.16 N,误差产生的原因可能是本文仿真忽略了气动推力的贡献;推力器的效率同样表现为随磁感应强度B增大而增大的趋势.

3.3.3 与实验测试数据对比分析

Myers[20]实验测量了不同外加磁场情况下的MPDT 推力器的比冲、效率等性能参数,为了进一步验证本文仿真程序的准确性,将仿真结果与Myers[20]的实验结果进行对比,如图9 所示.

图9 仿真与实验结果对比Fig.9.Comparison of simulation and experimental results.

由图9 可知,仿真获得的比冲、效率等性能参数随外加磁场强度的变化趋势与实验结果一致.计算比冲值低于实验结果约20%,其原因是计算比冲由Isp=F/m˙ 计算得来,在本文仿真中未考虑气动力对总推力的贡献,使得计算推力性能低于真实情况,并导致计算比冲低于实验结果.计算效率值相比实验结果有约50%的误差,其原因可能是计算效率与推力和比冲的乘积成正比,计算推力和计算比冲均低于实验结果,导致计算效率误差进一步增大.

虽然推力、比冲和效率的仿真结果因忽略气动加速贡献略小于真实值,但误差不是很大,并且能准确反映工况对性能的影响,故认为,本文建立的模型和方法在高磁场强度工况下有较好的准确性和可信度.

3.4 阴极电流对MPD 工作过程及性能的影响分析

阴极电流是影响MPDT 推力器性能的另一个重要参数,在外磁场B=0.025 T,流量=0.1 g/s工况下,开展阴极电流I=750,1000,1250,1500,1750 和2000 A 情况下的MPDT 工作过程仿真,分析阴极电流对电子数密度分布和推力器性能的影响.

3.4.1 电子数密度分布

不同阴极电流情况下,MPDT 稳态工作时的电子数密度分布如图10 所示.

由图10 可知,推力器内部的电子数密度随着电流的增大而增大;对于B=0.025 T,流量=0.1 g/s 工况,750—2000 A 的阴极电流情况,都在阴极尖端下游位置都出现了电子高密度分布区域.阴极电流增大,高密度区域也逐渐扩大,具体表现为随着电流增大,外磁场对放电等离子体束的约束效果逐渐减弱,低电流工况的等离子体束流发散角明显小于高电流工况,其原因是电流增大,电子数密度增大,电子间的库仑碰撞增多,会使得更多的电子沿跨越磁力线沿径向扩散.

图10 不同阴极电流下的电子数密度分布Fig.10.Distribution of electron number density under different cathode current.

3.4.2 阴极电流对推力器性能的影响

根据不同阴极电流条件下的仿真结果,计算了MPDT 稳态工作时的放电电压U、推力F和效率η,如图11 所示.

图11 MPDT 性能随阴极电流的变化Fig.11.MPDT performance changes with cathode current.

从图11 可知,放电电压随阴极电流的增大而增大,原因是阴极电流越大,阴极表面的电子数密度越大,阴极表面电荷增加,放电电压增大;另外,自感磁场随阴极电流增大而增大,并进一步约束电子向阴极聚集,使得阴极表面电荷增加,放电电压增大.

在750—2000 A 阴极电流范围内,推力随电流增大而近似线性增大.在I=750—2000 A 的电流工况下,计算区域内的自感磁场强度为0.001—0.040 T,量级上与外磁场相当,但自感磁场主要集中于阴极附近,在其余大部分区域外磁场明显强于自感磁场,故此时MPDT 的推力由涡旋加速机制产生,正比于电流大小.

推力器的效率同样表现为随阴极电流I增加而增大的趋势,这与MPDT 低功率工况效率低,高功率工况效率高的实验结果一致.

4 结论

本文建立了基于缩比模型的MPDT 工作过程的粒子仿真模型,研究了MPDT 点火启动时的放电电弧建立过程,以及外磁场强度和阴极电流对MPDT 性能参数的影响,并通过仿真结果与已有实验数据对比,验证了仿真模型的准确性.点火启动过程研究表明,放电电弧建立后,在阴极下游和阳极附近形成两处高电子数密度区并连通以维持稳定的放电电弧,离子主要集中在阴极下游并沿轴向喷射形成等离子体喷流.不同外磁场强度下的结果表明,提高外磁场强度有利于构建阴阳极放电电弧,随着外磁场强度由0.008 T 提升至0.069 T,MPDT 的推力产生机制逐渐由自身场加速转变为涡旋加速,推力和比冲性能随外磁场强度增加近似线性增大.不同阴极电流的结果表明,在B=0.025 T,I=750—2000 A 条件下,在推力器内部绝大部分区域,外磁场强度大于自感磁场,推力产生机制为涡旋加速,推力随电流增加而近似线性增加.MPDT 效率总体表现为随外磁场强度和阴极电流增大而增大的趋势,符合其适用于高功率工况的特性.