磁控热防护系统在天地往返运载器上的应用仿真

2021-08-03丁明松刘庆宗江涛董维中高铁锁傅杨奥骁

航空学报 2021年7期

丁明松,刘庆宗,江涛,董维中,高铁锁,傅杨奥骁

中国空气动力研究与发展中心计算空气动力研究所, 绵阳 621000

从1981年4月12日美国“哥伦比亚”号(OV-102)航天飞机进行首次载人轨道飞行到2011年7月8号“阿特兰蒂斯”号航天飞机第135次升空，航天飞机完成了一次又一次航天壮举[1]。尽管持续30年的航天飞机时代已经宣告结束，但人们从未放弃滑翔返回式天地往返运载器这一尖端领域，相关研究正不断推动航天科学关键技术发展[2]。

技术瓶颈是制约航天项目跨越式发展的关键。如气动热防护这类技术瓶颈问题得不到有效解决，航天型号的使用依然代价巨大。天地往返系统高速飞行时，面临的热障问题尤为突出，从未来发展趋势来看，运载器性能要求越来越高，气动热环境越来越严酷，给飞行器热防护系统设计带来极大的挑战[2]。

通过配置合适的机载电磁场系统，对高超声速流动进行控制，改善飞行器气动热环境，这就是高超声速飞行器磁控热防护技术[3]。对于高超声速流动，速度较高时，波后气体热电离产生等离子体具有足够的电导率，可直接应用于磁流体力学控制[4]；同时机载磁场通过电流响应控制，并安装于飞行器内部，无需改变飞行外部结构，因而具备响应快速、不改变气动外形、不影响其他被动热防护措施等优点[5]，满足了未来高超声速飞行器高效稳定、可重复使用的要求。

自20世纪90年代以来，得益于稀土磁工业的发展，易于制造便携的机载磁系统[4,6]，磁控热防护技术掀起了新的研究热潮。美国空军研究实验室Bisek团队基于低磁雷诺假设开展了带尾舵钝头椭圆锥体磁流体数值研究[7]，分析了外磁场对尾舵前缘热流的影响，但由于前缘局部电导率较低(<0.1 S/m)，磁控热防护效果不明显；俄罗斯高温研究所Bityurin和Bocharov[8]基于完全气体模型和偶极子磁场配置方案，开展返回舱简化外形磁控热防护数值模拟，发现通过配置磁场可有效降低驻点区域热流；日本Otsu[9]、Fujino[10]等开展了多种条件下简单钝体外形磁流体力学控制数值模拟研究，发现绝缘壁面条件霍尔效应影响几乎可以忽略，磁控热防护效果较为显著；德国宇航中心和欧洲航天局共同开展了简单钝体外形磁控热防护试验[11]，发现采用螺线管磁场，可使平头圆柱热流下降85%；国防科技大学李开和刘伟强[12]针对简单钝体外形开展了磁控热防护系统建模研究，发现增大螺线管半径有利于提高磁控热防护效果，而螺线管长度对磁控热防护效果影响相对较小。

尽管国内外对磁控热防护系统进行了不少研究，热防护效果得到了有效验证，但很少见到磁控热防护系统在实际复杂外形飞行器的工程应用。这主要是由于高超声速流动磁流体力学控制涉及高温气体化学反应效应、热力学效应、等离子体碰撞与迁移机制、带电粒子与电磁场相互作用、电磁场交叉感应等多种物理效应或现象[13-16]，每一种效应或现象包含复杂的变化规律和多种影响因素，对磁控热防护效果造成不容忽视的影响；各种物理机制的特征尺度或特征时间存在的差异较大，多种物理现象耦合时会出现强烈的刚性问题和收敛性问题，极大地增加了仿真模拟难度。基于此，国内外高超声速流动磁流体力学控制研究大多针对简单外形，如球头或钝体外形的一维/二维或轴对称磁流体数值模拟；物理模型方面也多以简化模型为主，如完全气体或平衡气体模型等；磁场方面多采用简化磁场，如偶极子或均匀磁场等，极少见到复杂外形飞行器磁控热防护系统仿真研究。

