航天器在轨运行过程中需要对自身姿态进行调整和稳定控制，常用的方式是利用动量交换装置与航天器本体实现角动量交换^[1]，动量轮是该工作中常用的设备. 传统的航天器位姿控制主要通过控制力矩陀螺或惯性动量轮机构实现，其中控制力矩陀螺一般用于输出较大的力矩和角动量，适合大型航天器的姿态控制，惯性动量轮则更适用于对姿态控制精度的调整^[2]. 传统的动量轮以机械式为主，滚珠轴承磨损、机械振动和冲击等因素会显著影响其控制精度和使用寿命^[3-4]. 基于磁流体动力学（magnetohydrodynamics，MHD）原理的新型动量轮设备通过电磁场耦合驱动金属流体在封闭环管内运动，可以输出相应的角动量和力矩. 流体可以弥补机械组件的摩擦损耗、抑制扰动和延长使用寿命，因此在微小卫星的位姿控制问题中适用性良好.

1958年，Haviland^[5]设计用于飞行器方向控制的磁流体致动器，通过电磁泵驱动三坐标轴面内的环形导电流体管道，输出力矩控制飞行器的位姿. 类似的研究成果还包括Davis^[6-7]在1968年发表的利用太阳能电源的磁流体位姿方向控制器以及1969年的三坐标轴磁流体角度稳定装置. 1988年，NASA的Maynard^[8]发明流体动量控制器，提出通过在管道环路内泵送流体以输出较大动量的设计方案. 2007年，美国ATA公司的Laughlin^[9]提出MHD致动/传感复合装置的设计原理和输出指标计算方法，初步得到性能标度因数，但关于流体的致动原理分析不够完善. 2013年，Noack等^[10]设计基于洛伦兹力的新型MHD致动器，使用电磁泵加速环形圆管内的金属流体获得力矩输出，并在气浮轴承平台上进行测试实验，与原有致动器对比，简要说明输出性能的改进效果. 2015年，法国空间研究中心与图卢兹国立综合理工学院的Casteras等^[11-12]也对MHD致动器开展相关的研究，建立电矢量模型和二维有限差分模型并进行仿真分析，对概念样机进行初步测试. 2017年，Curti^[13]研究MHD致动器的结构布局，搭建有限差分混合模型和集总参数模型，针对关键参数进行仿真分析，探究输出性能的几项影响因素.

近年来，国内的相关研究主要集中于机械式动量轮在卫星位姿控制中的应用，对磁流体致动器的研究起步晚，相关成果较少. 2017年，王志远^[14]对皮纳卫星姿态控制所需的微型机械动量轮组设计方法及控制算法进行相关研究. 同年，吴启东^[15]针对微小卫星机械式动量轮的机电能量转换问题进行有限元分析和驱动方案设计. 2018年，程旭^[16]通过COMSOL有限元仿真软件对基于电磁效应的电磁流体环进行初步仿真分析，对原理样机进行粗略测试. 总体而言，国内关于磁流体动力学动量轮的原理研究和致动特性分析仍须进一步深入.

本文结合磁流体动力学基本原理，对电磁场作用下矩形管道金属流体的流动规律和哈脱曼边界层对流动特性的影响原理进行探究，对电流和电压控制模式下动量轮的频域输出特性进行对比分析，通过有限元仿真软件COMSOL实现整机多物理场建模和耦合仿真，并对理论和仿真结果进行对照.

MHD动量轮的主要原理是通过施加相互正交的电场与磁场产生安培力，驱动金属流体在管道内流动，流体产生的角动量和转矩可传递到负载，用于控制航天器的位姿状态. 以矩形截面的环管模型为例，施加轴向磁场、径向电场时的基本模型如图1所示. 由图可知，磁场和电场分别由底部的永磁铁和内外电极施加，通过控制电极内的电流或流体内外壁电压，可以实现对流体输出角动量和力矩的控制. 动量轮中流体环受力为电磁场施加的安培力与流动产生的黏滞力的合效果，根据力的相互作用原理，负载端会得到相应的力矩输出.

MHD动量轮涉及的电磁场方程包括：

(1) $ \nabla \times {\boldsymbol{E}}{\text{ = }} - \frac{{\partial {\boldsymbol{B}}}}{{\partial t}}, $

