张凯宇-AFMPDT涡旋流动解析研究

　AF-MPDT 涡旋流动的解析研究 张凯宇，汤海滨，王一白 （北京航空航天大学）

　摘要：洛伦兹力驱动的磁流体角向涡旋流动是附加磁场磁等离子体推进器（AF-MPDT）内一种重要流动现象，先前的研究认为这种涡旋流动是 AF-MPDT 内部重要的能量储存形式和主导加速机制。但受限于磁流体偏微分方程的求解困难，针对 AF-MPDT 涡旋流动的解析研究只停留在一维，而二维的研究仅有数值仿真。本文基于磁流体力学偏微分方程组和 Sturm-Liouville 理论，用分离变量法结合压降平衡积分给出了同轴等长圆柱电极的 AF-MPDT 内的流场-电场自洽的涡旋流动的二维的解析解。该解析解被表达为无穷级数的形式，分析表明级数收敛迅速，收敛速度只与通道的内外径之比有关，且该解析结果与数值仿真结果吻合良好。该解析解为 AF-MPDT 涡旋加速的计算提供了一个高效的方法，为磁流体涡旋流动的机理研究提供了有力的支撑。

　关键词：大功率电推进器，AF-MPDT，涡旋流动，磁流体力学，流动解析解 1. 引言 AF-MPDT 凭借其高比冲、大推力以及容易实现大功率和小型化结合的特点，在深空探测任务中有巨大优势。AF-MPDT 通道内轴向附加磁场的存在和同轴电极的放电结构带来了角向的洛伦兹力 j r B z ，洛伦兹力驱动等离子体在角向上涡旋运动。这种放电通道内的涡旋运动对 AF-MPDT 中的比冲形成和功率构成有着重要的影响。首先在比冲方面， AF-MPDT有四种等离子体加速机制：气动加速、感生磁场加速、霍尔效应加速、涡旋加速。涡旋加速就是指利用磁喷管效应，把放电通道内积累的角涡旋流动的能量转移成轴向的比冲。Tang的仿真计算表明涡旋加速在这四种加速机制中贡献占比高达 83 % [1] ；Choueiri 认为涡旋加速是 AF-MPDT 的主导加速机制 [2] 。其次在功率方面，AF-MPDT 的电极通道内有四种电压构成：等离子体电阻电压、霍尔效应电压、涡旋运动和附加磁场诱发的动生电动势、轴向运动和感生磁场诱发的动生电动势。Turchi 的仿真计算表明由涡旋运动带来的动生电动势在总电压中占比最高可达 70 % [3] ；Fradkin 认为涡旋运动是 AF-MPDT 通道内重要的能量储存形式[4] 。故无论从比冲还是功率看，等离子体的角向涡旋流动在 AF-MPDT 中都扮演十分重要的角色，具有重大的研究价值。

　此前研究者基于磁流体动力学针对 AF-MPDT 的涡旋流动提出了多种解析模型和数值模型，逐渐加深了对这种涡旋流动现象的认识。代数解析模型方面，Fradkin [4] 与 Albertoni [5]对涡旋动量矩采用径向积分平均的方法，得到了角涡旋速度的代数的解析模型，该模型所预测出的涡旋速度在数量级上和实验结果 [6] 较好吻合。一维解析模型方面，Sasoh 认为通道内的等离子体类似刚体转子在洛伦兹力作用下旋转，旋转角速度与径向位置无关，忽略了粘性力作用，得到涡旋速度的径向一维分布，该模型被用来评估涡旋加速的推力 [7] ；后来Mikellides 在研究中指出粘性作用是限制等离子体涡旋速度增加的重要因素，假定涡旋流动在通道出口处得到了充分的发展，并假定电流密度沿电极均匀分布，针对涡旋速度建立了一个忽略惯性力、考虑洛伦兹力和粘性力平衡的常微分方程系统，得到通道出口处涡旋速度的解析一维分布 [8] ，该一维解充分把握了粘性作用对涡旋速度径向分布特征的影响。数值模型方面，Kimura 从广义欧姆定律出发，忽略了流动在径向上的分布，数值计算了电流密度分布，分析了仅在洛伦兹力作用下涡旋速度沿轴向的一维加速过程 [9] ；Mikellides 利用 MACH2

