1 引言

静电离子推进器和霍尔推进器以其高效率和大比冲优势, 成为目前在轨空间飞行器的主要推进器^[1]。螺旋波放电等离子体源被提出用于推进器中, 是由于其非同寻常的高电离效率:在相近的压力和输入功率下, 等离子体密度比其他放电模式高一个数量级^[2]。螺旋波等离子体源已被研发用于不同的空间电推进系统中^[3], 包括美国NASA JSC研制的可变比冲等离子体火箭(VASIMR)^[4~6]、华盛顿大学研制的高功率螺旋波等离子体推进器(HPHT)、澳大利亚国立大学的螺旋波双层推力器(HDLT)^[7]、麻省理工大学的微型螺旋波等离子体推进器(mHPT)^[8]等。国内大连理工大学从2009年起开展螺旋波等离子体推力器相关研究^{[9, 10]}, 国防科技大学、北京航空航天大学、中国科学院物理所、苏州大学等^[11~14]也都开展了螺旋波放电特性原理研究。

早期的简化螺旋波理论忽略了有限电子质量, 只针对螺旋波的功率耦合作用进行研究, 当考虑电子质量时, 观察到了TG波(Trivelpiece- Gould wave)的存在^[15~19]。Arnush和Chen对圆柱型等离子体内传播的螺旋波开展了大量的研究, 指出, 无阻尼标准模态螺旋波和TG波在强磁场条件下波构型具有明显的差异, 但是在低磁场条件下波构型相似, 因此耦合较强。Arnush提出TG波的作用可能仅对低直流电磁场和轴向波长有着重要影响。低磁场下的标准模态是TG波和螺旋波的混合态。计算结果显示天线向m=+1模式螺旋波的耦合比m=-1模式的更强, 可能是因为m=-1模式的螺旋波分布较窄, 因此在径向边缘处耦合到TG波的振幅较小。Shoji观察到m=+1型螺旋天线在中心轴附近产生了集中的高密度等离子体, 而对于m=-1型螺旋天线, 产生的低密度等离子体充满了放电腔。可以简单地理解为, 射频场向m=-1模式的螺旋波耦合非常弱以至于仅产生了非共振电离^[16]。

对于螺旋波激发高密度等离子体的机理, 已有大量的学者进行研究, 但因其复杂的耦合关系, 至今未得出一致的结论。Chen提出了加速电离电子的朗道阻尼作用^[20], 但是这部分电子数量太少不足以诱发实验中获得的高密度等离子体, 之后Chen和Blackwell否定了这一假设^{[2, 21]}; Shamrai和Taranov提出不同磁场强度下TG波-螺旋波的模式耦合机制, 以及射频功率吸收从体分布到面分布的变化特性^[22]; Blackwell在实验中证明了TG波的存在^[23], Arnush的计算结果支持天线下游远场区域TG波模式的优势^[19]; Lee等提出, 径向密度梯度引起对射频功率的共振吸收, 径向局部螺旋波(RHL)在轴向形成驻波形态, RHL本征模式的频率与射频天线频率相匹配^[24]。

螺旋波放电中的能量沉积机制虽然复杂, 但也是该领域研究的热点。本文立足于螺旋波与TG波线性耦合模式, 探究了不同磁场条件下二者的耦合情况, 以及对波电场、磁场、电流密度的影响规律, 分析了功率沉积密度沿径向和轴向的变化特征。相对于Shamrai的结论^[22], 根据其功率沉积在不同磁场条件下由体分布向面分布的径向变化特点, 给出了功率沉积的轴向分布特征以及本征模式的影响。

2 计算模型 2.1 波场方程

对于径向非均匀等离子体, 麦克斯韦方程的频域表达式为

$\nabla \times \mathit{\boldsymbol{E}}={\rm{i}}\omega \mathit{\boldsymbol{B}}$

(1)

$\nabla \times \mathit{\boldsymbol{B}}={\mu _0}\left( {\mathit{\boldsymbol{j}} - {\rm{i}}\omega {\varepsilon _0}\mathit{\boldsymbol{E}}} \right)= - {\rm{i}}\omega {\varepsilon _0}{\mu _0}\mathit{\boldsymbol{\varepsilon }} \cdot \mathit{\boldsymbol{E}}$

(2)