(2) $ \nabla \times {\boldsymbol{H}} = {\boldsymbol{J}}, $

(3) $ \nabla \cdot {\boldsymbol{B}} = 0, $

(4) $ \frac{{\partial {\boldsymbol{B}}}}{{\partial t}} = \nabla \times ({\boldsymbol{v}} \times {\boldsymbol{B}}) + \eta {\nabla ^2}{\boldsymbol{B}}, $

(5) $ {\boldsymbol{J}} = \sigma ({\boldsymbol{E}} + {\boldsymbol{v}} \times {\boldsymbol{B}}). $

式中：v为流体流速，η为磁扩散系数，σ为电导率，E为电场强度，B为磁感应强度，H为磁场强度，J为电流密度.

涉及的流体力学方程为流体的连续性方程和纳维−斯托克斯方程：

(6) $ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho {\boldsymbol{v}}) = 0, $

(7) $ \rho \left(\frac{{\partial {\boldsymbol{v}}}}{{\partial t}} + ({\boldsymbol{v}} \cdot \nabla ){\boldsymbol{v}}\right) = - \nabla {\boldsymbol{p}} + \mu {\nabla ^2}{\boldsymbol{v}} + {\boldsymbol{f}}. $

式中：ρ为流体密度， $\nabla $p为流体压强梯度，μ为流体的动力黏度，f为流体附加的体积力. 在电磁场驱动下，附加体积力项为 $ {\boldsymbol{f}}{\text{ = }}{\boldsymbol{J}} \times {\boldsymbol{B}} $，导电流体内的电场可以通过控制输入电流或内外壁间电压施加.

为了简化模型，流体须进行以下假设：1）流体的磁雷诺数Re_m远小于1，式（4）对应的磁感应方程退化为扩散方程的形式，流体流动产生的影响项 $ \nabla \times ({\boldsymbol{v}} \times {\boldsymbol{B}}) $可以忽略，金属流体表现为磁扩散效应，式（4）简化为

(8) $ \frac{{\partial {\boldsymbol{B}}}}{{\partial t}} = \eta {\nabla ^2}{\boldsymbol{B}}. $

2）导电流体设定为不可压缩流体，流体密度和黏度均为常值，在该假设条件下流体连续方程（6）可改写为

(9) $ \frac{{\partial \rho }}{{\partial t}} = 0, $

(10) $ \nabla \cdot {\boldsymbol{v}} = 0. $

3）流体运动状态按照层流规律进行设定计算，雷诺数初步设定满足 ${{Re}} \leqslant 2\;300$. 4）导电流体充满管道，并与壁面充分接触，管道内壁设置为无滑移条件，近壁面流体与壁面的相对速度为0. 5）流体的重力项和其他物理场（如温度场）对流体特性的影响暂不考虑. 流体运动过程中产生的角动量L和转矩T计算公式分别为

(11) $ {\boldsymbol{L}} = {I_Z}{\boldsymbol{\omega }}, $

(12) $ {\boldsymbol{T}} = {I_Z}{\boldsymbol{\beta }}. $

式中： ${I_Z} = 0.5{{\rho {\text{π}} }}h(R_2^4 - R_1^4)$，为环筒结构绕转轴的转动惯量，h为流体环高度，R₁与R₂为环的内外半径；ω、β分别为流体的转动角速度、角加速度. 由于动量轮的壳体与微小卫星等载体固连，壳体与负载一同作为流体环的转矩输出对象，流体环和负载所构成的系统合外力矩为0，满足角动量守恒条件. 因为流体部分转动惯量远小于负载的转动惯量，所以外壳的角速度要远小于流体角速度，因此二者的相对速度近似等于流体的流速.

将式（7）中的纳维-斯托克斯方程在柱坐标系 $ (r,\theta ,z) $下进行分量展开，得到沿径向、周向和轴向的分量方程分别为