本文综合考虑高超声速流动磁流体力学控制涉及的多种物理效应或现象，通过耦合求解电磁场泊松方程和带电磁源项的高温热化学非平衡流动控制方程组，搭建磁控热防护仿真模拟平台。结合美国航天飞机 “哥伦比亚”号近似外形和多种磁场配置方案(均匀磁场、偶极子磁场、螺线管磁场及多个磁场组合)，较为系统地开展磁控热防护系统在滑翔返回式天地往返运载器上的应用仿真研究。

1 等离子体生成机制及其电导率计算

(1)

基于电子/离子/中性粒子扩散、碰撞与迁移的原理，多电离组分等离子体电导率可写为[13]

(2)

(3)

当s为中性组分O2、N2、NO、O和N时，其表达式为[10,14]

(4)

式(3)和式(4)中电子温度Te的模拟，需要考虑等离子体流动过程中的热力学效应，气体分子平转动温度和振动-电子温度之间的能量输运可写为[19]

(5)

2 电磁场对流场动量和能量输运过程的模拟

具有弱导电性的等离子体在磁场中流动时，会产生感应电流，使流体微元受到洛伦兹力作用，这会改变气体的动量和能量。与此同时，感应电流还会产生焦耳热作用，进而改变气体内能等状态参数。由于在绝缘壁面条件下，霍尔效应对磁控效果的削弱影响相对较小[9-10,20]，因此为了增强磁控效果，飞行器表面可采用绝缘材料(或涂层)，此时忽略霍尔效应，电磁场对流场动量和能量输运可写为

(6)

式中：F为洛伦兹力；J为感应电流；B为磁感应强度；U为速度矢量；F·U为洛伦兹力做功；σ-1J2为焦耳热。由于电磁焦耳热的部分能量作用于振动-电子温度模态，因此电磁场对振动-电子温度模态的能量输运可写为

EMV=γσ-1J2

(7)

式中：γ为焦耳热振动-电子能量配比。

感应电流J可通过离散求解各向同性电场泊松方程计算[16]：

(8)

磁感应强度按照磁场矢量叠加原理计算：

B=Bext+Bin

(9)

式中：Bext为外加磁场，由机载磁场发生装置产生；Bin为感应磁场(或称诱导磁场)，由感应电流诱导产生。对于磁雷诺数较低的状态有|Bin|≪|Bext|，因此可近似认为B=Bext。而对于磁雷诺数较大的状态，Bin通过离散求解磁矢量势泊松方程计算[10]：

(10)

3 机载磁场Bext的模拟

机载磁场的模拟，最常见的有均匀磁场、偶极子磁场和螺线管磁场。其中均匀磁场模拟较为简单，一般是在飞行器流场某一区域，直接附加强度和方向相同的磁场。偶极子磁场是螺线管磁场的简化模型，其表达式一般可写为[15]

(11)

式中：(r,α)为极坐标单位矢量；B0为极轴上距离偶极子中心r0处磁感应强度。

对于实际的螺线管线圈磁场，其计算较为复杂，计算量较大。但如果忽略电流分布的非均匀性和线圈轴向的螺旋性，不仅可以极大地减少计算量，而且计算与实际测量结果之间的差异极小[21]。此时，通电电流I、线圈密度n、直径2a、长度2b的螺线管磁场可写为[21-22]

(12)

式中：坐标原点为螺线管中心；Br、Bφ和Bz为磁场径向、周向和轴向分量，部分变量表达式为

其中：A为特征强度；k、k1、k2和h均为拟合函数中间变量；f(k)、fz(r,z)分别为r方向和z方向磁场拟合函数；K(k)、F(h,k)为第一类和第三类椭圆积分。

4 高超声速电磁流动控制方程及离散

高超声速电磁流动控制方程与一般的高超声速热化学非平衡Navier-Stokes(N-S)方程的主要区别在于包含电磁作用源项，其无量纲形式为[14]

(13)

式中：守恒变量Q=[ρi,ρ,ρu,ρv,ρw,ρEt,ρEv]T,ρi为组分i气体密度,u、v、w分别为x、y、z3个方向速度分量，Et为气体内能，Ev为气体振动-电子能；Re为雷诺数；F、G、H与FV、GV、HV分别为x、y、z方向无黏向量与黏性向量；W和WMHD分别为非平衡源项和电磁作用源项，其形式为

(14)

其中：fx、fy、fz为洛伦兹力F3个方向分量。