$\mathit{\boldsymbol{\varepsilon }}=\left[ {\begin{array}{*{20}{c}} S&{ - {\rm{i}}D}&0\\ {{\rm{i}}D}&S&0\\ 0&0&P \end{array}} \right]$

(3)

式中E和B分别为波电场和波磁场, j是天线电流密度, ω是天线驱动频率(即螺旋波频率), ε为冷等离子体介电张量。

$\left. {\begin{array}{*{20}{l}} {S=1 - \sum\limits_\alpha {\frac{{\omega + {\rm{i}}{\nu _\alpha }}}{\omega }\frac{{\omega _{{\rm{p}}\alpha }^2}}{{{{\left( {\omega + {\rm{i}}{\nu _\alpha }} \right)}^2} - \omega _{c\alpha }^2}}} }\\ {D= - \sum\limits_\alpha {\frac{{{\omega _{{\rm{c}}\alpha }}}}{\omega }\frac{{\omega _{{\rm{p}}\alpha }^2}}{{{{\left( {\omega + {\rm{i}}{\nu _\alpha }} \right)}^2} - \omega _{c\alpha }^2}}} }\\ {P=1 - \sum\limits_\alpha {\frac{{\omega _{{\rm{p}}\alpha }^2}}{{\omega \left( {\omega + {\rm{i}}{\nu _\alpha }} \right)}}} } \end{array}} \right\}$

(4)

式中下标α表示粒子种类, ${\omega _{{\rm{p}}\alpha }}={({n_\alpha }{q_\alpha }^2/{\varepsilon _0}{m_\alpha })^{1/2}}$是等离子体频率、${\omega _{{\rm{c}}\alpha }}={q_\alpha }{B_0}/{m_\alpha }$是粒子回旋频率、${\nu _\alpha }$是粒子间碰撞频率。假设等离子体中离子是单电荷, 即${q_{\rm{i}}}= - {q_{\rm{e}}}=\left| e \right|$。根据傅里叶变换, 得到波场耦合微分方程组

$\frac{{\partial {E_\varphi }}}{{\partial r}}=\frac{{{\rm{i}}m}}{r}{E_r} - \frac{{{E_\varphi }}}{r} + {\rm{i}}\omega {B_z}$

(5)

$\frac{{\partial {E_z}}}{{\partial r}}={\rm{i}}k{E_r} - {\rm{i}}\omega {B_\varphi }$

(6)

${\rm{i}}\frac{{\partial {B_\varphi }}}{{\partial r}}=\frac{m}{r}\frac{k}{\omega }{E_\varphi } - \frac{{{\rm{i}}{B_\varphi }}}{r} + \left( {P - \frac{{{m^2}}}{{k_0^2{r^2}}}} \right)\frac{\omega }{{{c^2}}}E$

(7)

$ {\rm{i}}\frac{{\partial {B_z}}}{{\partial r}}=\frac{\omega }{{{c^2}}}{\rm{i}}D{E_r} + ({k^2} - k_0^2S)\frac{{E\varphi }}{\omega }{\rm{ + }}\frac{m}{r}\frac{k}{\omega }{E_z} $

(8)

式中k₀=ω/c, m为角向模数。式(5)~(8)组成了封闭集合, 可以用标准方法数值求解。本文用C++程序Helic对上述波场径向耦合微分方程组数值求解, 结合边界条件, 对于每一个k_z值得到两个独立的波, 求得等离子体中波磁场、波电场、电流密度、相对功率沉积密度等特征参数。该程序相对于ANTENA程序具有更快的计算速度, ANTENNA程序将非均匀等离子体剖分成单层并分别在每一层对边界条件匹配运算。

对于波场的每一个分量V_i, 有

${V_i}={A_1}{v_{1, x}} - {A_2}{v_{2, x}}$

(9)

式中A₁和-A₂为待确定幅值, x表示z, r, Φ三个方向分量。结合边界条件, 可以求得这些代数方程的解, 其具体表达式为

${V_i}\left( r \right)=\frac{{{H_2}{v_{1, x}}\left( r \right) - {H_1}{v_{2, x}}\left( r \right)}}{{{F_1}{G_2} - {F_2}{G_1}}}{\rm{i}}{\mu _0}{K_\phi }$

(10)

其中

${F_n}={b_{n, r}}\left( a \right) + {\rm{i}}\frac{k}{{{T_0}}}{p_m}\left( a \right){b_{n, z}}\left( a \right) + \frac{m}{{aT_0^2}}{\mu _0}\left[ {\omega {\varepsilon _0}{e_{n, z}}\left( a \right)} \right]$

