不同尺度飞行器周围等离子体分布及电磁波传输效应

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

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

0 引言

高超声速技术是未来先进飞行器发展的核心技术之一。在对高超声速飞行器优化设计时,除了分析评估飞行器的气动特性和飞行性能,还需要评估飞行器周围等离子体流动对目标特性和通信性能的影响[1-3]。高超声速飞行器在大气中飞行时,与周围气体发生强烈作用产生强激波效应,使得波后气体温度急剧升高,气体分子发生振动激发、分解和电离等复杂气动物理化学现象,在飞行器周围形成高温电离气体层俗称等离子体鞘套。电磁波通过等离子鞘套传播时,被等离子体反射、折射和吸收,发生强度衰减、传播方向偏折和相位畸变等效应,导致电磁波作用距离缩短和信噪比下降等电磁性能退化现象,情况严重时,等离子体使电磁波传输中断,出现“黑障”现象,这些现象一般统称为气动电磁波传输效应[4]。气动电磁波传输效应特别是“黑障”现象严重影响飞行器与地面之间的微波通信性能,对飞行器实时控制和飞行安全构成严重威胁[5-6]。因此,研究气动电磁波传输效应对“黑障”现象的分析评估及高超声速飞行器通信系统的设计都具有重要意义。

随着临近空间飞行器飞行速度的不断提高,气动电磁波传输效应对电磁通信的影响问题愈加突出。考虑到飞行试验模拟的高成本及地面试验对实际飞行条件模拟能力限制等因素,对气动电磁波传输效应问题的研究目前仍以数值模拟为主要手段。1995年,Lin等[7]总结分析了飞行器等离子体鞘套的计算和试验研究状况,给出了均匀等离子体中各种电磁波传输效应的解析表述;1998年,Nusca等[8]采用7组分化学非平衡N-S方程和波动方程,计算分析了平面电磁波在无线电衰减测量C(Radio Attenuation Measurement C,RAM-C)再入等离子体鞘套中衰减效应;2002年,Funaki等[9]结合固体火箭发动机地面试验,预估分析了发动机羽流等离子体对不同波段电磁波的衰减特性;2003年,Starkey等[10]采用Park化学模型和解析方法,分析了航天飞机等多种飞行器等离子体鞘套对电磁波通信的影响;2006年,White等[11]基于高阶数值方法,分析了RAM-C再入等离子体鞘套的电磁波传输效应;2009年,Kim等[12]采用非平衡等离子体数值模拟方法及传输效应解析方法,研究了高超声速飞行器轨道再进入试验(Obiter Re-entry Experiment,OREX)绕流等离子体对电磁波传输的影响,探索外加磁场降低通信窗口等离子体传输效应的方法;2018年,龚旻等[1]对临近空间飞行器“黑障”现象的数值和地面模拟方法进行了回顾总结;2019年,左光等[13]采用数值及解析方法,分析了大钝头返回舱和类X-37B升力体飞行器再入等离子体鞘套对电磁波传输的影响特征。上述研究工作的共性特点是把流体力学和电磁学结合起来,分析高超声速飞行器等离子体鞘套中电磁波传输效应,但对不同特征几何尺度高超声速飞行器等离子体鞘套及电磁波传输效应的产生机制和规律分析较少。

本文基于求解非平衡流场N-S控制方程和求解电磁场波动方程的模拟方法,采用自主研发的气动物理流场数值模拟软件和再入黑障预测分析软件,研究高超声速再入体周围等离子体的产生机制及气动电磁波传输效应,重点分析不同频率电磁波在再入等离子体鞘套中传播的衰减效应,研究不同球头半径尺寸再入体的周围等离子体的分布规律及其对电磁波衰减的影响特征。

1 高超声速等离子体流场数值模拟研究

1.1 控制方程