(13) $ \left. {\begin{array}{*{20}{l}} \begin{array}{l} \dfrac{{\partial {v_r}}}{{\partial t}} + {v_r}\dfrac{{\partial {v_r}}}{{\partial r}} - \dfrac{{v_\theta ^2}}{r} + \dfrac{{{v_\theta }}}{r}\dfrac{{\partial {v_r}}}{{\partial \theta }} + {v_z}\dfrac{{\partial {v_r}}}{{\partial z}} = \dfrac{{{f_r}}}{\rho } - \dfrac{1}{\rho }\dfrac{{\partial p_r}}{{\partial r}} + \\ \;\;\;\;\;\;\;\;\dfrac{\mu }{\rho }\left(\dfrac{{{\partial ^2}{v_r}}}{{\partial {r^2}}} + \dfrac{1}{r}\dfrac{{\partial {v_r}}}{{\partial r}} - \dfrac{{{v_r}}}{{{r^2}}} - \dfrac{2}{{{r^2}}}\dfrac{{\partial {v_\theta }}}{{\partial \theta }} + \dfrac{1}{{{r^2}}}\dfrac{{{\partial ^2}{v_r}}}{{\partial {\theta ^2}}} + \dfrac{{{\partial ^2}{v_r}}}{{\partial {z^2}}}\right), \end{array}\\ \begin{array}{l} \dfrac{{\partial {v_\theta }}}{{\partial t}} + {v_r}\dfrac{{\partial {v_\theta }}}{{\partial r}}{\rm{ + }}\dfrac{{{v_r}{v_\theta }}}{r} + \dfrac{{{v_\theta }}}{r}\dfrac{{\partial {v_\theta }}}{{\partial \theta }} + {v_z}\dfrac{{\partial {v_\theta }}}{{\partial z}} = \dfrac{{{f_\theta }}}{\rho } - \dfrac{1}{{\rho r}}\dfrac{{\partial p_\theta}}{{\partial \theta }} + \\ \;\;\;\;\;\;\;\;\dfrac{\mu }{\rho }\left(\dfrac{{{\partial ^2}{v_\theta }}}{{\partial {r^2}}} + \dfrac{1}{r}\dfrac{{\partial {v_\theta }}}{{\partial r}} - \dfrac{{{v_\theta }}}{{{r^2}}} + \dfrac{2}{{{r^2}}}\dfrac{{\partial {v_r}}}{{\partial \theta }} + \dfrac{1}{{{r^2}}}\dfrac{{{\partial ^2}{v_\theta }}}{{\partial {\theta ^2}}} + \dfrac{{{\partial ^2}{v_\theta }}}{{\partial {z^2}}}\right), \end{array}\\ \begin{array}{l} \dfrac{{\partial {v_z}}}{{\partial t}} + {v_r}\dfrac{{\partial {v_z}}}{{\partial r}} + \dfrac{{{v_\theta }}}{r}\dfrac{{\partial {v_z}}}{{\partial \theta }} + {v_z}\dfrac{{\partial {v_z}}}{{\partial z}} = \dfrac{{{f_z}}}{\rho } - \dfrac{1}{\rho }\dfrac{{\partial p_z}}{{\partial z}} + \\ \;\;\;\;\;\;\;\;\dfrac{\mu }{\rho }\left(\dfrac{{{\partial ^2}{v_z}}}{{\partial {r^2}}} + \dfrac{1}{r}\dfrac{{\partial {v_z}}}{{\partial r}} + \dfrac{1}{{{r^2}}}\dfrac{{{\partial ^2}{v_z}}}{{\partial {\theta ^2}}} + \dfrac{{{\partial ^2}{v_z}}}{{\partial {z^2}}}\right). \end{array} \end{array}} \right\}$

在环形管道中，周向流动为主要的流动方向，径向和轴向速度分量远小于周向速度. 实际上，二次流现象是引入径向和轴向速度分量的主要原因，当矩形管道满足 $D{e^2}D/\left( {H{a^4}R} \right) << 1$条件时，可认为二次流的影响可以忽略^[17]，其中Ha为哈脱曼数， $Ha = {\left( {\sigma /\mu } \right)^{0.5}}B_0D/2 $， $ D $为管道的特征长度， $R$为弯管的曲率半径，B₀为平均磁感应强度， $ {{De}} $为迪恩数， ${{De}} =(0.5{{D}/{{R}}})^{0.5} {{\rho v_0}}D/{\mu }$，v₀为平均流速. 在上述条件下，速度和压强分量关于回转轴Z轴对称，则满足 $ {{\partial v_\theta/\partial \theta }} = 0 $且 $ {{\partial p_\theta/\partial \theta }}{\text{ = }}0 $，周向分式可初步简化为

