成玉国 夏广庆

1)(中国人民解放军91550部队91分队,大连 116023)

2)(大连理工大学,工业装备结构分析国家重点实验室,大连 116024)

1 引 言

目前预研的具有较大应用前景的无电极烧蚀、无羽流中和器的感应式脉冲等离子体推力器中,一类以FARADAY为代表[1],其将无电极等离子体源(如螺旋波源等[2,3])作为预电离级,在源室下游进一步利用脉冲电流圈加速;另一类以PIT为代表[4]的推力器则直接利用特殊构型的线圈,如图1(a)所示,电离其表面的中性气体在等离子体中形成环形电流片,与线圈排斥产生推力,实现中性气体的电离和加速.这种推力器中,中性气体通过脉冲流量阀的引导,均匀地分布在线圈表面,电容器能量经线圈激励产生电磁场,传递至气体,使其电离和加速(图1(b)).后一类脉冲工作方式具有推力器结构简单、方便调节、比冲高等优点,对其工质加速过程和推力产生机理的研究仍处于不断深化和发展之中.

图1 (a)PIT实物图;(b)PIT结构简图[5]Fig.1.(a)The picture of the actual object of the PIT;(b)the Simpli fied con figuration of the PIT[5].

Jahn[6]提出上述第二类推力器推进应用中脉冲气体电离和加速的两个主要问题:1)气体击穿特征时间与脉冲电场回路放电时间的匹配问题;2)如何在气体运动远离线圈表面之前实现绝大部分脉冲能量的注入,以提高电离和加速效率的问题,并进一步提炼出线圈等离子体解耦距离这一影响加速的主要参数.

针对上述问题,研究人员分别发展了线圈等离子体脉冲感应耦合机电模型和等离子体加速的磁流体力学模型.机电模型实现了线圈回路和等离子体负载回路通过感应线圈的耦合[5].Polzin和Choueiri[7]总结出动态匹配系数α,将回路放电和等离子体加速过程联系,分析了等离子体的解耦距离.在近似正弦的脉冲电流中,等离子体在前1/4周期内达到解耦距离时,等离子体可实现高效率的电离和加速.

二维磁流体数值研究中,Mikellides等[8]利用MACH2计算发现了在气体质量＞3 mg时,推力器的效率保持稳定.这一结果与实验中观测到的临界质量现象一致.在一些实验中,测量结果表明PIT推进效率随着能量的变化而改变[9],且Mikellide等[10]将MACH2与电路模型结合的计算进一步表明了超过临界质量之后效率的下降是由于内能的沉积所致.这些矛盾的结论表明,感应式脉冲加速推力效率这一问题值得进一步的探索.已有的模拟中,对于影响推力产生的等离子体运动过程、运动参数的空间分布以及磁场强度变化对等离子体加速和电离过程的影响研究较少.

本文研究感应式脉冲等离子体的电离和加速过程,计算在单脉冲、短放电周期(10µs)条件下等离子体流场参数分布及其特性,并着重分析脉冲强度变化时等离子体的加速特性和推进性能的变化情况.

2 感应式脉冲等离子体流动控制方程

2.1 物理模型

在µs脉冲作用下,等离子体中的感应电磁场变化较快,电子与重粒子之间的耦合强烈,形成整体的等离子体团运动,采用文献[11—13]处理等离子体组分的方法,对等离子体在脉冲电磁场激励下的运动采用单流体磁流体力学描述,认为电子和离子处于热力学平衡状态.其连续、动量和能量方程分别为:

连续方程

式中,ρ为等离子体密度,V为等离子体速度;动量方程

式中,B为等离子体流场中的磁场强度,p为等离子体压力,J为等离子体电流,J=∇×B/µ0,µ0为真空磁导率;

能量方程

式中,et为总比能,,q为等离子体热流,σe为等离子体电导率,f(p,ρ)为真实气体状态函数,其计算方法在2.4节中讨论,Qrad,tn为辐射能量损失,采用文献[14]的处理方法.

对等离子体流场中的电磁场运动采用参考文献[15]的处理方法,用如下形式的磁扩散方程关联磁场和速度场:

式中,∇×(V×B)=∇·(V B−BV).

为了解决以上方程组中的无黏项Jacobi矩阵不满秩和磁场散度由于离散造成的不为零的问题,本文采用文献[16,17]给出的双曲型散度清除代数方程,其已经在经典算例中得到了验证.方程如下:

