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

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

高超声速飞行器飞行过程中，强烈的气动加热使激波层气体发生热电离反应，气体能量模态激发，会形成具有弱导电性等离子体流场[1-2]。利用机载磁场发生装置向等离子体流场注入适当的动量和能量，可对高超声速飞行器等离子体流场进行控制，这就是高超声速磁流体(Magneto Hydro Dynamic, MHD)控制。它可以有效地提升飞行器的气动性能，在高超声速飞行器气动力操控、气动热减缓和等离子体分布调节等方面具有广阔的应用前景[3-4]。

影响高温电离混合气体电导率的因素很多，要得到较为准确的气体电导率分布，首先依赖于高超声速等离子体流场的准确获得,包括气体电离程度、组分、压力和热力学温度状态(主要是平转动温度和振动-电子温度)等参数精确模拟。由于等离子体是在高超声速流动中产生的，这就存在化学反应与热力学松弛的非平衡现象。在高超声速磁流体控制数值模拟时，考虑高温气体热化学非平衡效应得到详细的等离子体流场参数，是电导率准确获得的基础。

从高超声速磁流体控制研究发展来看，由等离子体流场分布得到等离子体电导率，国内外常见的处理方法主要有3种形式：定电导率方法、试验拟合经验公式(或半经验公式)和电离组分迁移碰撞模型。

定电导率方法[4]假设流体的电导率均为试验测得或人为假定的某一定值(或无穷大)，不随流场变化而变化。这种处理较为简便，缺点在于它与真实的高超声速等离子体流场电导率分布差别较大，大多用于磁流体控制定性分析或某些特殊流动中，如电导率变化不大的低速流动、超导流体(电导率无穷大)或近似超导流体的高电导率流体等。

试验拟合经验公式方法[6](或半经验公式)，一般是通过试验得到随温度、压力等参量变化的电导率测量值，进而得到电导率的拟合计算公式或半经验公式，这种方法极大地依赖试验数据的覆盖面和准确性。目前，电导率的试验数据，主要来源于国外平衡态气体的试验测量，这与高超声速热化学非平衡状态的等离子体电导率存在差别，且不同试验给出的测量结果差异较大[7]，导致了拟合公式的多样性，其准确性也需要进一步讨论分析。

电离组分迁移碰撞的电导率模型，主要基于分子运动碰撞理论，考虑了电子、离子等在等离子体中碰撞和迁移，得到混合气体电导率。相比于固定电导率和拟合经验公式方法，它不仅更贴近等离子体混合气体导电作用机理，而且能综合考虑混合气体电离度、密度、压力、气体组分、热力学状态等多种因素，因而能较好地适用于高超声速飞行器高温气体非平衡等离子体流场，得到广泛应用[7-10]。但由于等离子体混合碰撞导电机理尚不完备，尤其是模拟所需的某些特征量(如电子与中性粒子有效碰撞截面等)很难准确得到[11-13]，这客观影响了电导率计算的准确性。

近30年来，随着高超声速飞行器磁流体控制数值模拟的兴起，人们逐渐认识到等离子体电导率模拟的准确性问题。例如，2009年，Bisek和Boyd[14]对比分析了Raizer、Otsu、Spitzer-Harm、Chapman-Cowling等提出的电导率模型以及电离组分迁移碰撞电导率模型，发现各种模型计算得到的电导率差异很大，在高速流场中，其分布变化规律存在很大的不同，在5 000～10 000 K区间，与非平衡流动试验测量结果相差很大。2010年，Bisek和Boyd[15]采用3种电导率模型计算的钝体激波脱体距离，不同电导率模型差异明显，且均与试验结果差别较大。2012年，Cristofolini和Borghi[16]对意大利CIRA使用的两款数值模拟软件CAST和EMC3NS中使用的3种电导率模型进行了对比分析，发现不同电导率模型计算结果差异明显，对磁流体流场分布的影响十分显著。