高超声速飞行器从外层空间再入大气层过程中,由于热化学过程与流动过程的特征时间尺度效应,飞行器绕流等离子体一般要经历热力学和化学非平衡过程。文献[14]的研究发现,热力学非平衡效应对高超声速飞行器非平衡绕流中等离子体分布的影响很小,此时可以不单独考虑气体分子的振动能量模态,而是基于单温度气体模型,通过求解含化学反应源项的化学非平衡流动N-S方程,数值模拟高超声速行器周围等离子体鞘套,三维化学非平衡流动的N-S方程的无量纲化形式[14]如下：

(1)

式中：Q为守恒变量,Q=(ρi,ρ,ρu,ρv,ρw,ρE)T,ρi和ρ是组分的分密度和混合气体的总密度,u、v、w为直角坐标下3个方向的速度,E为混合气体的总能;Re是来流雷诺数;F、G、H和FV、GV、HV分别对应3个方向的无黏和黏性通量项;W为化学非平衡源项,W=(wi,0,0,0,0,0)T,wi为i组分生成源项,

(2)

(3)

1.2 化学模型

高超声速飞行器与来流空气发生作用形成脱体激波,波后压缩空气温度急剧增高,使得来流压缩空气不仅发生离解反应和置换反应,而且发生各种电离反应,在不计空气中微量组分的情况下,认为来流空气由O2和N2组成,此时高温空气中的电离反应主要包括缔合电离和碰撞电离反应机制[14]：

2 数值计算方法

2.1 等离子体流场数值模拟方法

基于有限差分方法对流动控制方程式(1)进行数值求解：对于方程中的无黏通量项,采用上下对称的高斯-赛德尔(Lower-Upper Symmetric Gauss Seidel,LU-SGS)隐式处理分法[18],以解决无黏通量雅可比矩阵直接求逆带来的运算量大的问题;为了解决流动和化学过程特征时间尺度效应导致的刚性问题,采用全隐式耦合方法对方程中化学反应源项进行处理[19],即用同一时间尺度同时求解流动方程和化学反应方程以保证计算的稳定性和收敛性;为了准确模拟流动中激波结构和边界层特性,同时保证计算的鲁棒性,采用高分辨压力权函数修正的迎风型矢通量分裂(Advection Upstream Splitting Method by Pressure-based Weight functions ,AUSMPW+)格式[20]离散无黏通量,采用中心差分格式离散项黏性通量项;为了模拟高空稀薄效应的影响,引入壁面参数滑移修正模型[21];为了模拟多组分混合气体输运特性,混合气体的黏性系数和热传导系数用Wilke半经验公式计算,各组分的输运系数基于Blotter曲线拟合公式和Eucken关系式计算,扩散系数采用等效二元扩散模型计算,具体计算方法详见文献[18-19]。

2.2 等离子体中电磁波传输特性预测方法

考虑一维情况,即等离子体参数仅在z轴方向非均匀分布,且平面电磁波沿z轴正向传播,此时电场只是随着z轴发生变化,设电场平行于y轴,此时波动方程[22-24]为

(4)

式中：Ey为电场强度;k为波数。对于缓变介质,即介质的电磁参数在z轴方向变化较小,该方程的Wentzel-Kramers-Brillouin(WKB)解[24]为

(5)

式中：Ey0为z=0 m处的电场强度。设电磁波从z=0 m处垂直入射到等离子体内部,并在z=d界面处透射出来,电磁波在此处的能量衰减为

(6)

在各向同性非磁化等离子体介质中波数k为复数,其表达式为

(7)

式中：ω与ωp分别表示电磁波与等离子体的角频率;c为光速;r为等离子体介质的相对介电常数;ν为等离子体的碰撞频率,碰撞频率基于工程经验关系进行计算[4]。

3 数值计算结果

3.1 不同天线位置等离子体分布及其电磁传输特性