式中,ψ为广义拉格朗日乘子;ch,cb分别为控制∇·B误差的对流和耗散系数.

2.2 广义拉格朗日乘子形式的磁流体动力学(MHD)方程组

由以上方程,得到广义拉格朗日乘子形式的MHD方程组的二维轴对称形式为

式中,µ为黏性系数.

以上方程组由等离子体的状态方程(11)封闭,

式中,ne为电子数密度;na为中性粒子数密度;j=1,2···,表示离子电离阶数;T为等离子体温度,与电子温度Te相同;kB为玻尔兹曼常数.

2.3 定解条件

2.3.1 气动边界条件

图1所示推力器的计算简化结构如图2所示,其主要由电路和等离子体部分组成,电容器经线圈激励在其表面产生高频振荡电磁场,并逐步渗透至等离子体内部.文中计算主要关注等离子体部分,电容器放电对等离子体的影响在计算中体现为:在等离子体边缘(线圈表面,AD边)施加时变磁场.气体流经阀门,在线圈表面沉积,在脉冲电流线圈的作用下,加热并电离,向近似真空的环境中膨胀加速.等离子体的入流/出流、壁面以及轴对称边界等三类边界条件根据文献[18]的讨论给出.

图2 (网刊彩色)感应式脉冲等离子体计算域构型Fig.2.(color online)Calculation con figuration of the inductive pulsed plasma.

2.3.2 磁场边界条件

入口处磁场的条件由电容器的放电频率和强度决定,假设其为正弦变化,线圈电流变化为J(t),径向上的磁感应强度Br随时间变化为

即Br(t)=B0sin(ωt),B0为峰值强度,ω=2πf,f为放电频率,文中设置放电周期为10µs,router,rinner为分别为线圈外径和内径.线圈表面的磁场主要在径向,在角向为零、径向均匀的情况下,忽略轴向磁场,此时表面的磁场散度为零.

文中计算的放电波形为正弦变化,根据COMSOL软件的计算,对于2.3节中所述的线圈的几何构型,单周期正弦脉冲能量和峰值磁场强度B0的关系如图3所示.

图3 单脉冲正弦能量Esingle-pulse与峰值磁场强度B0关系曲线Fig.3. The relationship between single sinusoidal pulse energy Esingle-pulseand peak field intensity B0.

壁面处的磁场采用无穿透条件.流场下游区域足够大,认为在脉冲作用的时间内,尚未影响下游的计算区域,边界处的轴向磁场Bz(t)采用外推条件.

实际计算中,发现在脉冲作用周期内,在电流圈斥力的作用下,等离子体轴向运动距离在1—2个线圈半径的长度.因此,为了突出主要运动区域,研究其中的流动机理,文中径向起算点为A,B点,轴向将AB线段设置为壁面条件,暂不考虑拐角处热膨胀对流动的影响.

对于文中的非定常MHD流动,空间格式应用M-AUSMPW+[15,19],采用结构网格离散,在时间上采用具有真实垂直深度性质的三阶Runge-Kutta法.文中计算的构型中,线圈内径为0.05 m,外径为0.2 m,轴向方向取5倍内径长度.文中对单脉冲周期中等离子体的电离、加热以及加速过程开展研究.

2.4 热力学性质

计算采用Ar气,高温情况下,其将会出现三阶以上电离.本文考虑氩气的五阶电离.文中采用的简并度及其能级数值由美国标准与技术研究院数据库(National Institute of Standards and Technology,NIST)给出.

混合组分的比焓表达式为[20]

式中,zj,Ej为各能级配分函数和电离能.

比焓与比内能的关系为

由(13)式,等离子体的比焓由重粒子分量、电子分量、电子激发能分量以及电离分量构成.比焓计算中对文献中报道较少的温度范围(2.5—5.0 eV)进行了拓展,如图4所示,在低温范围内(0—2.5 eV),结果与文献符合[20].

图4 100.0 Pa条件下Ar等离子体比焓随温度变化Fig.4.The speci fic enthalpy of Argon plasma varies with Tefor P=100.0 Pa.

3 计算结果与分析

3.1 脉冲作用周期内等离子体的运动分析

图5(a)—(c)中给出了在峰值磁场强度为0.1 T情况下,单脉冲周期内等离子体流场参数随时间的变化.可以看出,随着时间的推进,等离子体的能量和速度逐渐增加,且在全部脉冲作用时间内等离子体的速度始终在增加.