国内也有不少高超声速飞行器磁流体控制方面的研究。2005年，赫新等[17]数值模拟了理想磁流体一维MHD激波管流动和二维喷管MHD流动，其磁感应方程磁扩散项为零(即电导率为无穷大)；2007年，潘勇[3]采用定电导率方法对完全气体磁流体控制进行数值模拟研究，考虑化学反应对无黏磁流体(磁扩散项为零，电导率无穷大)进行了数值模拟研究；2008年，田正雨[5]采用试验拟合电导率对钝体化学平衡气体和完全气体磁流体热流控制进行了数值模拟分析；同年，陈刚等[18]采用Bityurin给出的温度拟合电导率模型对完全气体二维磁流体进行了数值模拟；2009年,黄富来和黄护林[19]采用Kantrowitz给出的电导率模型和7组分化学反应模型对高超声速磁流体控制进行了数值计算研究；2013年，黄浩等[20]采用完全气体模型和定电导率法对电子束电离的高超声速磁流体发电机进行了研究；同年，何淼生等[21]采用电子束电离电导率模型和完全气体模型研究了超高声速进气道磁流体控制；2016年，李开和刘伟强[22]采用完全气体模型和Raizer给出的电导率经验公式对高超声速飞行器常规螺线管磁控热防护系统进行了研究，2017年李开等又采用Fujino给出的电导率计算公式开展了磁控热防护霍尔电场数值模拟分析[23]；2018年，姚霄等[24]采用完全气体模型和Raizer给出的电导率经验公式对5种外磁场作用的磁流体增阻特性进行了分析。

可以看出，尽管国外高超声速磁流体控制研究得到了很大的发展，但等离子体电导率准确性及其影响问题，一直都困扰着高超声速MHD数值模拟。国内高超声速磁流体控制研究，在气体模型方面大多采用简化气体模型，仅有少量的考虑高温非平衡效应(主要是化学非平衡，很少考虑热力学非平衡效应)的高超声速磁流体控制研究；在电导率模型方面大多直接采用国外某一种电导率模型，忽略其准确性问题带来的影响，缺乏电导率准确性及其对磁流体控制影响的分析研究。

作者所在研究团队对高超声速飞行器非平衡等离子体流场及磁流体控制，进行了较为广泛的研究[25-27]。本文在此基础上，主要针对电导率准确性及其影响问题，基于高温热化学非平衡效应和低磁雷诺数假设[3-4]，开展高超声速磁流动控制数值模拟，结合国内外常用的电导率处理方法及模型，分析气体电导率对等离子体电磁控制流场特性及其气动力热特性作用机理和影响规律，结合国外部分试验数据探讨电导率准确性问题。

1 数值计算方法

高超声速飞行器绕流流场中，混合电离气体一般符合低电导率特征，满足低磁雷诺数(Rem≪1)假设[3-4]。此时感应磁场相比于外磁场来说，基本可以忽略，无量纲化控制方程可写为

W+WMHD

(1)

其守恒变量为

Q=[ρj,ρ,ρu,ρv,ρw,ρEt,ρEv]T

式中：ρ为混合气体密度；ρj为混合气体中第j组分的密度；u、v、w分别为直角坐标系x、y、z3个方向的速度；Et为混合气体内能；Ev为混合气体振动能；Re为雷诺数；F、G、H与FV、GV、HV分别为直角坐标3个方向的无黏通量与黏性通量；W和WMHD分别为热化学非平衡源项和电磁作用源项，其中电磁作用源项可写为

(2)

式中：J为电流密度；B为磁感应强度；E为电场强度；γ为焦耳热振动能量配比；σ为混合气体电导率。兼顾计算精度和效率，式(1)无黏项采用AUSMPW+(Advection Upstream Splitting Method by Pressure-based Weight functions)格式，黏性项为中心差分格式，时间离散为LU-SGS(Lower-Upper Symmetric Gauss Seidel)隐式方法。感应电流密度和电场强度，通过式(1)耦合广义欧姆定律与电场泊松方程数值求解：

(3)

式中：φ为电势函数。

2 物理化学模型

2.1 热化学非平衡模型