研究状态为RAM-C再入飞行状态[25]：再入体为球锥对称体外形,再入高度H=71 km,飞行速度V=7.65 km/s,壁面催化条件为非催化壁面(Non Catalytic Wall,NCW),壁面温度Tw=1 500 K,飞行攻角近似为0°。在零攻角飞行状态下等离子体流动为轴对称流动,等离子体分布沿轴向变化是主要的,因此通信天线的轴向安装位置将严重影响电磁波的传输效果。

图1 不同轴向位置法向峰值与临界等离子体分布Fig.1 Peak and critical plasma parameters distribution at different axial positions

实际上,再入体周围等离子体鞘套处于空间(定常情况下)非均匀分布状态,绕流等离子体厚度和大小分布均在变化。为了定量分析再入体周围等离子体鞘套对电磁波衰减的影响,以再入等离子体参数为基础,认为电磁波沿物面法向传播,采用上述WKB方法预测电磁波通过等离子体的衰减量,如图2所示。从图2(a)可以看出,对于C频段和X频段高频电磁波,其在等离子体中的衰减沿轴向总体上变小,而低频VHF电磁波的衰减量在球头附近区域出现极小值和极大值,使得电磁波衰减沿轴向呈现非线性变化。尽管越靠近头部峰值电子数密度越高,但头部等离子体鞘套的厚度较小,且等离子体参数沿物面法向和流向非均匀分布,电磁波衰减随等离子体参数(电子数密度和碰撞频率)非线性变化,导致不同频率电磁波衰减沿轴向变化出现不同变化规律。总之,不管是低频(VHF频段)还是高频(C、X频段)电磁波,其通过飞行器后部等离子体鞘套的衰减量均小于通过头部等离子体的衰减量,而且在同一轴向位置处,电磁波频率越高其衰减量越小,因此把天线位置安装在靠近飞行器后身部,或适当选择频率较高的电磁频段,有利于减少飞行器等离子体鞘套对电磁通信的影响。

图2 不同轴向位置的电磁波能量衰减和峰值碰撞频率Fig.2 Electromagnetic energy attenuation and peak collision frequency at different axial positions

为了分析等离子体碰撞频率的影响,图2(a)中还给出了理想无碰撞等离子体(假定碰撞频率ν=0 Hz)的预测结果。从图2(a)可见,碰撞频率对 C频段和X频段电磁波的衰减的几乎没有影响,而对低频VHF电磁波在球头附近区域等离子体中的衰减有较大影响。图2(b)给出了物面法向等离子体剖面的峰值碰撞频率随轴向的变化情况。由图2(b)可知：C频段和X频段的电磁波频率显著高于等离子体的碰撞频率,而VHF频段电磁波频率低于再入体球头附近等离子体的碰撞频率;碰撞频率对电磁波衰减的影响主要与碰撞频率与电磁波频率的相对大小有关[22],即当电磁波频率接近或小于等离子体碰撞频率时,碰撞频率对电磁波衰减才会表现出来,二者差异越大,等离子体碰撞频率对电磁波衰减的影响越大。

3.2 不同尺寸飞行器等离子体分布及电磁传输效应

以球头半径Rn=15.24 cm的RAM-C再入等离子体鞘套[25]为基础,改变球头半径尺寸,保持长度不变(L=1.295 m),NCW条件,壁面温度Tw=1 500 K,飞行攻角近似为0°。分析不同球头半径(Rn为15.24 cm、7.5 cm、2.5 cm)对再入体周围等离子体分布的影响。