由等离子体在各时刻的状态参数分布可以看出,在脉冲电流的作用下其运动图景大致描述如下.

放电初始阶段,在感应线圈产生的时变磁场作用下,等离子体中产生感应电磁场,电场强度达到气体电离阈值后,中性气体电子被剥离、加速,并进一步撞击中性气体.线圈表面等离子体迅速被加热、升温并向外膨胀,在瞬变磁场的持续作用下,表面的高速、相对高密度气体向前方运动(图5(a),(b)),同时,这一过程伴随着线圈的磁场能量在等离子体中的逐步渗透,且流场中的磁场强度及其影响范围也在逐渐增大.

随着放电和加速过程的持续进行,达到1/2周期以后,线圈电流的方向发生改变.气压足够大或者峰值磁场强度较小时,反向作用周期内表面残余气体可获得进一步加热,与前方等离子体团的电流片异向,等离子体中前后异向的电流片相互吸引,从而降低前向等离子体团的速度,并减弱等离子体的总体加速效果.

图5(c)给出了等离子体的磁场随时间的变化.脉冲电流超过1/4周期后,电流方向不变,但其变化率与前1/4周期相反.这段时间内,中性气体在脉冲电磁场作用下进一步电离,表面的高温等离子体向前缘移动,但气体密度较高的情况下,线圈表面将滞留带电粒子和中性气体的混合物.放电脉冲进入后半周期后,这部分气体受负变化率的影响产生反向等离子体电流环,这种新产生的电流环一方面对前缘电流环具有吸引作用,另一方面阻止了激励磁场进一步向前缘的渗透.

3.2 不同单脉冲能量下等离子体加速过程研究

图6给出了线圈峰值磁场强度为0.45 T时等离子体速度峰值时刻流场的参数分布(密度、轴向速度、径向磁场、温度、电子数密度和电流密度).在强脉冲电磁场的作用下,将有更多的电磁场能量用于加速等离子体.从图6的结果可以看到,此时流场运动与壁面黏滞阻力特征时间差异较大,尤以密度和电流的分布较为明显,表明在较大的磁场强度下,壁面黏性对等离子体造成的影响减小.

图6中,线圈表面等离子体密度与背景气体几乎相同,从电子数密度的分布图中可以看出,此时表面附近的电子数密度较前缘峰值降低了约1—2个量级(图6中的密度分布),但由于这一区域获得的电磁场能量最大,等离子体温度为流场的高值区域,如图6中的温度分布所示.在较高磁场强度下,放电一段时间后,表面的等离子体粒子数密度是极低的.因此,虽然温度较高,但基本不具备进一步产生可观冲量的价值.

图5 (网刊彩色)0.1 T时不同时刻等离子体流场参数变化 (a)密度随时间变化的分布;(b)速度随时间变化的分布;(c)等离子体磁感应强度随时间变化的分布Fig.5.(color online)The variations of the flow parameters at different time for 0.1 T:(a)The density distribution as the time is varied;(b)the velocity distribution as the time is varied;(c)the plasma inductive field intensity distribution as the time is varied.

图6 (网刊彩色)峰值磁场强度为0.45 T等离子体流场参数分布Fig.6.(color online)The plasma parameter distribution for peak field intensity of 0.45 T.

图7 (网刊彩色)峰值磁场强度为0.45 T时不同时刻等离子体速度的分布Fig.7.(color online)The plasma velocity distribution varies with the time for the peak field intensity of 0.45 T.

从图6中磁场和等离子体电流片的分布可看出,磁场的有效作用范围受到等离子体密度分布影响,加速过程中,远离线圈表面时,等离子体中磁场逐渐减弱.

图7中给出了不同时刻流场中等离子体速度随时间的变化(从零加速至峰值以及之后的减速).从图中可以看出,速度相对于3.1小节中的低强度磁场提升了约1—2个量级,其最大值出现于放电的1/4与1/2周期之间.随着放电时间的推移,等离子体团速度逐渐增加,但受到磁场与流体运动相位滞后以及惯性的影响,越过电流峰值之后,等离子体速度在1/4周期后达到最大.

从图7可以看出,在较高的脉冲能量情况下,由推力器几何构型决定的线圈等离子体解耦距离对等离子体加速影响显著,表明了高脉冲输入能量可获得较高的加速效率,但加速过早地达到解耦距离,也限制了能量的进一步有效利用.