　软件，从包含粘性的磁流体力学方程出发，数值计算出了在特定工况（电流1kA，功率50kW，附加磁场 0.1T）下涡旋速度的二维分布。

　总结来看，基于此磁流体力学的涡旋流动的数值研究已发展到二维，取得了丰硕的成果，但受限于计算量庞大，难以高效地开展变工况特性与机理研究。而求解磁流体动力学（MHD）偏微分方程的二维解析解往往面临巨大的困难，以至于先前的研究停留在了零维代数模型和一维常微分方程模型，不能全面考察惯性力、粘性力、洛伦兹力综合作用下的涡旋速度的轴向、径向二维分布问题。求解控制该涡旋流动的二维 MHD 方程组的困难在于，洛伦兹力的存在使得涡旋速度的动量方程是非齐次的，非齐次项通过广义欧姆定律与电场耦合在一起，且由于涡旋动生电动势的存在，这种耦合呈现出非线性特性。尽管困难巨大，但是二维解析解的建立对于理解该涡旋流动有着重要意义。因为在 AF-MPDT 放电通道内，涡旋速度的轴向分布问题和径向分布问题从本质上代表了两个不同的物理过程。轴向维度上看，从入口到出口，装置通过电流放电和附加磁场，将电能转化为角涡旋动能，涡旋速度会呈现逐渐增加的分布特性；径向维度上看，从通道中部到两侧电极，受粘性作用涡旋速度逐渐降低为零，涡旋动能逐渐耗散为粘性热。而且这两个物理过程不是相互独立的。涡旋速度在径向上的耗散特性直接决定了动生电动势在压降组分中的占比，而动生电动势通过电势平衡影响电流密度的轴向分布，进而影响洛伦兹力的轴向分布，最终耦合入涡旋沿轴向的启动特性中。而这种电场-流场在二维空间上的强耦合性质是零维或一维模型无法洞察的，因此获取二维的解析变的至关重要，这也正是本文研究的目标。

　图 1 涡旋流动的二维物理过程 本文针对同轴圆柱电极的 AF-MPDT 通道内的磁流体涡旋流动问题，基于偏微分方程的Sturm-Liouville 理论，首先利用分离变量法得到涡旋速度分布关于电流密度分布的无穷级数的泛函显式表达。再将该无穷级数带入广义欧姆定律中，取首次积分的级数主项，实现流场与电场的解耦，最终获得涡旋速度和电流密度的二维解析解。该解析解为 Fourier-Bessel 空间上的无穷级数形式。本文分析了该级数解的收敛特性、与二维数值仿真的结果及一维解析模型的结果的对比、所呈现出的流动特征及其背后的物理机制。

　2. 方法 图 1 展示了 AF-MPDT 电极流动通道的结构示意。如图所示，柱坐标系原点布置在通道入口中心，阴极半径为 R c ，阳极半径为 R a ，电极通道长度为 L，流经装置的总电流为 I，在通道内部，E 为电场强度，B 为附加磁场强度，J 为电流密度。采用低磁雷诺数假设 [10] ，在该坐标系中仅考虑轴向磁场 B z 和径向的电场 E r ，等离子体运动仅轴向的 u z 和角向的 u θ [11] :

　 图 2 AF-MPDT 流动通道概念图 所有物理量关于角向对称：

　0

　(1) 考虑稳态粘性磁流体动量方程，其中粘性系数 μ：

　p         u u J B u

　 (2) 结合基本假设得到其角向分量：

　2 22 2 21z r zu u u u u u J Br r r z z r               

　(3) 这是一个关于 u  的线性非齐次偏微分方程。关于边界条件的分析：在通道入口处，角向速度应为零；根据电极表面无滑移条件，壁面上的角向速度也为零；假定出口面无剪切。

　由洛伦兹力引入的非齐次项是空间坐标的函数，与附加磁场分布和电流密度分布有关。附加磁场分布 B z (r)=B 0 b(r)，其中 B 0 表征磁场强度大小，b 是无量纲的分布函数，均由附加磁场给定。再通过电荷守恒方程：

