物理化学模型对高温流场等离子体分布的影响

高铁锁，董维中，丁明松，江 涛

（中国空气动力研究与发展中心，四川 绵阳 621000）

0 引 言

高超声速飞行器在大气层中以高超声速飞行时，由于激波和粘性作用，波后气体温度超过10000K，这时高温空气发生振动能量激发、离解和电离等化学反应过程，成为部分电离气体，并在飞行器周围形成高温等离子体流场。由于飞行器周围等离子体分布直接关系到电磁波的传输特性［1－6］，因此研究高超声速飞行器流场中等离子体的产生机理和分布特性，对于高超声速飞行器的无线电通信技术的发展具有重要意义。

飞行器周围流场等离子体分布与飞行高度、速度、姿态和飞行器外形尺寸、表面材料催化特性等多种因素有关。采用求解N－S方程的数值方法对其进行模拟时，除了需要解决热化学非平衡源项引起的刚性问题外，还需要选取合适的化学反应模型，因为计算中考虑的组分越多计算越复杂，对计算资源的要求越高，而采用含组分少的简化模型有时不能全面反映高温流动中等离子体的产生机制，影响数值模拟的精准度［7］。数值模拟高温空气电离流动常用的有7组分和11组分化学反应模型［8－10］，如得到广泛应用的Blottner的7组分反应模型［8］和Dunn与Kang的11组分反应模型［9］等。由于高温电离流动问题的复杂性，目前还没有明确和统一的化学反应模型选取准则。另外，准确模拟高温电离流动还需要研究热力学非平衡效应及壁面催化和壁面温度等边界条件的影响［11－14］。高温气体热力学非平衡涉及分子振动能级之间和振动与平动能级之间的能量交换以及振动－离解耦合效应等复杂物理化学过程，对热力学非平衡效应的数值模拟一直是高温气体动力学领域的前缘课题之一［7，11，13－14］。本文针对复杂升力体外形，采用数值模拟方法，研究不同反应模型和非平衡气体模型以及壁面催化和温度条件对流场等离子体数值模拟结果的影响。

1 复杂外形战术弹头非平衡流场的计算方法

1.1 控制方程

控制方程是三维热化学非平衡N－S方程，其无量纲化形式如下［13］：

其中Q＝（ρi，ρEV，ρ，ρu，ρv，ρw，ρE）T，W＝（wi，wV，0，0，0，0，0）T，ρi是组分密度，EV是分子组分的总振动能量，wi和wV是化学非平衡源项和振动非平衡能量源项。

1.2 热化学模型

采用11组分（O2，N2，NO，O，N，NO＋，e，O＋2，，O＋，N＋）和7组分（O2，N2，NO，O，N，NO＋，e）化学反应模型，具体反应式及反应速率常数见文献［9，13］。对于正向和逆向反应速率常数为kfi与kbi，第i个反应表达式（2）对应的第j个组分的化学反应的生成源项见式（3）：

振动非平衡能量源项为两项之和，一项是平动能与振动能的交换项，另一项是化学反应中分子离解引起的振动能量的变化项，其表达式如下［13］：

1.3 壁面催化条件

图1 计算与飞行测量电子数密度分布比较Fig.1 Comparison of computation with flight experiment

1.4 数值方法

采用全隐式LU－SGS数值方法对方程（1）求解［13－15］，无粘项的离散采用Yee的TVD格式［15－16］，粘性项的离散采用中心差分方法。

2 计算结果分析

2.1 RAM－C模型绕流等离子体分布的比较

RAM－C模型为球锥体，球头半径Rn＝0.1542m，半锥角θ＝9°，长度L＝1.295m；飞行条件：U∞＝7650.0m／s，H＝61km、71km，壁面温度TW＝1500K。

采用7组分化学反应模型和两温度热力学非平衡模型，考虑完全非催化（NCW）和完全催化（FCW）壁面条件，主要考核电子数密度分布情况。图1给出了两个飞行高度的电子密度分布与文献［14］的飞行数据比较，其中图1（a）和图1（b）对应流场剖面电子数密度峰值沿轴向分布。可以看出，电子密度计算与飞行测量数据符合较好，误差只在3倍左右。飞行速度不变，随着高度的降低，流动越接近热化学平衡态，完全非催化和完全催化条件的计算结果越接近。从图1（c）可见，催化条件对X／Rn＝8.1处的电子数密度径向分布具有较大影响。总体上看，完全催化壁条件的计算结果更接近飞行测量值。