图3给出不同球头半径飞行器等离子体分布,图中每条等值线上的数字表示对应的电子数密度。可见随着球头半径减小,等离子体鞘套层的厚度明显减小,鞘套层的外边界越贴近壁面。图4进一步给出不同尺寸再入体头部驻点线的流场温度和电子数密度分布情况,图中R为壁面法向距离。随着球头半径增大,沿驻点线激波层的厚度也随着增大,激波电离层电子数密度的峰值也随之增大,如图4(b)所示。大尺寸球头激波层气体电离度的增强意味着较强电离效应,而气体离解和电离等吸热反应使得大尺寸球头来流气体过激波后沿驻点线温度下降得更快,如图4(a)所示。实际上,在同一飞行条件下,飞行器头部尺寸越大,脱体激波面的曲率半径越大。曲面激波可看作是多个微元段平面斜激波段的组合,曲率半径越大,对应每个微元面的激波角越大,头部接近正激波面元段就越多,激波对来流气体的压缩效应就越强,波后压缩气体温升使得化学电离效应就越强。此时高速来流气体通过激波滞止过程中动能减少,一部分动能转化为气体粒子的热运动能,另一部则转化为离解和电离产物的零点能。

图3 不同球头半径再入体周围等离子体鞘套特征Fig.3 Plasma sheaths around the reentry body with different sphere radii

图4 不同球头半径再入体壁面法向等离子体分布(x=0 m)Fig.4 Wall-normal plasma distribution (x=0 m)

图5给出再入体身部1.234 m处(通信天线附近)壁面法向剖面的NO+质量分数和电子数密度分布,比较图4(b)和图5(b)不难看出,从头部驻点区域到尾部天线附近1.234 m处,对于不同尺寸的再入体,等离子体鞘套剖面内平均电子数密度降低均超过两个量级,而且头部半径越小,等离子体中电子数密度衰减越快,这是因为沿再入体长度方向等离子体特征电子数密度与球头半径的幂次方呈正比关系[6,22],球头半径越大,从头部开始等离子体层厚度相对于上游增加的速度越快,等离子体电子数密度相对于上游衰减越慢。从图5还可以发现,在天线附近1.234 m处,等离子体剖面内平均电子数密度和NO+质量分数随头部半径的变化规律基本一致,二者均随着头部半径的减小而减小,表明此时NO的电离效应对等离子体电子数密度分布起主导作用,这可以从图6中各种离子组分的分布特征得到进一步印证。从图6可以看出,对于不同头部半径尺寸的再入体,在此尾部天线附近位置处,沿物面法向的所有离子组分分布中,NO+质量分数最大,是最主要的电离组分,说明有关NO电离反应机制对等离子体中电子数密度分布的贡献最大。

图5 再入体壁面法向等离子体分布(x=1.234 m)Fig.5 Wall-normal plasma distribution (x=1.234 m)

图6 再入体壁面法向离子组分质量分数分布(x=1.234 m)Fig.6 Ion mass fraction along wall normal (x=1.234 m)

图7给出球头半径对等离子体中电磁波传输特性影响,可见,随着球头半径增加,由于峰值电子数密度随之增加,等离子体对电磁波传输过程中产生的衰减量也随之增加。实际上,随着球头半径增加,等离子体的厚度也随之增加,即球头半径对等离子体参数大小与厚度分布均产生影响,从而影响电磁波在其中的传输特性。从图7中还可以看出,再入体的球头半径改变时,等离子体对不同频率电磁波衰减的影响的程度有所不同,但电磁波衰减量沿轴向的变化规律基本一致。

图7 球头半径对电磁波衰减特性的影响Fig.7 Effect of sphere radius on energy attenuation of electromagnetic wave

3.3 OREX通信中断预测分析