(11)

$\begin{array}{*{20}{l}} {{G_n}={j_{n, r}}\left( a \right) + {\rm{i}}\frac{m}{a}\frac{{k_0^2}}{{T_0^2}}\frac{1}{{{\mu _0}}}{b_{n, z}}\left( a \right) + }\\ {\omega {\varepsilon _0}\left[ {\frac{k}{{{T_0}}}{q_m}{e_{n, z}}\left( a \right) - {\rm{i}}{e_{n, r}}\left( a \right)} \right]} \end{array}$

(12)

${H_n}=\frac{{kb}}{{{T_0}a}}{p_m}\left( b \right){\mu _0}{G_n} - \frac{m}{a}\frac{{k_0^2}}{{T_0^2}}{\rho _m}{F_n}$

(13)

式中角标n为1时代表螺旋波模式, n为2时代表TG波模式。b_{n, x}, e_{n, x}, j_{n, x}表示波磁场、波电场和电流密度的基函数。

对于半径为a, 长度为2l的圆柱型等离子体放电区域, 功率吸收相对密度沿径向和轴向的分布计算表达式如下

$\begin{array}{*{20}{l}} {{P_r}\left( {r, m} \right)=\frac{{\omega {\varepsilon _0}}}{2}\mathop \int_{ - l}^l {\rm{Im}}\{ S\left( {\left| {{\rm{i}}{E_r}{|^2} + } \right|{E_\phi }{|^2}} \right) + }\\ {P|{E_z}{|^2} + 2D{\rm{Re}}\left\{ {{\rm{i}}E_r^{\rm{*}} \cdot {E_\phi }} \right\}\} {\rm{d}}k} \end{array}$

(14)

$ \begin{array}{*{20}{l}} {{P_z}\left( {z, m} \right)=\frac{{\pi \omega {\varepsilon _0}}}{2}\int_0^a {{\rm{Im}}\{ S\left( {\left| {{\rm{i}}{E_r}{|^2} + } \right|{E_\phi }{|^2}} \right) + } }\\ {P|{E_z}{|^2} + 2D{\rm{Re}}\left\{ {{\rm{i}}E_r^{\rm{*}} \cdot {E_\phi }} \right\}\} r{\rm{d}}r} \end{array} $

(15)

2.2 参数设置

计算模型如图 1所示, 等离子体位于内径为a的圆柱型石英管内, 缠绕于石英管外的螺旋波天线内径为b且忽略其厚度, 整个放电装置置于内径为c的筒形金属仓内。驱动频率、轴上等离子体密度、碰撞因子等参数取螺旋波放电的典型值, 选取Ar为工质气体, 天线类型为半波螺旋天线(Half helix antenna), 螺旋波角向模数m=+1, 具体参数以及计算区域尺寸见表 1。假设电子温度以及气体压强沿径向分布均匀, 等离子体密度沿轴向均匀, 沿径向非均匀, 分布情况如图 2所示, 服从三参数(s=2, t=1, f_a=0.1)方程^{[19, 25]}

$\frac{{{n_\alpha }\left( r \right)}}{{{n_{\alpha 0}}}}={\left[ {1 - {{\left( {\frac{r}{\omega }} \right)}^s}} \right]^t}, \;\omega =\frac{a}{{{{\left[ {1 - f_a^{1/t}} \right]}^{1/s}}}}$

(16)

3 结果和分析

在不同的磁场条件下, 由于螺旋波和TG波所受阻尼不同, 且径向波长差异较大, 其在等离子体内沿径向的传播距离或者功率耦合范围会受较大影响。TG波所受阻尼大, 高磁场下径向波长较短, 很难向等离子体中心处传播, 主要作用在径向边界处的功率耦合^{[16, 17]}; 而螺旋波所受阻尼小, 在等离子体内的穿透能力强, 可沿径向和轴向传播较远距离。本文计算发现, 随着磁场的增大, TG波和螺旋波耦合模式发生转变, 螺旋波的本征模式对功率的沉积作用逐渐加强。

3.1 波磁场分布特点