2.2 化学反应模型和非平衡模型对流场等离子体分布的影响

计算外形为升力体外形，计算高度60km，飞行马赫数为20，攻角12°。认为壁面组分满足完全非催化壁面条件，壁面温度取600K。P1和P2分别表示靠近升力体头部和底部的两个剖面位置。下面主要给出升力体对称面上迎风区的计算结果。

图2给出了采用7组分与11组分反应模型计算的壁面离子和电子数密度沿轴向分布的比较。可以看出，11组分化学反应模型的电子数密度计算值高于7组分的计算值，二者在升力体球头区域差别较大，最大相差一个量级，随着流动向下游发展，其差别越来越小。头部流场温度最高达10000K以上（见图3a），此时除了考虑NO＋离子，化学反应模型中还应考虑O＋的作用。NO＋主要由O与N原子的缔合电离反应机制产生，O＋则主要由O原子与电子及其它组分的碰撞反应机制产生。由于头部流动区域温度较高，使得O原子电离起主导作用，O电离生成O＋对电子数密度的贡献超过NO电离生成NO＋的贡献，而在下游区域流场温度逐渐降低，NO的化学电离逐渐起主导作用，7组分和11组分两种反应模型的计算结果也逐渐趋于一致。高温流动中化学电离过程除了与温度密切相关外，还与流动压力与非平衡驰豫过程有关，压力有抑制电离和离解的作用。

图3和图4给出P1和P2剖面的温度、离子组分质量分数和电子数密度分布。其中图3（a）和图4（a）给出完全气体模型（PG）和非平衡气体模型的流场温度的比较，可以看出，完全气体模型的流场温度和激波脱体距离的计算值明显高于考虑真实气体效应的非平衡气体模型的计算值。这是因为分子离解和电离等吸热反应过程使得气体流场温度降低，密度增加，而化学过程对流动压力影响很小，故使得激波脱体距离减小。从图3和图4还可以看出，在升力体的头部流场区域（P1剖面），由于流场温度很高，两种化学反应模型计算的剖面电子数密度差别很大，此时电子主要来源于O电离生成O＋的化学过程，然后以次是气体分子生成NO＋、O＋2和N＋2的化学过程。在升力体的后部（P2剖面），由于流场温度降低，二者的剖面电子数密度分布差别变得很小，此时O＋的作用较头部明显减弱。

图5给出P1和P2流场剖面的温度与电子数密度分布云图。可见，迎风面的等离子体厚度明显小于背风面的等离子体厚度，随着流动从升力体头部（P1剖面）发展到后部（P2剖面），剖面温度最大值从12144K降低到4118K，剖面最大电子数密度从5.98×1014／cm3衰减到2.49×1010／cm3，但等离子体厚度明显增加，并且在升力体对称面附近等离子体相对较厚，沿展向逐渐变薄。

图2 壁面离子及电子数密度分布Fig.2 Distribution of ions andelectron number densities on the wall

图3 P1剖面流场参数分布Fig.3 Distribution of flow parameter of P1profile

图4 P2剖面流场参数分布Fig.4 Distribution of flow parameters of P2profile

图5 剖面温度与电子数密度云图（Ns＝11）Fig.5 Temperature andelectron number densities contours

图6给出了采用化学非平衡（CNeq）和热化学非平衡模型（TCNeq）计算的头部轴线温度和电子数密度分布，可见在壁面附近区域，两种模型的计算结果一致，但在离开壁面较远的流动区域二者出现差异，化学非平衡模型计算的脱体距离较小，电子数密度衰减较快。热化学非平衡模型中考虑了分子振动内能模式的非平衡驰豫过程和振动－离解耦合效应，目前此模型还在发展之中［7，14］，对其可靠性验证还存在很大困难。

图6 非平衡模型对温度和电子数密度分布的影响Fig.6 Effects of nonequilibrimmodels on distribution of temperature and electron number densities

2.3 壁面条件对流场等离子体分布的影响

图7（a）给出了不同壁面催化条件对P2剖面电子数密度分布的影响。这里没有考虑NO＋等离子在壁面的催化复合反应，只是考虑了O和N原子在壁面的有限催化反应，并认为它们的壁面催化速率常数为5m／s。从图7（a）可以看出，有限催化（KW＝5m／s）条件下，电子数密度计算值接近完全非催化壁面的计算结果。

图7 壁面条件对电子数密度分布的影响Fig.7 Effects of wall condition on distribution of electron number densities