(14) $ \frac{{\partial {v_\theta }}}{{\partial t}} = \frac{{{f_\theta }}}{\rho } + \frac{\mu }{\rho }\left(\frac{{{\partial ^2}{v_\theta }}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {v_\theta }}}{{\partial r}} - \frac{{{v_\theta }}}{{{r^2}}} + \frac{{{\partial ^2}{v_\theta }}}{{\partial {z^2}}}\right). $

在相互垂直的电磁场作用下，不可压缩黏性导电流体沿均匀矩形截面管道的定常层流运动称为哈脱曼流动，该理论由丹麦物理学家朱利叶斯·哈脱曼于1937年通过实验证实^[18]. 以二维平板间的哈脱曼流动为例，流体在流动方向y上达到稳定时满足

(15) $ \mu \frac{{{d^2}{\boldsymbol{v}}}}{{d{z^2}}} - \frac{{\partial {\boldsymbol{p}}}}{{\partial y}} + \sigma ({\boldsymbol{E}} - {\boldsymbol{v}} \times {\boldsymbol{B}}) = {\boldsymbol{0}}. $

在z方向求得哈脱曼流动速度v_y的分布解^[19]为

(16) ${v_y}\left( z \right) = \frac{{{v_0}Ha\left[ {\cosh \;(Ha) - \cosh \;\left( {\dfrac{{zHa}}{d}} \right)} \right]}}{{Ha\cosh \;(Ha) - {{\sinh }}\;(Ha)}}.$

式中： $ d $为板间距的一半. 流速分布如图2所示. 由图可知，当 $ {{Ha = 0}} $时，截面流速分布规律近似为抛物线分布，与普通管道泊肃叶流动规律类似；随着Ha增大，中央位置的流速分布将逐渐变均匀，靠近上下绝缘壁面处的薄层内存在较大的速度梯度和黏滞切应力^[20]，称为哈脱曼边界层. 在高Ha条件下，该边界层厚度 ${\delta _1}$与 $d{{H}}{{{a}}^{ - 1}}$为同一数量级^[21-22]. 根据前文假设，环形管道中的周向流速v_θ近似等于总体平均流速v₀. Baylis等^[21]指出，哈脱曼边界层中的黏滞力项满足等数量级替换关系：

(17) $ \left| {\mu \frac{{{\partial ^2}{v_\theta }}}{{\partial {z^2}}}} \right|\sim \Theta \left[\mu {v_0}\frac{{{{H}}{{{a}}^2}}}{{{d^2}}}\right]. $

根据哈脱曼流动的速度分布曲线可知，位于管道中央的高流速区域的法向流速梯度非常小，因此流体的黏滞力高度集中于哈脱曼边界层^[23]，因为洛伦兹力项 $ {\boldsymbol{f}}{\text{ = }}{\boldsymbol{J}} \times {\boldsymbol{B}} $作为驱动力以克服黏滞阻力，哈脱曼边界层中电流密度J高度集中，相比于高流速区域而言，提供了电流的低阻抗通路^[24]，高流速区域因为式（5）中 $ {\boldsymbol{v}} \times {\boldsymbol{B}} $项引入感应电场抵消原电场效果，所以实际电流密度较低.

对于三维矩形管道模型，导电壁面附近也存在类似的侧壁边界层，厚度 ${\delta _2}$与 $ d{{H}}{{{a}}^{ - {0.5}}} $数量级相同，黏滞力项等数量级替换为

(18) $ \left| {\mu \frac{{{\partial ^2}{v_\theta }}}{{\partial {z^2}}}} \right|\sim \Theta \left[\mu {v_0}\frac{{{{Ha}}}}{{{d^2}}}\right]. $

该层内的黏滞力影响明显小于哈脱曼边界层.