从图 3中可以看出, 射频磁场角向分量B_phi在径向边界处为0, 满足边界条件。B_phi和B_z沿径向分布的最大值主要集中于近轴处(r=0m, r=0.02m)。随着磁场的增大, |B_z|峰值有降低的趋势, 而|B_phi|的峰值略有波动。B_r沿径向分布特征与B_phi较为接近, 且幅值相当, 本文不再列举。

r_p表示诊断点的径向位置, 即在径向r=r_p处求解方程组计算波场的轴向变化特征, 如图 1中虚线箭头所示。在诊断位置为r_p=0.02m和r_p=0.048m处, 对射频磁场分量B_phi进行对比分析, 其分布特点如图 4所示。在r_p=0.048m处(图(b)), B_phi沿z轴的分量在天线位置处存在较强的峰值, B₀≤500G时沿z轴逐渐衰减, 当磁场强度B₀增加到700G以上时, B_phi在天线下游沿z轴的分布逐渐趋于均匀且呈周期形态。由于射频磁场在径向边界处的作用较弱, 在此对r_p=0.02m的分布状态进行了对比, 同样, 当B₀从低于500G增加到高于700G时, B_phi由衰减态的振动分布转变为准周期形式的简谐振动分布, 且幅值增大。

3.2 电流密度分布特点

从图 5中可以看出, 低磁场条件下, 电流密度角向分量J_phi沿径向分别在r=0的中心轴和径向边界处取得最大值, J_z的径向分量则在r=0.015m左右以及径向边界处取得最大值。在r=0.045m附近, J_phi和J_z达到最小值, 之后向边界处迅速增大。对比不同磁场B₀下的分布特点, 可以看出, 随着磁场强度的增加, J_phi和J_z在径向的分布优势向边界处转移, 二者的幅值在径向边界处迅速增大而在内部则逐渐减小。

B₀=200G, 300G时, |J_z|在天线末端出现明显的峰值, 并且沿z轴迅速衰减。当磁场增加到400G和500G时, 出现了多峰值现象, 且在远离天线的位置处, 有向周期波动发展的趋势, 如图 6所示。随着磁场强度的进一步增加, 沿z轴的波动状态进一步加强, 从图 6(a)中可以看出, 磁场强度达到700G以上时, |J_z|呈准周期分布, 其波动周期近似相等, 只是幅值较为不同。J_z的实部和虚部分量如图 6(b)所示。J_phi沿z轴的分布特征以及随磁场的变化趋势与J_z相同。

3.3 波电场分布特点

E_r沿径向的分布特征为, 0 < r < 0.045m范围内, 幅值变化较为缓慢, 当r≥0.045m时迅速增大达到峰值, 且幅值远远高于近轴处(r < 0.045m), 如图 7所示。随磁场强度增加, 径向边界处的峰值迅速增大。E_phi与E_z幅值低于E_r近两个数量级。

根据E_r的径向分布特征, 选取r_p=0.048m的诊断位置, 对E_r沿轴向的分布进行分析, 正如前文提到的, r_p=0.02m处的分布特征以及随磁场的变化特点与r_p=0.048m诊断位置处相似, 只是幅值较低, 这也与沿径向的分布特点相吻合。

从图 8中可以看出, 低磁场条件下E_r的幅值在天线下游近距离处形成峰值, 之后沿z轴逐渐衰减, 当B₀=500G时, 这种衰减形态已不明显, 同时从远离天线位置的z=0.6m开始, |E_r|每隔0.45m有一个幅值增减的周期。当B₀≥700G时, 无论是实部或虚部的分量, 还是绝对值, 都呈现出了准周期分布状态。另外, 随着磁场强度的增加, E_r的幅值逐渐增大。

3.4 功率变化特点

外加磁场B₀=200G, 300G时, 在螺旋波和TG波的耦合作用下, 径向出现两个功率吸收峰值, 分别为中心轴处和径向边界处, 如图 9(a)所示。随着外加磁场强度的增加, 中心轴处功率沉积密度逐渐减小, 而径向边界处的密度逐渐增加, 功率沉积密度由中心轴向径向边界处转移, 当磁场B₀≥700G时, 径向边界(r > 0.045m)处的功率耦合密度已经远远高于其他位置, 占据了绝对优势。这种变化趋势与Chen和Shamrai的表述一致。Shamrai提出, 功率沉积机制为射频天线将输入功率耦合到螺旋波, 由螺旋波将功率线性耦合到等离子体和边界处的TG波中, 即功率耦合的要素首先是螺旋波与TG波在边界处的线性模式转化, 其次是TG波的强阻尼作用。这种机制能够导致近壁面处能量沉积的增强, 而不是在B_z²达到峰值的径向位置处。