图7（b）给出了不同壁面温度条件对头部轴线流场电子数密度分布的影响。在完全非催化条件下，壁面温度对电子数密度的影响主要集中在壁面附近很小的区域，贴近壁面的电子数密度以及电子数密度下降的梯度随壁面温度增加而减小。这里给出的结果是在同样计算网格下获得的，改变网格的数量和分布可能对计算结果产生不同程度的影响。

3 结束语

建立了复杂外形飞行器三维流场等离子体数值计算手段，实现了多块网格并行运算功能，具备模拟复杂外形高温流场电离效应、热化学非平衡效应和壁面催化效应的能力。RAM－C模型绕流等离子体分布数值计算与飞行测量结果具有较好的一致性，表明了模型和数值方法的可行性。通过研究热化学模型与壁面条件对升力体流场等离子体分布的影响，可得到以下两点认识：

（1）准确数值模拟高温流场等离子体分布，必须选择合适的化学反应模型。化学反应模型的选取主要取决于流场等离子体的温度和压力，另外还与流动的非平衡松驰过程有关。对于本文研究的升力体外形和飞行状态，头部流场温度最高达10000K以上，此时除了考虑NO＋离子，O＋离子对电子数密度的也有重要贡献，化学反应模型中应包括O＋。采用不同非平衡气体模型（热化学非平衡或化学非平衡气体模型）对高温流场等离子体分布计算结果也会产生一定影响。

（2）壁面条件特别是壁面催化条件对流场等离子体分布具有较大影响。为了准确数值模拟升力体飞行器绕流等离子体分布，需要结合飞行和地面试验测量结果，获得可靠的材料催化特性和壁面温度数据，在此基础上对其周围等离子体环境进行精细模拟。

［1］MUYLAERT J，CIPOLLINI F，et al.In－flight research on real gas effects using the ESA EXPERT vehicles［R］.AIAA 2003－6981，2003.

［2］MORABITO D，KORNFELD R，et al.The mars phoenix communications brownout during entry into the martian atmosphere［R］.IPN Progress Report，2009：42－179.

［3］HARTUNIAN R A，STEWART G E，et al.Implications and mitigation of radio frequency blackout during reentry of reusable launch vehicles［R］.AIAA 2007－6633，2007.

［4］钟育民，谌明，卢满宏，等.无线电波在再入等离子体中传输的衰减模型及仿真验证［J］.遥测遥控，2010，31（2）：1－6.

［5］MATHER D E，PASQUAL J M，SILLENCE J P.Radio frequency（RF）blackout during hypersonic reentry［R］.AIAA 2005－3443.2005.

［6］STARKEY R P.Electromagnetic wave／magnetoactive plasma sheath interaction for hypersonic vehicle telemetry blackout analysis［R］.AIAA 2003－4167，2003.

［7］SCALABRIN L C，BOYD I D.Numerical simulation of weakly ionized hypersonic flow for reentry configurations［R］.AIAA 2006－3773，2006.

［8］BLOTTNER F G.Viscous shock layer at the stagnation point with nonequilibrium air chemistry［J］.AIAA Journal，1969，7（12）：2281－2288.

［9］DUNN M G，KANG S W.Theoretical and experimental studies of reentry plasmas［R］.NASA－CR－2232，1973.

［10］GUPTA R N，YOS J M，et al.A review of reaction rates and thermodynamic and transport properties for an 11－species air model for chemical and thermal nonequilibrium calculations to 30000K［R］.NASA－RP－1232，1990.

［11］LIN T C，SPROUL L K，et al.Reentry plasma effects on electromagnetic wave propagation［R］.AIAA 95－1942.

［12］NUSCA M J，COOPER G R.Computational simulation of EM attenuation by plasmas formed in hypervelocity atmosphereic flight［R］.AIAA－98－1573，1998.

［13］董维中.热化学非平衡效应对高超声速流动影响的数值计算与分析［D］.［博士学位论文］.北京航空航天大学，1996.

［14］CANDLER G V，MacCORMACK R W.The computation of hypersonic ionized flows in chemical and thermal nonequilibium［R］.AIAA－88－051，1988.

［15］高铁锁，李椿萱，董维中，等.高超声速电离绕流的数值模拟［J］.空气动力学学报，2002，20（2）：184－191.

［16］YEE H C，KLOPFER G H，et al.High－resolution shock－capturing schemes for inviscid and viscous hypersonic flows［R］.NASA－TM－100097，1988.