流体环截面示意图如图3所示，根据任务书中的尺寸指标和动量轮小型化设计要求，设流体环高度h=3 mm，矩形截面宽度a等于环内外半径差，内外径分别为R₁=12 mm、R₂=38 mm，流体环的特征长度为 $ D = 2ah/\left( {a + h} \right) = 5.38\;{\rm{ mm}}$，内外壁面为电极面，上下壁面为绝缘面. 流体材料选择镓铟锡（GaInSn），该金属材料熔点仅为6~12 ℃，常温下（20 ℃）电导率σ=3.27×10⁶ S/m，动力黏度μ=2.159×10⁻³ Pa·s. 轴向磁场由钕铁硼（NdFeB）永磁铁提供，该磁铁是目前仅次于绝对零度钬磁铁的剩磁强度最高的实用永磁材料.

在哈脱曼边界层中，大部分位置到上下壁面的距离要远远小于到内外壁面的距离，因此满足关系式 ${\partial ^2}{v_\theta }/\partial {z^2} > > {\partial ^2}{v_\theta }/\partial {r^2}$，根据式（17）中的代换可知径向梯度项 $ {r^{ - 1}}\left( {\partial {v_\theta }/\partial r} \right) $和 $ - {{{v_\theta/r^2}}} $相比于 ${\partial ^2}{v_\theta }/\partial {z^2}$项也可以省略，式（14）可以进一步简化为

(19) $ \rho \frac{{\partial {v_\theta }}}{{\partial t}} = {f_\theta } + \mu \frac{{{\partial ^2}{v_\theta }}}{{\partial {z^2}}}. $

对于电流控制模式，由于绝大部分电流由哈特曼边界层通过， $ {f_\theta } $中的电流密度近似等于电流除以上下两层哈脱曼边界层对应的平均柱面面积，得到哈脱曼边界层内平均电流密度约为

(20) $\bar J = \frac{I}{{4{\text{π}} \bar R{\delta _1}}}.$

式中： $\overline R $为流体环的平均半径，I为电流. 对于哈脱曼边界层内的流体，当流速达到稳定时，层内的洛伦兹力与黏滞力平衡，联立式（17）、（19）、（20）并代入哈脱曼边界层厚度δ₁得到：

(21) $\frac{{{B_0}IHa}}{{2{\text{π}} \bar RD}}{\rm{ }} = {\rm{ }}\frac{{\mu {v_0}H{a^2}}}{{{{(D/2)}^2}}}. $

得到稳态流速解为

(22) $ {v_0} = \frac{I}{{4{\text{π}}\overline R \sqrt {\mu \sigma } }}. $

进一步对流体环整体进行受力分析，总的安培力大小为

(23) $ {F_{\rm{e}}} = B_0I({R_2} - {R_1}). $

上下壁面处哈脱曼边界层中的黏滞阻力为

(24) $ {F_{{\rm{s}}1}} = 2\mu {v_0}\frac{{{{H}}{{{a}}^2}}}{{{{(D/2)}^2}}}{\delta _1}{\text{π}}(R_2^2 - R_1^2). $

内外侧壁处边界层中的黏滞阻力为

(25) ${F_{{\rm{s}}2}} = \mu {v_0}\frac{{{{Ha}}}}{{{{(D/2)}^2}}}{\delta _2}2{\text{π}}({R_1} + {R_2})h.$

流体环整体在任意时刻的受力关系为

(26) $ m\frac{{\partial {v_0}}}{{\partial t}} = {F_{\rm{e}}} - ({F_{{\rm{s}}1}} + {F_{{\rm{s}}2}}). $

式中：m为流体质量， $m{{ = {\text{π}} }}\rho h(R_2^2 - R_1^2)$，用角速度替换线速度，联立式（23）~（26），并对式（26）进行拉普拉斯变换，得到电流控制模式下的传递函数

(27) $ \dfrac{{\omega (s)}}{{I(s)}} = \dfrac{{\dfrac{{2{B_0}}}{{\rho h{\text{π}} {{({R_1} + {R_2})}^2}}}}}{{s + \left[ {\dfrac{{4\mu Ha}}{{\rho Dh}} + \dfrac{{4\mu (Ha) ^{0.5}}}{{\rho D({R_2} - {R_1})}}} \right]}}. $

动量轮的频域特性为低通一阶惯性环节，其转折角频率