考虑高温流场中离解、电离、置换、复合等化学反应，采用较为常用的11组分空气化学反应Park模型[28]。考虑热力学非平衡效应时，采用热力学双温度模型[25]。

2.2 平衡气体模拟

与本文热化学非平衡模型计算方法一致，平衡气体模拟同样采用11组分空气化学反应Park模型。通过迭代求解化学反应平衡方程、电荷守恒方程、氧元素守恒方程和氮元素守恒方程，得到气体在当地温度、压力条件下的平衡组分。

2.3 电导率模型

采用定电导率方法时，全流场电导率为某一固定值，不受流动变化影响。

对于平衡态空气，电导率主要受温度和压力影响，其中气体温度对电导率影响程度远大于气体压力，这里给出几种基于温度T的电导率半经验公式：

1) 电导率模型M1，由Raizer基于试验测量给出，应用较为广泛[14-16,22,24]：

σ=8 300.0e-36 000.0/T

(4)

2) 电导率模型M2，Otsu等[29]给出基于温度的电导率半经验：

σ=300.0×(T/18 000.0)4

(5)

3) 电导率模型M3，Nagata等[30]采用基于温度的半经验公式：

σ=250.0×(T/25 000.0)2

(6)

对于等离子体态电离气体，电导率不仅受温度、压力等状态参数影响，而且与气体电离度密切相关，这里给出几种考虑气体电离度的电导率模型:

1) 电导率模型M4， Macheret和Shneider认为在电离度较低时，电导率与电离度Xe呈近似正比关系[13]：

σ=2.7×105XeXe≤10-2

(7)

2) 电导率模型M5，Chapman和Cowling等基于玻尔兹曼分布和二元扩散模型，发展了弱电离气体电导率模型[14-15]：

(8)

3) 电导率模型M6，Spitzer-Harm发展了强电离气体电导率模型[14-15]：

(9)

4) 电导率模型M7，Kantrowitz结合Chapman-Cowling弱电离气体电导率模型和Spitzer-Harm强电离电导率模型，发展了Kantrowitz模型[5]：

(10)

式中：Xe为电离度；Qs为有效碰撞截面；ne为电子数密度。

对于含多种离解与电离成分的高温非平衡混合气体，即使温度、压力和电离度相同，混合气体的电导率也会存在差别。国外常用多组分等离子体电导率模型[7-10]是基于电子和离子在等离子体混合气体中碰撞和迁移模型，该模型不仅综合考虑了弱电离和强电离导电机理，而且考虑了不同气体组分的差异。其表达式为

(11)

(12)

当s为中性组分，其表达式为

(13)

2.4 磁流体控制气动力特性计算方法

一般情况下，气动力系数计算只需考虑飞行器表面的压力和黏性应力积分，而高超声速磁流体控制气动力系数计算时，还需考虑洛仑兹力反作用力的空间体积积分，以二维磁流体阻力系数CD计算为例：

(14)

式中：FD为飞行器受到的阻力；F为表面压力和黏性应力造成的阻力；FM为洛仑兹力反作用力造成的阻力(以下文中简称洛仑兹磁阻力)；α为飞行迎角；q∞和Sref分别为动压和参考面积。

3 数值方法验证

在以往工作中，作者所在的团队对数值模拟方法，从非平衡等离子体流动[2,25]、气动热环境模拟[31]以及磁流体控制[27]等多个方面进行了考核验证。这里仅补充高超声速磁流体控制气动力特性方面的校验。

采用球柱外形[30]，半径为1 m，柱长10 m。计算条件为：飞行高度69 km，速度6 500 m/s。采用磁偶极子磁场，磁场分布为

(15)

式中：(r,θ)为极坐标单位矢量，磁偶极子中心位于头部球心(坐标原点)；B0为磁感应特征强度，B0=0～0.2 T；r0=1.0 m。