高磁场下射频功率沉积主要集中于径向边界处, 而实验测得电子数密度通常呈中心轴处高、径向边界处低的分布形式, 本文假设的等离子体密度径向分布满足该特征。Chen等^[26]认为这种现象的成因是端板处的Short-circuit效应, 该效应使得理论上沿磁感线螺旋运动的电子实际沿径向垂直于磁场方向运动, 之后电子服从玻尔兹曼关系, 形成电场, 驱动在径向边界处激发的离子向内运动, 当中心密度峰值形成后电场反向。

从图 10中可以看出, 功率沉积沿轴向分布特征随磁场强度变化非常明显。当外加磁场B₀=200G~400G时, 功率在天线下游近距离内耦合较多, 在z=0.43m附近出现了明显的峰值, 且幅值沿z轴迅速降低; 当磁场强度达到500G时, 已无较为明显的峰值。同时可以看出, 随着磁场强度的增加, 功率沉积密度的幅值沿z轴的减小趋势变缓, 且曲线波动变大。当外加磁场强度达到700G以上时, 功率沉积密度沿轴向的分布特征发生了根本性的改变。功率沉积密度从天线中心z=0.2m处迅速增大达到峰值, 随后沿z轴以呈准周期性波动趋势发展, 幅值表现为缓慢的衰减趋势, 如图 10(b)所示。

Chen提出^[16], 低磁场条件下, TG波和螺旋波的构型相似, 因此耦合较强, 并且TG波沿径向的传播距离大, 图 9(a)中可以明显观察到径向的两个吸收峰值; 从轴向功率沉积特征可以看出, 天线末端区域功率沉积较多, 这是因为螺旋波将大部分功率耦合到TG波中, 同时TG波的阻尼较强, 沿轴向的传播距离短, 因此在靠近天线处功率耦合出现明显的峰值, 且沿z轴迅速衰减。随着磁场强度的增加, 螺旋波耦合到TG波的功率降低, 因此螺旋波将一部分功率输运到下游, 在传播的过程中持续向TG波耦合, 这样沿轴向的功率沉积密度就表现为螺旋波的本征模式, 即图 10(b)中沿轴向的准周期分布。随着磁场强度的增加, TG波被限制在径向更薄的区域内, 基本上算是表面波^[22], 所以导致功率沉积密度向径向边界处积累, 并远高于内部。

随磁场强度增加, 波磁场、波电场和电流密度的轴向分布由衰减趋势变为周期波动趋势, 螺旋波在等离子体内的轴向传播能力增强, 其反射波与入射波形成共振的本征模式, 影响着TG波和螺旋波的耦合, 进而形成这种沿轴向的功率沉积准周期形式。从图中可以看出, 磁场增加的同时, P_z在下游区域的波动频率发生微弱变化。|J_z|, |E_r|沿z轴的分量与P_z随磁场的变化特征很接近。

4 结论

针对不同磁场条件下的螺旋波-TG波耦合模式进行了数值模拟研究, 分析了磁场对功率沉积、波电场、磁场和电流密度的影响。获得主要结论如下:

(1) 随着外加磁场强度的增大, 波磁场、电场的部分分量沿z轴的分布由幅值衰减状态变为准周期性波动状态, 螺旋波的轴向传播能力增强, 表现为轴向的本征模式, 电流密度的变化特征与功率沉积密度较为相似。

(2) B₀≤500G的低磁场条件下, 螺旋波与TG波的构型相似, 线性耦合较强, 螺旋波将大部分射频功率耦合到TG波内, 并由TG波的强阻尼作用, 在天线下游0.2~0.4m距离内沉积到等离子体中, 该条件下TG波在径向的传播距离远, 相对功率沉积在径向有两个吸收峰值。

(3) B₀≥700G的高磁场条件下, 螺旋波向TG波的耦合效率降低, 螺旋波将一部分射频能量输运到下游并持续向TG波耦合, 由TG波的阻尼作用沉积到等离子体中, 轴向的功率沉积密度分布特征就表现为螺旋波的本征模式, 即沿z轴的准周期性分布, 而径向功率沉积集中于边界r > 0.045m范围内。