在对飞行器OREX通信中断预测分析之前,先对本文计算模型方法及软件的可靠性进行验证分析。沿飞行弹道的高超声速飞行器周围等离子体及通信中断数据可以通过飞行测量获得,但目前公开发表飞行测量数据很少,而地面试验还不能完全模拟实际飞行条件,且缺乏天地换算相似准则。RAM-C飞行器再入试验提供了相对全面包括等离子体分布和通信中断的测量数据[12,25],该数据是目前少有的能够同时验证等离子体流场及其对电磁波通信影响预测方法的测量数据,这里基于此测量数据及弹道条件开展验证分析。数值模拟等离子体流场采用7组分Park化学反应模型[16],壁面采用完全催化壁面(Full Catalytic Wall,FCW),壁面温度为1 500 K,通信天线位置在x=6.4Rn处,飞行器以零攻角再入。图8(a)给出典型再入高度条件下剖面峰值电子数密度沿轴向分布,可见数值模拟结果与飞行测量结果符合很好。在沿弹道对RAM-C飞行器再入流场等离子体分布数值模拟的基础上,计算了对应飞行弹道条件下电磁波通过天线附近等离子体的能量衰减,计算中认为电磁波沿天线壁面位置的法向传播。一般认为电磁波衰减30 dB以上时发生通信中断[23],这样就可以通过预测的衰减量判断是否发生通信中断,获得沿弹道的通信中断区间,如图8(b)所示,其中测量数据包括VHF波段通信中断区间,中断时间持续约半分钟[6],而飞行测量只获得C与X波段通信中断区间,因此图中只比较了这两个波段通信中断的起始高度。从图(8)可见,RAM C-II沿弹道再入过程中,VHF、C和X波段电磁波通信中断预测与测量数据具有较好的一致性,预测的通信中断起始或结束高度最大误差在5 km以内。

图8 RAM-C飞行器再入预测结果与飞行测量对比Fig.8 Comparison between predictional and experimental data for RAM-C vehicle reentry

下面针对OREX条件[12,26]开展计算分析,其外形尺寸、等离子体探针及通信天线布置情况见图9。在OREX飞行器以零攻角再入[12]飞行条件下,飞行器绕流对称,探针附近等离子体分布可以反映通信天线附近的等离子体特征。图10给出了FCW和NCW条件下电子数密度的数值模拟结果与电子探针测量结果的比较,其中横坐标表示探针分布方向,离散点表示飞行器肩部分布探针的测量结果。由图10可以看出,数值模拟结果在总体上和飞行测量结果具有较好的一致性,高度84.01 km的符合程度更好,高度79.9 km的计算结果相对测量结果偏高,但计算误差在3倍左右,这与之前对RAM-C飞行器再入等离子体预测精度相当[14],表明了数值模拟等离子体模型算法及软件的可行性。至于高度79.9 km的数值模拟结果偏高的原因,很可能是高温空气反应特别是电离反应速率常数数据的不确定性所致,特别是高温条件下,不同化学动力学模型对等离子体分布的数值模拟结果可存在量级上差异[27],而高度79.9 km下流场中气体电离效应相对84.01 km更强,有关电离反应化学动力学数据对数值结果的影响更大。图11给出OREX飞行器在高度59.6 km条件下绕流等离子体及头部驻点线上离子组分的分布情况。从图11(a)可见,飞行器头部激波层内等离子体效应很强,从头部到飞行器接近尾部的肩部位置,等离子体浓度沿流向虽然有所减弱,但衰减速率比较缓慢,这与RAM-C飞行器的等离子体沿流向的变化特征存在明显不同(参见图1和图3),RAM-C飞行器周围等离子体浓度沿流向的衰减效应相比OREX飞行器更强,这是由于二者在飞行条件接近的情况下,OREX飞行器属于大钝头体外形,其球头半径和锥角都明显大于RAM-C飞行器球头半径和锥角,从3.2节的分析可知,飞行器头部球头尺寸及锥角的大小造成了飞行器周围等离子体分布及沿流向衰减程度的不同。从图11(b)中电离组分的分布可见,头部激波层等离子体流场中,NO的电离机制是最主要的,NO+离子对等离子体浓度的贡献最大,其次是O的电离效应,其他组分电离效应较弱。图12给出了OREX飞行器天线附近不同高度和不同催化壁条件下的壁面法向等离子体的分布情况,图中标识H84.0V7416表示再入高度84.0 km与对应再入速度7 416 m/s,其他类同。从图12不难看出,FCW和NCW的计算结果只是在靠近壁面附近有所差别,在离开壁面的大部分流场区域二者非常一致,沿再入弹道发展过程中,等离子体强度在高度 67 km 左右达到最高。图13是以此等离子体流场数据为基础获得的不同频段电磁波能量衰减沿再入弹道高度的变化情况,计算时认为电磁波的传播方向沿壁面安装天线位置的法线方向,不难发现,沿再入弹道等离子体对电磁波衰减沿再入高度降低逐渐增大并在高度67 km左右达到峰值然后开始减小。这是因为对于特定飞行器外形和飞行姿态,来流大气密度和再入速度对飞行器周围等离子体的分布其主导作用[12,22],而等离子体强度的沿高度分布规律决定了通信电磁波的衰减规律,在电磁波衰减出现峰值以前,再入速度沿高度降低减速较慢,大气密度沿高度降低增大,对等离子体分布和衰减起主导作用,而在电磁波衰减出现峰值以后,再入飞行速度沿高度快速降低成为等离子体衰减的主导因素。从图13还可以看出,FCW与NCW的衰减计算值相差不大,在认为电磁衰减30 dB以上发生通信中断情况下[22-23],X波段电磁波在高度84 km左右开始出现中断现象,而L和C波段出现中断的高度更高,由于再入段末端飞行减速很快,3个波段电磁波在高度50 km左右时衰减迅速降低而恢复通信,可见随着通信电磁波频率提高,等离子体对电磁波衰减及通信中断的影响减弱,发生通信中断的高度区间缩小。