对于式(13)，采用AUSMPW+(Advection Upstream Splitting Method by Pressure-based Weight functions)格式离散无黏项，采用中心差分格式离散黏性项，采用全隐式全耦合LU-SGS(Lower-Upper Symmetric Gauss Seidel)进行时间离散。采用MPI(Message Passing Interface)并行技术，应用负载平衡、非阻塞发送/阻塞接收、整合打包/解包以及对等/主从模式相组合等方法，建立流场分区和功能分解结合的并行算法，提升计算效率。

5 数值验证

采用钝锥体外形[23]对数值计算进行校验。头部半径为0.09 m，椎体的锥角为15°，全长为0.122 4 m；计算飞行高度40.0 km、马赫数15.0的状态，等温壁面条件。采用与文献[23]相同气体模型和磁偶极子配置，偶极子中心为头部圆心，B0=2.0 T,r0=0.09 m。M1偶极子方向为X轴负方向；M2偶极子方向为Y轴正方向，见图1(a)。

图1(b)给出了磁场配置对钝锥体表面热流Q影响，S为表面弧长。可见，3种磁场配置条件下，计算结果均与文献[23]符合，这说明本文数值模拟方法可以较为准确地预测和评估高超声速磁控气动热防护效果。

图1 磁场矢量和钝锥表面热流

6 磁控热防护数值仿真

6.1 计算状态、网格及气动热环境模拟

以航天飞机(OV-102)近似外形为研究对象，以飞行器头部尖点为坐标原点，飞行器长为X轴方向，飞行器高为Y轴方向，飞行器横向为Z轴方向，建立直角坐标系。采用多块结构网格，在壁面附近以及形变较大区域，如竖直舵、水平翼等区域，进行局部加密，如图2(a)所示。计算状态采用典型轨道参数条件，见表1，这些状态均有飞行试验热流结果。

表1 航天飞机飞行状态

由于航天飞机使用表面隔热瓦，相对于原子复合和有关的离解能量调节是非催化的[24]，因此本文数值模拟时飞行器表面采用辐射平衡温度条件和壁面非催化条件。

图2(b)给出了两套不同疏密度网格计算得到表面热流分布与试验结果的比较，飞行器全长L=30.24 m。Grid1为稀网格，壁面法向第1层网格间距为10-5m，网格总量约为1 200万；Grid2为密网格，第1层间距为10-6mm，网格总量约为5 000万。可见，计算结果与飞行试验[24]符合较好，这说明本文方法能较好模拟航天飞机的气动热环境；两套网格计算差异较小，满足网格无关系要求。为了保证流场的分辨率，采用Grid2进行计算。

图2 网格示意图和航天飞机迎风中心线热流

6.2 数值模拟结果分析

由于航天飞机头部迎风面和翼前缘等表面区域热流较高，为了降低局部区域热流，可在这些区域加载磁场发生装置，包括偶极子磁场、螺线管磁场和均匀磁场以及多个磁场组合。这里将位于X-Y平面内与X轴负方向的逆时针夹角定义为θ；位于Y=Constant平面与X轴负方向的逆时针夹角定义为φ。外加磁场配置如下：

Case1采用单个螺线管，主要用于头部驻点区表面热防护。螺线管中心点位于坐标(0.4,-0.273 8,0)m，该点与头部顶点的横坐标相差0.4 m，与对应横坐标的下表面点距离也为0.4 m，利于螺线管磁场物理上调整方向。螺线管轴线正方向位于X-Y平面内，与X轴负方向的逆时针夹角θ=0°,20°,40°,60°。由于增大螺线管直径有利于提高磁控热防护效果，而螺线管长度对磁控热防护效果影响相对较小[12]，因此这里尽量增大螺线管直径，结合防热层设计和螺线管磁场灵活转向要求，这里取螺线管直径为0.4 m、长度为0.4 m。通过调节螺线管通电电流和匝数，可调节磁感应强度，使螺线管轴线正方向距离中心点0.4 m处，磁感应强度为0.3 T，这在工程上可以实现[6,11]。

Case2单个偶极子磁场，主要用于与Case 1螺线管磁场对比。其中心点和方向设置均与Case 1相同。特征长度r0为0.4 m，特征磁场感应强度B0为0.3 T。