图1为磁阻力伞形成示意图。图2为不同磁感应强度球柱阻力系数，标号CD为阻力系数，CD1为压力和黏性应力阻力分量，CD2为洛仑兹力反作用力阻力分量。可以看出，对于球柱外形，弱导电性的等离子体在磁场中流动产生环形电流，环形电流受到与流动方向相反的洛仑兹力，使激波外推，形成磁阻力伞；磁感应强度增大，洛仑兹力磁阻力分量系数增大，使总的阻力系数增大；本文计算结果与文献[30]符合较好。

图1 磁场阻力伞示意图Fig.1 Schematic of magnetic field drag

图2 阻力系数随磁场强度变化Fig.2 Variation of drag coefficient with magnetic strength

4 数值计算分析

计算采用RAM-C钝锥外形[32]，头部半径为0.152 4 m，全长为1.295 m，锥角为9°。计算飞行高度71 km，速度7 650 m/s，壁温1 500 K。采用磁偶极子磁场(式(15))，磁感应特征强度B0=0.5 T,特征长度r0=0.152 4 m。磁偶极子方向为直角坐标横轴负方向。

为了排除网格对计算的影响，图3给出了两套不同密度网格计算结果。Q为表面热流，X为横坐标；Ci为组分质量分数，R为离开驻点的距离；图中曲线标识O Grid1为采用Grid1计算得到氧原子组分，标识N Grid1为Grid1计算得到氮原子组分，其他标识含义类似。Grid1为稀网格，壁面第1层法向间距为0.005 mm；Grid2为密网格，第1层法向间距为0.001 mm。可以看出，两套网格计算结果重合，这说明网格对计算结果影响较小，文中如无特殊说明，均采用密网格Grid2。

图3 不同网格计算结果Fig.3 Results under different grids computation

图4给出了热化学非平衡流场电子数密度ne数值模拟结果与飞行试验结果[32]对比，曲线NCW为采用完全非催化壁面条件的计算结果，FCW为采用完全催化壁面条件的计算结果，Rn为球头半径。可见，数值模拟结果与飞行试验结果吻合较好，这说明本文采用的计算方法能较好地模拟热化学非平衡等离子体流动，满足磁流体控制电导率准确模拟的基础需求。

图4 等离子体电离环境模拟Fig.4 Simulation of plasma ionization environment

4.1 电导率对磁流体控制影响

采用定电导率方法，开展磁场作用对飞行器流场特性及其气动力热特性作用影响分析。

图5给出了不同电导率条件下流场驻点线温度和主要组分N2质量分数分布，图中标识No Mag.为无磁场作用结果。可以看出，磁场作用使激波脱体距离外推，其外推距离与电导率呈近似比例关系，电导率越大，激波脱体越远；磁场作用下，流场高温区(6 000 K以上)随电导率增大而显著增大，但温度峰值随电导率增大变化规律不明显；高温区增大，化学离解反应有效时间变长，主要气体组分N2离解程度增大。

图6(a)给出了不同电导率条件下阻力系数变化，由图可以看出随电导率增大，洛仑兹力反作用力阻力系数分量CD2逐渐增大，进而使总的阻力系数CD增大；由表面压力和摩擦力形成的阻力系数CD1随电导率增大，呈略微下降趋势，下降幅度较小。为了进一步分析洛仑兹力反作用力阻力系数分量CD2随电导率变化，图6(b)和图6(c)分别给出了流场环形电流密度大小与洛仑兹力反作用力阻力分量云图，图上半部分电导率为100 S/m，下半部分电导率为300 S/m。可以看出，电导率越大，环形电流密度越强，磁阻力伞越大，高洛仑兹力阻力区域越大，综合使CD2增大。

图5 驻点线温度和氮气质量分数分布Fig.5 Distribution of temperature and mass fraction of N2 along stagnation line