图9 OREX飞行器外形尺寸、静电探针及天线位置Fig.9 Geometry with electrostatic probes and antenna of OREX vehicle

图10 OREX等离子体数值模拟与飞行测量值的比较Fig.10 Comparison between computational and experimental data of OREX

图11 OREX等离子体组分质量分数(H=59.6 km)Fig.11 Plasma and ion mass fraction for OREX (H=59.6 km)

图12 天线附近壁面法向电子数密度分布Fig.12 Electron number density near antenna

图13 OREX等离子体对电磁波衰减影响Fig.13 Electromagnetic wave attenuation due to OREX

4 结论

本文针对RAM-C和OREX钝头体高超声速试验飞行器,开展不同特征尺度下飞行器周围等离子体分布特性及气动电磁波传输效应的数值模拟和分析。得出以下主要结论：

1)在同一再入条件下,对于钝头体飞行器,保持长度不变而改变其球头尺寸,随着球头尺寸的增加,球头脱体激波内压缩气体的化学效应更强,来流气体动能更多地转化为波后气体离解和电离的化学能,导致飞行器周围等离子体鞘套厚度、等离子体的浓度以及对传播电磁波的衰减均随之增大,NO电离反应对飞行器天线附近等离子体分布特性起主导作用。

2)再入飞行器等离子体鞘套空间分布决定了对电磁波传输效应的不同特征：飞行器头部等离子体对电磁波衰减效应更强,随着流动向下游发展等离子体效应减弱;而在同样等离子体分布条件下,随着电磁波频率的提高,等离子体对电磁波的衰减越弱;可以通过合理选择天线安装位置和适当提高通信电磁波频率的方法减缓或消除等离子体鞘套对微波通信的影响,天线安装最好选择飞行器的后身部位。

3)大钝头高超声速飞行器(OREX飞行器)周围等离子体鞘套厚度大且电离度强,沿流动方向衰减慢,在沿再入弹道飞行过程中,随着飞行高度降低,等离子体鞘套对电磁波衰减逐渐增强出现峰值,然后随之逐渐减弱;随着通信电磁波工作频率提高,等离子体鞘套对电磁波通信中断的影响区间缩小。

4)针对自主发展的计算模型方法和软件开展验证分析,典型飞行状态下等离子体分布的数值模拟结果与飞行测量数据符合较好;典型飞行器沿弹道再入过程中,通信中断预测与测量结果较为一致。该计算分析工具可为高超声速飞行器的通信系统设计提供技术支持。