Case3多个螺线管磁场组合，主要用于大面积防热。螺线管1～3中心点分别位于坐标(0.4,-0.273 8,0)m、(0.8,-0.445 6,0)m和(1.6,-0.703 8,0)m，与对应横坐标的下表面点的距离为0.4 m；其轴线正方向均位于X-Y平面内，与X轴负方向的逆时针夹角分别为θ1～θ3。

Case4采用均匀磁场，在迎风面X=0.4～1.0 m、Y<0 m、Z=0～0.6 m的区间内，设置均匀磁场，磁感应强度方向位于X-Y平面内，与X轴负方向的逆时针夹角θ=0°～180°。需指出的是，该方案磁场边界上不满足磁场无源、无散度条件，这里仅作洛伦兹力的对比分析之用。

Case5采用单个螺线管，主要用于水平翼前缘热防护。中心点位于坐标(21.228,-0.685 9, 5.84)m。螺线管轴线正方向单位矢量为(-0.707,0,0.707)；考虑翼前缘尺寸，这里螺线管直径和长度均设为0.1 m；通过调节螺线管通电电流和匝数，使螺线管轴线正方距离中心点0.1 m处，磁感应强度分别为0.3 T和1.0 T。

1)不同飞行状态和磁场方向差异

图3(a)给出了有/无螺线管磁场条件下不同飞行高度迎风面中心线热流比较。这里飞行器全长取为L=32.35 m；采用Case 1的螺线管磁场配置，螺线管轴线正方向与X轴负方向的逆时针夹角θ=0°。可以看出，在头部设置螺线管磁场可以有效减低驻点附近的热流；磁场还会使飞行器迎风面靠近底部区域的热流略有上升；飞行高度越高，磁场使头部热流下降幅度越大，高度H=70.2 km 时驻点热流下将约26.7%。

图3(b)～图3(d)给出了不同飞行高度电离气体电导率分布。可以看出，电导率较高区域均在头部迎风区域附近，这说明将Case 1磁场设置于头部迎风区域是合理的；高度H=57.2 km流场中气体电导率明显低于其他状态，这与图3(a)该状态的磁控热防护效果较弱相符，这可能是57.2 km时飞行速度相对较低造成的(表1)。

图3 不同飞行高度表面热流和流场电导率分布

图4(a)给出了Case 1螺线管磁场不同磁场方向时迎风面中心线热流。可以看出，单一螺线管磁场方向大幅变化，对磁控热防护效果存在一定影响；随θ增大，热流凹陷低点后移，但整体的定性规律变化不大，均表现为驻点区热流显著下降，下游迎风面局部区域出现不同程度的热流反冲现象[16]。

图4(b)给出了螺线管磁场Case 1不同磁场方向时Z方向电流(云图颜色)和洛伦兹力矢量(黑色箭头，长短表征洛伦兹力大小)分布，红色箭头为来流方向。可以看出，螺线管磁场方向大幅变化，环形电流和洛伦兹力分布存在一定差异；随θ增大，洛伦兹力交汇点略微后移，但其整体特征变化不大，这与图4(a)的热流变化规律相符。

图4 不同螺线管方向表面热流和流场洛伦兹力

2)偶极子和螺线管磁场差异

图5(a)给出了螺线管磁场(Case 1)和偶极子磁场(Case 2)磁场矢量图和等值线图。可以看出，当特征长度和特征磁场强度相同时，两者存在一定差异，但差异幅度较小。图5(b)进一步比较了螺线管磁场(Case 1)和偶极子磁场(Case 2)的磁控热防护效果差异。可以看出，两种磁场使表面热流下降的整体规律基本一致，在局部细节上存在一定差异。这说明采用偶极子磁场，可在一定程度上近似模拟螺线管磁场的分布及其影响；但如果要精细地分析磁流体力学控制的效率，则需采用接近物理实际的磁场模型。

图5 磁场差异及其对热防护影响

3)多个螺线管磁场组合

由图3和图4可以看出单个螺线管磁场可以使头部热流有效下降，但下游区域热流下降幅度不大，甚至出现热流反冲现象。为了进一步减小迎风面较大区域热流。可采用多个螺线管组合磁场的方式。