　0   J

　 (4) 直接积分得到待定解：

　( )2rI zJLr 

　 (5) 考虑广义欧姆定律：

　1e       J E u B J B

　 (6) 定义霍尔参数 

　=  B 0 /  e ，其中  为电导率，  e 为电荷密度。解得其径向分量：

　2 2r zr rJ Jb E u B   

　 (7) 对其从阴极到阳极径向积分可得电阻压降、霍尔压降、涡旋动生电动势压降：

　constantTotal Hall emfV V V V   

　 (8) 做无量纲处理：

　02/2 /2ca a a aθaR z r Lx y l kR R R Rπlμu V Jδ j wI πσL I πLR IB     

　(9)

　引 入 Reynolds 数 /z aRe u R    ， 表 征 惯 性 力 和 洛 伦 兹 力 之 比 。

　Hartman 数0/aM R B    ，表征洛伦兹力与粘性力之比。整理(3)和(8)可得到：

　2 22 2 210( , )( )( , ) ( )100y k yxx lx ybw w w w wRe iy y y y x xiwyx y xywwwx             

　(10)

　2021 1( )( ) ( ( )ln(1/ ) )())d, ; (Tklkb yy w x x k x dy b dyyx x ly      

　(11) 由于方程为(10)Sturm-Liouville 型问题，利用分离变量法寻找解 w = f (x) g (y)，先通过齐次方程和非齐次的边界条件寻找特征函数系：

　221 1" 0( ) (1) 0g"" g gy yg k g          

　(12) 解得特征函数系：

　1 1 1 1( ) ( ) ( ) ( ) ( )n n n n ng y J y Y J Y y      

　 (13) 其中 J 1 、Y 1 分别为一阶贝塞尔函数。其中 λ n 为正特征值，满足 g(λ n )=0。固有值的特征函数g n 具有正交性和完备性，为满足带有非齐次项的方程，考虑非齐次项在贝塞尔空间上展开，利用 Fourier-Bessel 展开式的正交性质求得系数：

　1121( , ) ( ) d( , ) ( )( ) dnknnnki x y g y y yi x y g yg y y y

　(14) 则得到关于 f 的方程：

　1212( ) ( ) d( ) "( ) ( )( ) d(0) 0"(( )) 0nknnky g y y yf"" x Re f x fbxg y y yxff l      

　 (15) 设其通解为 f n (x)，最终求得方程的解析解：

　1( , ; ) ( ) ( )n nnw x y f x g y  

　 (16)

　它是电流分布函数 

　(x)的泛函，为了求得方程(11)中对于动生电动势的压降积分系数，只取级数的主项：

　 1 11 11( ) ( , ; ) ( ) ( ) ( ) ( )n nk kny w x y f x y b dy b G d f x y y     

　(17) 第三章结果里的后验分析表明这种做法的精度是绝对足够的。然后将表达式(17)带入方程(11)，得到关于 f 1 的自洽方程，最终求解得到确定的 f n ：

　  2 21 122 211( ) 1 / ( )1; ,/1( )1n nnnnnnx xKKKf x ReK           

　(18) 其中各种变量错综复杂但又无法避免，现整理如下：

　 21,1 ,1 ,2 ,2 ,2 ,1,1 ,1 ,2 ,1122 21212 212 221 10( ) ( )( )exp( )exp( ) exp( )exp( )exp( )ln(1/ )( )exp( )4 124211111 ( )dn n n n n nn nkn nknnn n n nnb dydy g yg y yMKx xRe ReRe Reykl lxl lKnnx xl                    22 1011 ( ) dllx xl    最终的涡旋速度解：

　1( , ) ( ; , ) ( )n nnw x y f x Re g y K 

　 (19) 3. 结果 3.1 解的收敛特性 涡旋速度的解析表达(19)为无穷级数的形式，该级数的收敛速度对实际能否高效计算涡旋速度至关重要，同时也是近似(17)是否合理的后验证明。为了比较不同阶级数量级的相对大小，我们用级数的首项作为基准计算各级的相对大小：

　1 1( ) ( )( , )( ) ( )n nMax MaxMax Maxg y f xMAX k ng y f x

　(20) 令人惊讶的是，数学分析表明，级数的收敛特性仅仅与阴阳极半径比 k 有关，而与电流、磁

　场、密度等物理参数均无关！这是与数值方法迥然不同的一个结论，也是该解析解一个非常长的优点。这意味着在实际计算中，该解析解对于不同的工况条件具有强适应性。图 3 展示了在不同内外径之比 k 下，级数的各阶（n 表示阶数）的相对大小。从该图我们可以看出，级数收敛十分迅速，且 k 越大，级数收敛越快；对于 k=0.1-0.9，取级数首项的精度是绝对足够的，这证明了近似处理(17)是合理且可靠的。在实际应用该解析解计算涡旋速度时，k=0.1或 0.2 取级数前三项即可达到 95 %的精度，k>0.3 时仅取首项能达到 95 %的计算精度。

　图 3 级数解的收敛特性 总结来看，仅与内外径比值 k 有关且十分迅速的收敛特性，使整个解析解成为一个计算涡旋速度高效实用的工具。由于级数解是显式数学表达，且不存在迭代求解的过程，实际计算的时间完全可以忽略不计，相比数值方法会节省大量的时间成本。

　3.2 解的对比验证 解析解(19)的可靠性是一个值得关注的问题：简化建模与冗长的数学处理是否把握了涡旋流动的物理本质特征，能否展示出足够多的流动细节？由于 NASA Lewis 的 50kW 级AF-MPDT 具有十分丰富的实验数据 [12] 和详细的数值结果 [3] ，本文以该推力器（总电流 1000A，通道内电流840A [3] ，附加磁场0.1T，入口流量0.1g/m 3 ，通道外径0.051m，内外径比值k=0.25）为具体对象展示该解析解与先前数值模型和一维模型的对比结果，如图 4 所示。

　图 4 采用不同模型计算的出口端面涡旋速度的径向分布对比

　为了便于对比，图 4 呈现的是采用不同模型计算的出口端面涡旋速度的径向分布，蓝点为 Mikellides 采用 MACH2 软件得到的数值仿真结果 [3] ，绿线为基于常微分方程的一维模型计算结果 [8] ，红线为本文基于偏微分方程得到的二维解析解。对比来看，一维模型和二维模型和数值仿真均呈现出了粘性作用对涡旋速度的限制，导致速度型呈现出中间高，靠近电

　极逐渐降低为零的特性。但一维模型由于忽略了轴向的洛伦兹力和粘性力的作用，相比于二维模型和数值仿真计算出的涡旋速度偏大。二维解析解与数值仿真结果整体上吻合较好。差异在于右半边区域里（阳极区），数值解呈现出微弱的色散性质。这种色散性是由于数值算法导致还是真实的物理特性还有待详细考察。但值得注意的是，更进一步的研究表明这种色散性在高雷诺数时同样出现在了二维的解析解中，但鉴于不是本文重点，此处不予展开阐述。

　总结来看，对比验证表明该解析解与数值结果较好地吻合，具有相当的可信度。

　3.3 解的物理特征 仍然在 3.2 节描述的 NASA Lewis 的 AF-MPDT 的物理参数下，图 5 和图 6 分别展示了涡旋速度的二维分布结果和径向电流密度的分布结果。

　从径向上看：

　（1）涡旋速度在通道中间最高，靠近两侧电极逐渐降低为 0，如图 5 所示；这是由于粘性效应和壁面的无滑移边界条件导致。（2）电流密度越靠近阴极越大，如图 6所示；这是由于电荷守恒和同轴圆柱的几何构型导致。（3）靠近阴极更大的电流密度导致更大的洛伦兹力，克服粘性作用的能力更强。进而使得涡旋速度的最大值更靠近阴极，如图 5所示。从轴向上看：（1）受洛伦兹力加速效果，涡旋速度从入口到出口逐渐增加，如图 5所示。（2）涡旋速度的增加导致涡旋动生反电动势的增加，但由于电极等电势的边界条件，电阻压降不得不降低，从而电流密度从入口到出口逐渐减小，如图 6 所示。（3）电流密度沿轴向的递减导致洛伦兹力从入口到出口逐渐减小，从而涡旋加速的梯度逐渐降低，在出口处涡旋速度等高线近似水平，如图 5 所示。该解析解呈现的物理特征鲜明合理。

　图 5 角向涡旋速度二维分布

　图 6 径向电流密度二维分布

　4. 结 结 论 与展望 本文基于 MHD 方程得到了 AF-MPDT 流动通道内涡旋流动的二维解析解。该二维解析解被表达为 Fourier-Bessel 空间上的无穷级数形式。数学分析表明级数收敛速度仅与阴阳极半径之比有关，且收敛十分迅速，取级数首项即可满足绝大多数情况下的精度需求，计算时间成本可忽略不计。相比于一维解析结果，该二维解析解与二维数值仿真结果具有更好的吻合程度。该解析解充分反映了涡旋速度和电流密度的分布特征的自洽，这背后是对流场-电场闭环耦合特征的完整把握。

　该二维解析模型为 AF-MPDT 涡旋流动的计算分析提供了一个强有力的工具。在此基础上可继续开展对通道内电磁流动的能量转化特性、附加磁场分布的影响、流场电场非线性耦合的详细研究。同时也从理论上为普遍意义下环形通道正交电磁场作用下磁流体涡旋流动提供了解析方法。

　参考文献 [1] Tang H-B, Cheng J, Liu C, 等. Study of applied magnetic field magnetoplasmadynamic thrusters with particle-in-cell and Monte Carlo collision. II. Investigation of acceleration mechanisms[J]. Physics of Plasmas, 2012, 19(7): 073108. [2] Kodys A, Choueiri E. A Critical Review of the State-of-the-Art in the Performance of Applied-Field Magnetoplasmadynamic Thrusters[A]. 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit[C]. Tucson, Arizona: American Institute of Aeronautics and Astronautics, 2005. [3] Mikellides P G, Turchi P J, Roderick N F. Applied-Field Magnetoplasmadynamic Thrusters, Part 1: Numerical Simulations Using the MACH2 Code[J]. Journal of Propulsion and Power, 2000, 16(5): 887–893. [4] Blackstock A W, Fradkin D B, Liewer K W, 等. Experiments using a 25-kw hollow cathode lithium vapor MPD arcjet[J]. AIAA Journal, 1970, 8(5): 886–894. [5] Albertoni R, Paganucci F, Andrenucci M. A phenomenological performance model for applied-field MPD thrusters[J]. Acta Astronautica, 2015, 107: 177–186. [6] Tobari H, Ando A, Inutake M, 等. Characteristics of electromagnetically accelerated plasma flow in an externally applied magnetic field[J]. Physics of Plasmas, 2007, 14(9): 093507. [7] Sasoh A, Arakawa Y. Thrust Formula for Applied-Field Magnetoplasmadynamic Thrusters Derived from Energy Conservation Equation[J]. Journal of Propulsion and Power, 1995, 11(2): 351–356. [8] Mikellides P G, Turchi P J. Applied-Field Magnetoplasmadynamic Thrusters, Part 2: Analytic Expressions for Thrust and Voltage[J]. Journal of Propulsion and Power, 2000, 16(5): 894–901. [9] Tanaka M, Kimura T. Current distribution and plasma acceleration in MPD arcjets with applied magnetic fields[J]. Journal of Propulsion and Power, 1988, 4(5): 428–436. [10] Suslov S A, Pérez-Barrera J, Cuevas S. Electromagnetically driven flow of electrolyte in a thin annular layer: axisymmetric solutions[J]. Journal of Fluid Mechanics, 2017, 828: 573–600. [11] Baylis J A, Hunt J C R. MHD flow in an annular channel; theory and experiment[J]. Journal of Fluid Mechanics, 1971, 48(3): 423–428. [12] Myers R, Mantenieks M, Sovey J. Geometric effects in applied-field MPD thrusters[A]. 21st International Electric Propulsion Conference[C]. Orlando,FL,U.S.A.: American Institute of Aeronautics and Astronautics, 1990.

　联系作者：

　张凯宇 zkyae@buaa.edu.cn 17801004122