(28) $ {\varOmega _{\text{0}}}{\text{ = }}\frac{{4\mu {{Ha}}}}{{\rho Dh}} + \frac{{4\mu (Ha) ^{0.5} }}{{\rho D({R_2} - {R_1})}}. $

求得 $ {{\varOmega }_{\text{0}}} = 5.36{\text{ rad/s}} $，对应的特征时间参数 $ \tau {\text{ = }}1/{{\varOmega }_{0}} \approx $ $ 0.187{\text{ s}} $. 采用电压控制模式将式（5）沿半径方向积分，得到壁面间电压

(29) $ U = B_0{v_0}({R_2} - {R_1}) + I{r_z}, $

(30) $ {r_z} = \frac{{\ln \;({R_2}/{R_1})}}{{2{\text{π}} \sigma h}}.$

式中： $ {r_z} $为流体环径向电阻. 联立式（26）、（29）、（30），整理得到电压控制模式下的传递函数

(31) $\dfrac{{\omega (s)}}{{U(s)}} = \dfrac{{\dfrac{{2B_0}}{{\rho {\text{π}}hr{{({R_1} + {R_2})}^2}}}}}{{s + \left[\dfrac{{4\mu {{Ha}}}}{{\rho Dh}} + \dfrac{{4\mu (Ha) ^{0.5} }}{{\rho D({R_2} - {R_1})}} + \dfrac{{{B_0^2}({R_2} - {R_1})}}{{\rho {\text{π}}hr({R_1} + {R_2})}}\right]}}. $

对于电压控制模式，其一阶惯性环节形式与电流模式一致，但转折角频率增大至178 rad/s，特征时间缩短至0.005 6 s，因此电压控制模式的响应速度更快. 但导电流体电阻率很低，实际加在流体环内外壁间的电压微小、仅为毫伏量级，在导线和电极电阻影响下，该电压值难以直接准确控制，因此在实际工作中，电压控制模式难以实现，电流控制方案更符合实际测试要求.

为了深入分析磁流体动量轮导电流体环的运动特性，可以通过数值仿真工具对流体的运动情况进行计算并验证理论的正确性. 常用的数值仿真计算方法包括有限差分法、有限体积法、有限元法和有限分析法等，其中有限元仿真方法在多物理场耦合问题中具有很好的适用性，应用广泛且实用高效. 有限元法使用数学近似方法模拟真实物理系统，利用有限数量的离散单元近似代替复杂的实际多物理场模型^[25]，理论基础是变分原理和加权余量法. 使用基于有限元法的仿真软件COMSOL对MHD动量轮整机进行电−磁−流多物理场建模和耦合仿真，探究金属流体受驱动后的运动特性和流速分布规律，并与理论计算结果进行对照.

在COMSOL中调用磁场和电场（mef）和层流（spf）模块，控制通过电极端面的电流，磁钢通过设置剩磁强度参数后产生相应的磁场，模块之间通过洛伦兹力项 ${\boldsymbol{J}} = \sigma ({\boldsymbol{E}} + {\boldsymbol{v}} \times {\boldsymbol{B}})$实现多物理场耦合计算. 将流体流动的离散化状态设置为二阶离散化，使用二阶基函数描述速度和压力. 针对耦合时存在较强的非线性问题，采用分离式求解设置有效降低内存消耗，采用相应的优化迭代求解器进行稳态解计算.

数值仿真中模型的网格划分效果对仿真结果具有重要影响，网格细化使运算结果更接近准确值^[26]，但网格密度过大会导致计算量过大. 对完整样机几何模型进行适当简化处理，装配体截面如图4所示，涉及的主要材料参数（20 ℃）如表1所示. 表中，B_r为剩磁强度. 考虑计算精度、零件尺寸的因素，将流体区域设置为流体动力学极细化网格，电极设置为普通物理超细化网格，其余位置网格适当粗化以减少计算量，得到导电区域和磁钢结构的网格划分效果如图5所示. 图中，网格质量符合计算精度和网格收敛性要求.