图8给出了根据电感计算理论得到的线圈结构与解耦距离的关系[21,22],对于文中内径0.05 m、外径0.2 m的线圈,其解耦距离约为0.04 m,这与0.45 T峰值磁场情况下等离子体在1/4周期内运动距离相近,表明在这一脉冲能量下,等离子体的运动开始达到解耦距离.

图8 (网刊彩色)解耦距离与线圈内外径的关系Fig.8.(color online)The relationship between the decouple distance and the coil outer-inner diameter.

图9 不同时刻等离子体放电参数沿轴向的分布(0.45 T)Fig.9.The axial plasma discharge parameters vary with the time for peak intensity of 0.45 T.

图10 不同时刻等离子体密度沿轴向的分布(0.45 T,径向截面平均)Fig.10.The axial distribution of the plasma density as the time is varied(0.45 T,radial cross section averaged).

图11 (网刊彩色)不同磁场强度的峰值速度时刻电流密度分布Fig.11.(color online)The current density distribution of different field intensity at the time of peak velocity.

由上文中峰值强度0.1 T和0.45 T的结果对比可以看出,磁场强度较低的情况下,等离子体速度低,加速距离短,脉冲作用时间内可始终获得加速;但电磁场能量低,等离子体的能量以及感应强度有限,加速性能提高有限.磁场强度增大后,一部分能量增加等离子体内能,另一部分增加等离子体团的动能,于是等离子体中诱导电流片吸收的能量增加[23],从而导致等离子体团推进速度增加.主等离子体的运动距离超出解耦距离后,受脉冲电流及磁场的影响减弱,同时部分能量传递通过传导以及辐射等方式传递给中性气体,且随着磁场强度的增强,诱导电流片速度进一步增加,与电流线圈解耦的时间缩短.

图9给出了不同时刻等离子体放电参数沿轴向的分布(参数对径向进行了平均).随时间变化的曲线表明,等离子体不断向前推进过程中,受影响的气体区域不断扩大.可以看出磁场的影响区域明显大于等离子体流场的扰动区域,这表明磁场在等离子体中有效渗透,有利于等离子体的电离.

图10给出了不同时刻等离子体数密度沿轴向的分布,密度对径向进行了平均.结果表明,随着时间的发展,等离子体数密度逐渐增加,高阶电离粒子数密度组分增大,线圈表面附近等离子体几乎完全电离.由上文的结论,密度的这一分布结果主要由于表面是加热的主要区域造成.

图11和图12给出了不同峰值磁场强度下最大速度时刻等离子体流场中电流密度和流体密度的分布.从图中可以看到,磁场强度越高,等离子体获得的能量越多,其加速越显著,磁场强度提高的极限是在放电经过峰值时.进入负半周期之后,由于没有中性气体的补充,等离子体团远离耦合距离,这段放电时间的能量几乎得不到合理的利用.因此,提高磁场强度可能会造成能量利用效率增长的缓慢.从图13速度的变化趋势可看出,磁场强度增大后,速度幅值增量逐渐减小.

图12 (网刊彩色)不同磁场强度的峰值速度时刻密度分布Fig.12.(color online)The plasma density distribution of different field intensity at the time of peak velocity.

图13 不同峰值磁场强度的最大激励速度变化曲线Fig.13.The maximum velocity varies as the peak filed strength is changed.

3.3 比冲和推力效率分析

图14给出了感应式脉冲等离子体比冲和推力效率随脉冲能量的变化,结果表明,在前半周期运动能够达到解耦距离的情况下,等离子体的比冲和效率将有较大的提升.比冲从低脉冲(0.1 T)情况下约700 s提高到高脉冲能量(0.55 T)的约为2000 s.感应式加速中能量的提高对推进性能的提升显著,显示出感应式脉冲加速的优势,尤其是在四分之一周期电离充分的情况下,等离子体的加速效果明显.

图14 比冲和推力效率随单脉冲能量的变化Fig.14. The speci fic impulse and thrust efficiency varies as the pulse energy is changed.

4 结果与讨论

本文建立了二维轴对称的感应式脉冲等离子体作用的磁流体力学模型,针对计算中出现的高温情况,对Ar等离子体的热力学性质进行了计算,着重分析了不同脉冲电流能量情况下等离子体的加速特性、磁场形态的变化.