图7(a)给出了不同电导率条件下表面热流沿轴线的分布，可以看出，磁场作用使表面热流下降，其下降幅度受电导率影响，但影响程度远小于电导率对阻力系数的影响。当电导率大于200 S/m时，表面热流几乎不受电导率影响，出现电导率主导的“磁控热饱和”现象，其中驻点区域这一现象更明显。为了进一步分析该现象产生的原因，图7(b)～图7(e)分别给出了驻点线N原子质量分数、O原子质量分数、化学焓和气体总焓分布。结合图5和图7可以看出，随电导率上升，波后高温区增大，有效化学反应时间逐渐变长，气体离解度增加，趋于化学反应平衡态。壁面附近N、O原子质量分数逐渐升高，其变化幅度逐渐变缓，趋于某一定值，这导致壁面化学焓的增大幅度逐渐减小，进而使壁面总焓趋于定值。由驻点热流公式Fay-Riddell公式[2]可知，在类似条件下，热流与焓值差成正比：Q∝Hs-Hw(Hs为边界层外缘总焓，Hw为壁面总焓)，结合图7(a)和图7(e)可以看出，磁场使热流的减小幅度，基本与焓值差的减小幅度相当。由此可见，磁场作用使化学离解程度增强，焓值差减小，可能是磁控热流下降的原因之一；而随电导率增加，流场中的化学反应逐渐趋于平衡态，离解度趋于定值，可能是热流下降幅度趋于饱和的原因之一。这里只根据本文低电导率特征条件下的计算状态进行分析，其深层次机理，还需结合更多状态的数值模拟和理论分析给出。

图6 不同电导率条件下阻力特性Fig.6 Drag characteristics under different electrical conductivity conditions

图7 不同电导率条件气动热环境分析Fig.7 Aerodynamic thermal environment analysis under different electrical conductivity conditions

4.2 电导率模型差异的影响及其准确性分析

针对本文给出的9种电导率计算模型，展开数值对比分析。图8为无磁场作用时采用不同电导率模型计算的驻点线电导率、温度、振动温度(Tv)和电离度分布。由图可以看出，不同电导率模型计算得到的流场电导率差别非常大。模型M1～M3虽然均采用温度作为气体电导率的唯一变量，但局部区域计算结果相差一个数量级以上，这种不确定性与模型的适用范围存在关联(例如M1仅适用于1 atm(1 atm=105Pa)下弱电离的空气、氮气和氩气)，因此这类模型在高超声速流动中(压力随流场结构变化较大、电离度受非平衡效应影响具有不确定性)的适用性存在很大限制；模型M4～M9，考虑了气体在高温下电离效应，其差异主要来源于对气体导电机理的认知不同，例如M4为弱电离条件下电导率与电离度近似正比的经验公式；M5主要考虑弱电离条件下电子与中性粒子的碰撞对电导率的影响占主导；M6主要考虑强电离条件下电子与离子的碰撞对电导率的影响占主导；M7综合考虑弱电离(模型M5)和强电离(模型M6)产生的电阻效果的串联；模型M8和M9，可看作M7模型的修正形式，进一步考虑了电子与多种中性粒子碰撞的有效截面差异和热力学平动-振动温度松弛对电子温度影响的差异。由于高超声速外流场中气体电离度一般较低，因此M6模型的适用性相对较差；M5和M7考虑单一中性粒子与电子的碰撞，一般仅适用于单一组分气体电导率模拟；M8和M9模型考虑的导电机理相对完备，因此能较好地适用于高温非平衡多种组分混合气体。

图8 不同电导率模型驻点线参数Fig.8 Parameters along stagnation line under different electrical conductivity models

图9给出了采用不同电导率模型数值模拟的驻点线温度分布和表面热流轴向分布。表1给出了不同电导率模型计算的驻点热流Q0和阻力系数及其磁控效率。可以看出，电导率模型差异，直接影响了磁流体的数值模拟控制效果：采用不同电导率模型，磁场控制定性效果相似(激波脱体距离增加、热流下降、洛仑兹力导致的磁阻力增强)；其中激波脱体距离和磁阻力的增大程度，受电导率模型影响较大，表面热流下降程度受电导率模型影响相对较小，这与4.1节磁控热饱和现象的分析一致；值得指出的是，电导率模型导致的电导率分布差异，不等同于磁控效果差异，例如采用M2、M3和M7计算得到的电导率分布差异明显，但其驻点热流、阻力系数控制效果差别较小。