在仿真测试中，选取0.25~6 A范围内的若干电流值作为输入，仿真流体的运动状态. 以6 A输入电流为例，流体环纵向和横向切面的流速分布如图6所示. 由图可知，流体在管道中部流速分布较为均匀，流速随半径增大逐渐递减，在近壁面位置有明显的高流速梯度薄层，轴向流速梯度明显高于径向梯度，金属流体运动符合哈脱曼流动规律. 进一步绘制0.25~6 A范围内流体平均流速−电流曲线图，并对照理论模型对应的曲线如图7所示. 定义流速相对偏差 $ \varepsilon = \left| {\left( {{v_2} - {v_1}} \right)/{v_1}} \right| \times 100{\text{%}} $，表征理论流速 $ {v_1} $与仿真流速 $ {v_2} $的偏差程度，相对偏差随电流的变化曲线如图8所示. 由图可知流体环的输入输出线性关系良好，仿真结果与理论推导整体吻合，相对误差均小于10%，计算得到动量轮的角动量输出标度因数约为9.68×10⁻⁵ N·m·s/A，该指标与文献[10]、[13]相比具有一定的优势. 在电流较小时相对偏差较大，原因是低流速时仿真结果偏差量的相对占比更明显，随着电流提升、流速加快，情况逐渐改善.

保持其他条件不变，输入电流设定为6 A，进一步探究磁场强弱对流体速度的影响. 目前已知的钕铁硼材料各型号中剩磁强度最高可达到1.5 T，在0.25~1.5 T设置模型磁钢剩磁参数，得到流体环区域实际磁感应强度在0.12~0.74 T. 平均流速随磁感应强度变化曲线如图9所示，相对偏差变化如图10所示. 随着磁场增强，流体流速先逐渐增大并在0.36 T附近达到最大流速，理论推导的稳态流速解与最大流速值吻合，之后继续增大磁感应强度，流速将不断减小. 进一步分析原因，当磁感应强度较小时，驱动流体的洛伦兹力较弱，黏滞阻力对流速的限制作用明显；继续增强磁场，流速先随驱动力的增大而加快，当磁感应强度增大到一定程度后，哈脱曼边界层的迅速变薄将导致近壁面黏滞力影响显著增强，同时中央高流速区域反向感应电流产生的洛伦兹力对流动的抑制作用更加明显，导致流速下降. 在上述因素影响下，磁场过强或过弱均会导致动量轮输出性能降低、理论推导吻合度下降，最优磁感应强度值与流体属性、结构尺寸参数有关.

流体平均磁感应强度保持0.36 T，对流体电导率和动力黏度参数分别进行参数化扫描，得到理论和仿真平均流速随二者变化的对照曲线分别如图11、12所示. 由图可知，磁流体平均流速仿真解与前述理论推导具有高度一致性，流速大小 $v_0$分别与 ${\sigma ^{ - 0.5}}$和 ${\mu ^{ - 0.5}}$保持正比关系，由于二者阶次相同，因此输出量对这2项参数的敏感程度基本一致. 当电导率和黏度较小时，流速较高并向无序的湍流状态过渡，导致理论计算值与层流模式下仿真结果的相对偏差逐渐增大.

综上分析，为进一步提升MHD动量轮的动量输出指标，在充分考虑材料实际工作状态和属性的前提下，可以通过增大输入电流幅值、选择剩磁强度适中的磁铁材料、选择电导率和动力黏度较低的流体材料等措施实现.

（1）针对封闭环形矩管内金属流体的哈脱曼流动问题，探究纳维−斯托克斯方程的柱坐标系分量展开、流动平衡模型简化方法，围绕流动边界层的力学特性建立电流和电压控制模式下，MHD动量轮的简化平衡方程和传递函数模型，并进行频域特性的对比分析.

（2）建立MHD动量轮的整机有限元仿真模型，通过有限元仿真软件COMSOL实现电−磁−流多物理场耦合仿真，得到流体的稳态速度分布解，仿真结果与理论解具有良好的一致性. 针对电流、磁感应强度、电导率和黏度进行参数扫描并明确各参数对流速的影响原理，对后续的动量轮性能测试和改进优化具有重要的指导作用.

（3）本文主要针对层流状态下磁流体受驱动后的运动规律进行分析，当驱动电流继续增大、流速过渡到湍流阶段后，流体中漩涡、二次流的现象增强，将加剧流动的不稳定性，流动方程和求解模型会更加复杂，文中模型的适用性有待进一步探究.