计算中发现在峰值较低的磁场强度(～0.1 T)情况下,电容器的放电能量主要用于中性气体的电离以及气动加热,造成感应电磁场对粒子的加速效果有限,等离子体吸收的能量不足以支持高阶电离,从而不能形成强度较高的感应电流片,与线圈之间的斥力较弱,因此,等离子体在低强度脉冲电流作用周期内,速度始终较低.

当峰值磁场强度增大(＞0.45 T)时,推进工质以及压力不变的情况下,更高比例的输入脉冲能量将被用于加速等离子体,剧烈变化并且幅值较大的电磁场有利于等离子体重粒子和电子之间的耦合,提高了等离子体感应电流片的电流密度和强度,可在等离子体-线圈之间产生较大的斥力,增加等离子体的速度,使得脉冲加速的效果有较大的提升.

磁场强度较高的情况下,等离子体速度可加速至20 km/s以上,从文中的计算结果可以看出,等离子体运动达到解耦距离的时间缩短,同时伴随着等离子体高阶电离、能量辐射、磁场影响迅速扩大等非线性过程,且随着强度的不断增大,等离子体速度提高,这一时间逐渐减小.获得这段时间内各种物理过程对能量吸收和加速的影响,以及等离子体的解耦距离与磁场之间的关系,是下一步提高感应式脉冲等离子体推进性能的重要研究方向.

[1]Choueiri E Y,Polzin K A 2006J.Prop.Power223

[2]Cheng Y G,Cheng M S,Wang M G,Li X K 2014Acta Phys.Sin.63035203(in Chinese)[成玉国,程谋森,王墨戈,李小康2014物理学报63035203]

[3]Cheng Y G,Cheng M S,Wang M G,Li X K 2014Chin.Phys.B23105202

[4]Mikellides P G,Neilly C 2007J.Prop.Power231

[5]Polzin K A 2011J.Prop.Power273

[6]Jahn R G 1968Physics of Electric Propulsion(New York:McGraw-Hill)p268

[7]Polzin K A,Choueiri E Y 2006IEEE Trans.Plasma Sci.343

[8]Mikellides P G,Villarreal J K 2007J.Appl.Phys.10210

[9]Dailey C L,Lovberg R H 1993NASACR-1993-191155

[10]Mikellides P G,Turchi P J,Roderick N F 2000J.Prop.Power165

[11]Xie Z H 2013Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[谢泽华2015博士学位论文(长沙:国防科学技术大学)]

[12]Polzin K A,Sankaran K,Ritchie A G,Reneau J P 2013J.Phys.D:Appl.Phys.46475201

[13]Tian Z Y 2008Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[田正雨2008博士学位论文(长沙:国防科学技术大学)]

[14]Li X K 2011Ph.D.Dissertation(Changsha:National University of Defense Technology)(in Chinese)[李小康2011博士学位论文(长沙:国防科学技术大学)]

[15]Kim S,Soogab L,Kyu H K 2008J.Comput.Phys.2278

[16]Dedner A,Kemm F,Kröner D,Munz C D,Schnitzer T,Wesenberg M 2002J.Comput.Phys.1752

[17]Dedner A 2003Ph.D.Dissertation(Freiburg im Breisgau:Alert Ludwigs-Universität Freiburg)

[18]Yan C 2006Computational Fluid Dynamic(Beijing:Beihang University Press)pp254–258(in Chinese)[阎超2006计算流体力学方法及其应用(北京:北京航空航天大学出版社)第254—258页]

[19]Kim K H,Kim C 2005J.Comput.Phys.2082

[20]Chen X 2009Thermal Plasma Heat Transfer and Flow(Beijing:Science Press)pp48–59(in Chinese)[陈熙 2009热等离子体传热与流动(北京:科学出版社)第48—59页]

[21]Kalantarov P L,Neitlqy L A(translated by Chen T M,Liu B A,Luo Y L,Zhang Y H)1992 Inductance Calculation Handbook(Beijing:Mechanism Industry Press)pp1–9(in Chinese)[卡兰塔罗夫P L,采伊特林L A 著 (陈汤铭,刘保安,罗应力,张奕黄 译)1992电感计算手册(北京:机械工业出版社)第1—9页]

[22]Che B X 2015M.S.Thesis(Changsha:National University of Defense Technology)(in Chinese)[车碧轩2015硕士学位论文(长沙:国防科学技术大学)]

[23]Polzin K A,Choueiri E Y 2006IEEE Trans.Plasma Sci.343