图9 不同电导率模型磁控效果Fig.9 Efficiency of magnetic control under different electrical conductivity models

表1 不同电导率模型计算得到的驻点热流和阻力系数
Table 1 Heat flux at stagnation point and drag coefficient under different electrical conductivity models

模型Q0/(kW·m-2)D1/%CDCD1CD2D2/%NoMag.880.40.2920.2920M1548.0380.4750.2670.20963M2683.0220.3340.2780.05614M3722.5180.3340.2860.04714M4517.7410.4370.2690.16850M5696.8210.3270.2840.04312M6433.2510.9890.2360.754239M7707.4200.3260.2840.04112M8464.2470.4890.2680.22167M9516.0410.4390.2730.16650

注：D1为磁控驻点热流下降百分比；D2为磁控阻力增加百分比。

为了进一步分析本文给出的电导率模型，图10(a)给出了1个大气压下采用不同电导率模型计算的平衡气体电导率与试验结果[33]的比较，图10(b)为该条件下电离度分布。可以看出：当电离度较大(Xe>10-2)时，M6、M8、M9电导率与试验结果符合较好，此时电子与离子的碰撞对电导率的影响占主体；当电离度较小(Xe≤10-2)时，M4、M8、M9与试验结果符合较好，此时电子与中性粒子的碰撞对电导率的影响占主体；M5和M7仅考虑单一中性粒子与电子的碰撞，这可能是其与试验结果差别较大的原因；M8和M9在较低电离度时，差异相对明显，这体现了电子与中性粒子碰撞截面试验拟合准确性对电导率的影响；综合来看，电导率模型M8与试验结果符合程度最高。

由图6(a)、图8(a)和表1还可以看出，采用电导率模型M8计算得到的电导率在激波层内很大区域可达200～400 S/m，其增阻效果为67%，

反而比采用定电导率100 S/m时的增阻效果123%低得多。为了分析这一现象，图11给出了这两种条件下环型电流密度大小和洛仑兹力阻力分量分布。由图可以看出，采用定电导率方法，激波外有较大区域电流密度和洛仑兹磁阻力分量较大，而采用电导率M8模型计算的该区域电导率、环形电流和洛仑兹磁阻力分量接近于0。由此可见，对于气动加热导致离解、电离产生的高超声速等离子体流场来说，采用定电导率方法，将人为放大磁场的洛仑兹磁阻力效果，从而使阻力系数预测值偏大。

图10 平衡空气电导率与电离度(1 atm)Fig.10 Electrical conductivity and degree of ionization of equilibrium air(1 atm)

图11 定电导率方法(σ=100 S/m)与电导率模型M8条件下电磁作用效应分析Fig.11 Analysis of electromagnetic effect under given electrical conductivity (σ=100 S/m) and electrical conductivity model of M8

5 结 论

1) 开展了定电导率条件下磁控数值模拟与磁控机理分析。气体电导率直接影响高超声速MHD控制效果，一定范围内电导率与激波脱体距离、阻力系数呈近似线性关系，电导率越大，激波脱体距离和阻力系数越大；磁控热流减缓效果与电导率呈非线性关系，电导率较大时会出现电导率的磁控热饱和现象；飞行器表面不同区域，出现热饱和现象的程度存在差别；热饱和现象产生的原因可能与化学反应趋向于平衡态存在一定关系；采用定电导率方法，会人为放大磁场的洛仑兹磁阻力效果，使阻力系数的预测值偏大。

2) 开展了常见电导率模型适用性及其磁控数值模拟分析。相同的等离子体流场，不同电导率模型计算得到的电导率分布差别很大，甚至存在数量级的差异，这与电导率模型的适用范围、参数选取原则存在很大关联；由于其对磁控效果的影响不容忽视，定量分析时必须选择适用的电导率模型，对于含多种离解、电离组分的高温气体流动来说，采用基于多电离组分迁移碰撞的电导率模型M8，计算结果与试验符合程度最好。