图6(a)给出了采用螺线管组合磁场Case 3得到的迎风面中心线热流。这里θ1=0°，由图4可知此时局部热流反冲较小；θ2和θ3的值相等(即图中θ2，3)，分别取为45°、90°和135°。可以看出，当θ2，3=135°时，可降低头部迎风面较大区域的热流，且热流反冲现象较小。

图6 组合磁场条件下的表面热流和洛伦兹力

图6(b)给出了采用螺线管组合磁场Case 3时压力和洛伦兹力(黑色箭头，长短表征洛伦兹力大小)分布。可以看出，当θ2，3=45°时，在X=0.5 m 附近，洛伦兹力使高温气流向壁面附着，形成高压、高热流区；而当θ2，3=135°时，由于磁场分布变化，洛伦兹力使高温气流向壁面附着的区域前移至X=0.4 m附近，且该洛伦兹力矢量相对较小，因此热流反冲区域略微前移，且反冲幅度较小。

4)磁场方向对洛伦兹力分布和流场结构的影响

为了进一步分析磁场方向对流动的影响，图7(a)给出了均匀磁场条件下(Case 4)的磁场矢量(红色箭头)、压力云图和洛伦兹力矢量(黑色箭头)分布。可以看出，这种情况下，无论磁场方向怎么变化，洛伦兹力矢量永远是垂直于磁场矢量方向的，其作用效果永远是阻碍流动方向的。当θ∈[0°, 60°)时，洛伦兹力使高温气流附着，随θ增大，附着的角度减小；当θ∈(60°, 160°)时，洛伦兹力使高温气体远离壁面外推，随θ增大，外推角度增大。

进一步结合洛伦兹力强度分析，在忽略感应电场影响时，JZ∝(U×B)Z，因此洛伦兹力|F|=|J×B|∝JZ∝(U×B)Z。当θ∈[0°，60°)时，随θ增大，速度矢量U与磁场矢量B夹角增大，因此JZ增大，洛伦兹力强度增大。但由于附着角度减小，洛伦兹力使其流向壁面的压缩程度反而减小。当θ∈(60°，160°)时，随θ增大，外推的角度增大，洛伦兹力越来越倾向于使高温气流远离壁面方向。同时，速度矢量与磁场矢量之间的夹角逐渐减小，会使相同条件下的洛伦兹力强度减小，因此在θ≈120°时，磁场的外推作用最强。

图7 磁场方向对流动结构的影响

5)水平翼前缘热环境的磁控

图8(a)给出了高度H=68.2 km时的表面热流分布。可以看出头部驻点区域并不是表面热环境最严酷的区域，水平翼部分区域的表面热流更高，因此针对飞行水平翼的高热流区设计了磁场配置方案Case 5。水平翼的厚度较小，螺线管磁场发生装置尺寸不能太大，其长度和直径均取为0.1 m；螺线管中心位于高热流峰值所在的区域附近Z= 5.84 m平面内，见图8(b)。

图8 全表面热流分布和Case 5磁场设置

图9给出了磁场对水平翼前缘热流的影响。可以看出，当磁场特征强度为 1.0 T时，表面热流下降很明显，这说明采用合适的磁场可以有效地降低水平翼前缘的热流；而当磁场特征强度为0.3 T时，热流受磁场影响不明显。

图9 不同磁场条件下水平翼前缘热流

为了进一步分析图9的结论，图10分别给出了磁场配置Case 5时3种情况下水平翼附近的电导率分布。这里显示的平面经过方案Case 5螺线管中心点，平面方向为(0.707，0，0.707)m。

由图10可以看出，当B0= 0.3 T时，水平翼前缘附近气体电导率较低；当B0=1.0 T时，电导率显著升高。这可能是由于随磁感应强度增大，磁场使气体减速的效果增强，气体速度下降，更多的气体动能转化为内能，从而使气体温度上升，高温区增大，气体离解和电离程度增加，电导率增大；而电导率增大，又反过来增强了磁场对气体滞留和减速等磁控效果，由此形成“正反馈”。当磁场强度较小时，这种正反馈机制强度较小，表面热流变化不明显；当B0=1.0 T，正反馈效应较强，电导率明显升高，因而使表面热流显著下降。由此可见，磁场的作用效率与流动中等离子体生成过程紧密相关，准确模拟时必须予以综合考虑。