临近空间高超声速巡航飞行器红外特征

2022-05-13石安华李海燕石卫波梁世昌

兵工学报 2022年4期

石安华，李海燕，石卫波，梁世昌

(中国空气动力研究与发展中心超高速空气动力研究所，四川绵阳 621000)

0 引言

临近空间高超声速飞行器近年来受到各军事强国的热捧，在临近空间高超声速武器的攻防对抗研究中，目标辐射特性是重要的研究内容，除了试验研究外，建模研究作为主要方法之一也得到大力发展。临近空间高超声速飞行器可分为滑翔类和带动力类，滑翔类高超声速飞行器辐射特性主要受气动热和高温绕流气体影响，而带动力类高超声速飞行器辐射特性主要受气动热和发动机喷焰辐射影响。针对带动力类高超声速飞行器辐射特性，尽管已有一些学者对此开展了相关研究工作，但鲜有同时考虑时变气动加热、壳体结构与热传导、喷焰中气体分子能级跃迁辐射对带动力类高超声速飞行器辐射特征影响的研究，使得分析结果较为粗糙。

本文针对上述问题，从考虑流动过程化学反应、燃烧化学反应、气动热、热传导、本体辐射、分子跃迁辐射等辐射产生过程及机制出发，开展以超燃冲压发动机为动力的类X-51A(见图1)高超声速巡航飞行器在典型飞行状态下的红外辐射特征分析研究，分析结果可为临近空间高超声速导弹的相关探测技术研究提供数据参考。

图1 类X-51A高超声速巡航飞行器数值模型Fig.1 Digital model for hypersonic cruise vehicle like X-51A

1 计算方法

1.1 流场及气动热特性计算方法

流场包括飞行器绕流场和绕流场与吸气式超燃冲压发动机燃气耦合的喷焰流场。

临近空间高超声速巡航飞行器所在飞行区域的流动是连续流动，飞行器燃料通常采用航空煤油，由于处于高超声速飞行条件，气体温度高、化学反应强烈，绕流和喷焰流场中考虑O、O、N、NO、OH、NH、H、CO、NO、NO、HO、HO、CO、HO、CH、H和N 17种组分，化学动力学模型使用包括46种反应的燃烧化学动力学模型。假设飞行器外部绕流和超燃冲压发动机内部的流动是可压缩湍流流动。对于飞行器外部绕流、超燃发动机尾喷流与飞行器外部绕流相互作用下的流场特性模拟，采用三维可压缩化学反应湍流的Favre平均方程计算流体动力学方法模拟。在笛卡尔坐标系下，三维守恒形式的Favre平均方程为

(1)

式中：为守恒变量向量；为时间(s)；、和分别为笛卡尔坐标系中轴、轴和轴方向的无黏通量向量；、和分别为笛卡尔坐标系中轴、轴和轴方向的黏性通量向量；为化学反应源项向量。

为封闭平均方程，引入Boussinesq涡黏性假设，并采用-两方程模型来模拟湍流效应，为湍动能(J/kg)，为湍动能耗散率(1s)。首先对三维化学反应湍流控制方程进行无量纲化，然后采用有限体积法对无量纲化后的控制方程进行离散。对离散后的无量纲控制方程组，采用上下对称Gauss-Seidel (LU-SGS)隐式方法进行迭代求解。气动热采用表面热流公式计算分析：

(2)

式中：为壁面热流(W/m)；为热传导系数(W/(m·K))；为温度(K)；为飞行器表面局部法向距离(m)；为化学组分种类数量；为当地气体密度(kg/m)；和分别为第种组分的扩散系数(m/s)和焓值(J/kg)；为第种组分的质量分数。

1.2 燃烧室流动特性计算方法

超燃冲压发动机尾喷管前的燃烧室流场采用准一维流动模型，将准一维Euler方程与有限速率化学动力学方程耦合，在方程中考虑壁面摩擦、面积变化和燃料质量添加的影响。然后对无量纲化的准一维Euler方程应用有限体积法进行离散，采用隐式方法求解离散后的方程。

在笛卡尔坐标系下，准一维Euler方程和组分连续方程的强耦合形式为

(3)

式中：为燃烧室截面积(m)。

在准一维Euler方程和组分连续方程中，质量添加项出现在方程的右端源项中，此外，右端源项也考虑了壁面摩擦、面积变化以及化学反应的贡献。燃料采用以CH为简化分子式的JP7航空煤油，其燃烧化学动力学模型使用CH燃烧的9种组分(H、H、O、O、OH、HO、CO、CO和CH)、13个化学反应方程模型，高温条件下空气与煤油反应混合气体中各组分热力学参数由Janaff表获得。CH的当量油气比为，化学恰当比为，根据油气比关系求得喷射的煤油流量：

(4)

考虑到燃料的混合过程，质量添加项采用指数分布进行处理。采用参考温度法来计算摩擦系数，通过雷诺比拟可以计算传热系数。对于燃烧室前方的隔离段而言，准一维方法无法准确捕捉到其中的预燃激波串这种复杂非一维结构。对于激波串的处理，采用工程预燃激波串模型进行计算。将预燃激波串模型耦合到准一维非稳态方法中，得到一种改进的准一维非稳态计算方法，结合相应的碳氢燃料化学反应动力学模型和燃料混合模型，来获得燃烧室出口的参数。

1.3 流场辐射特性计算方法

流场组分光谱辐射使用分子振转光谱模型，以振转谱线标准吸收强度为基础，根据爱因斯坦辐射与吸收关系计算分子光谱辐射强度及吸收强度。

辐射传输使用窄带吸收模型，对传输路径上的辐射按组分种类分别计算再求和，计算方向从观测点起至视线到达飞行器表面或穿出视场。

流场辐射使用光学薄传输模型计算会带来较大误差，在流场光辐射强度计算中使用吸收模型，在无粉尘颗粒的气体辐射与传输环境中忽略光散射影响。假设流场中第个体积元d内存在的化学组分共有种，d中第种气体组分的数密度为,(个m)、辐射强度为,、吸收系数为,，则

(5)

(6)

式中：、分别为流场体积元d的辐射强度(W/Sr)、吸收系数(1m)。

在光谱分辨率要求不高的情况下，分子的光谱吸收系数通常采用近似算法，在本文中采用窄带平均吸收系数模型。窄带平均吸收系数模型是从吸收系数的严格公式出发，计算出每条单独谱线的吸收系数后，将研究的波段划分为多个窄带Δ，计算落于此Δ内所有谱线的平均吸收系数，采用(7)式计算：

(7)

体积元d照射到观测点处的辐照度为

(8)

(8)式中的求和项是针对辐射传输路径上所有体积元d进行的。

1.4 飞行器本体辐射特性计算方法

141 表面温度求解方法

飞行器表面温度是关系到飞行器辐射计算的关键物理量，其受到热环境、飞行器结构、材料物性参数等多种因素的影响。为精细描述包含表面温度在内的飞行器结构温度场变化，一般需要求解三维瞬态热传导方程。本文采用有限元方法对热传导过程进行模拟，求得飞行器表面温度。

瞬态温度场的温度在直角坐标中应满足的微分方程为

(9)

式中：为飞行器壳体材料密度(kg/m)；为材料比热容(J/(kg·K))；、、分别为材料沿轴、轴、轴方向的热导率(W/(m·K))；为物体内部热源的热流密度(W/m)。

求解域内温度场采用的边界条件为

(10)

式中：、、为边界外法线方向余弦；为边界上给定热流量或热流密度(W/m)。

求解瞬态温度场问题是求解在一定初始条件下满足瞬态热传导方程及边界条件的温度场，为空间坐标与时间的函数。采用四节点四面体单元对空间进行离散，依据泛函理论，将传热控制方程离散为代数方程组，写成矩阵形式为

(11)

142 本体辐射特性计算方法

假定飞行器是朗伯辐射体，半球发射率、方向发射率和光谱发射率相等，且与方向无关，则本体红外辐射可采用灰体辐射模型计算，表面网格使用温度场网格的外表面部分。

利用普朗克公式，对波长从到积分，得到黑体面源在～波段范围内的辐射出射度为

(12)

式中：为黑体辐射出射度(W/m)；()为黑体的光谱辐射出射度(W/m)；为第一辐射常数，=3741 8×10W·m；为第二辐射常数，=1438 8×10m·K。

飞行器表面材料发射率为、温度为的小面源的辐射出射度=·，光谱辐射亮度()=()π。

由于投影到观察方向的面源辐射强度具有独立性，相互之间不干扰，并具有可叠加性；分析中考虑飞行器离观察点很远，面元中心至观察点连线的矢量几乎相互平行，对飞行器表面各点辐射积分即可得到飞行器本体的总辐射。

1.5 总辐射计算方法

考虑本体辐射穿过流场的吸收效应、本体对部分流场辐射的遮挡效应以及本体之间的遮挡效应，即可计算得到目标总辐射。

2 计算结果

2.1 计算条件

计算中使用的飞行器外形见图1，飞行器壳体各部位使用的材料和厚度取值见表1，壳体材料物性参数见表2，计算飞行状态为飞行器以马赫数=6、攻角=8°在海拔为30 km、40 km两种高度下长时间飞行150 s、400 s和800 s共6个状态。计算中，燃料使用JP7航空煤油(简化分子式为CH)，油气当量比假设为0.8；飞行器初始温度假定为298 K。

表1 飞行器壳体材料及厚度

表2 飞行器壳体材料物性参数

本文光谱辐射强度计算光谱范围为1～14 μm，光辐射强度及光辐射亮度计算波段为3～5 μm、8～12 μm。计算辐射方向为从头部开始往背风面方向绕到尾部共17个方向(=10°～170°，间隔10°)。

2.2 典型计算结果

图2给出了飞行器在30 km高度飞行时的典型流场参数和飞行150 s时的本体温度。

图2 典型飞行状态下飞行器流场参数、表面热流及温度(H=30 km、Ma=6、α=8°)Fig.2 Flow field parameters, heat flux and surface temperature (H=30 km, Ma=6, α=8°)

图3、图4给出了飞行器在30 km高度飞行150 s、400 s、800 s时，飞行器在90°辐射方向上波段分别为3～5 μm、8～12 μm的辐射亮度。从图3和图4中可以看出，飞行器喷焰辐射亮度随飞行时间基本无变化，而本体辐射亮度随飞行时间增加而增强。

图3 飞行器在典型飞行时刻在3～5 μm波段的辐射亮度(H=30 km, Ma=6, α=8°, θ =90°)Fig.3 Luminance in 3～5 μm wave band at typical flight time (H=30 km, Ma=6, α=8°,θ =90°)

图4 飞行器在典型飞行时刻在8～12 μm波段的辐射亮度(H=30 km, Ma=6, α=8°, θ =90°)Fig.4 Luminance in 8～12 μm wave band at typical flight time (H=30 km, Ma=6, α=8°, θ=90°)

图5、图6给出了飞行器在30 km高度飞行 800 s 时，飞行器在50°和130°辐射方向上波段分别为3～5 μm、8～12 μm的辐射亮度。从图5和图6中可以看出，不同辐射方向上的辐射图像是不一样的，50°方向上的辐射亮度强于130°方向上的辐射亮度。

图5 飞行器飞行800 s时典型辐射方向3～5 μm波段的辐射亮度(H=30 km、Ma=6、α=8°)Fig.5 Luminance in 3～5 μm wave band at 800 s and different θ (H=30 km, Ma=6, α=8°)

图6 飞行器飞行800 s时典型辐射方向8～12 μm波段的辐射亮度(H=30 km、Ma=6、α=8°)Fig.6 Luminance in 8～12 μm wave band at 800 s and different θ (H=30 km, Ma=6, α=8°)

图7给出了飞行器在30 km高度飞行800 s时飞行器典型辐射方向的光谱辐射强度随波段的变化曲线。从图7中可以看出：飞行器辐射总体体现为在本体热辐射上叠加喷焰气体的特征辐射，主要为HO、CO、CO气体的特征辐射，特征辐射为峰值波长在1.52 μm、2.68 μm、4.39 μm处的3个辐射带。其中，1.52 μm辐射带主要为HO辐射，2.68 μm辐射带主要为CO、HO、CO辐射，4.39 μm辐射带主要为CO辐射；特征谱波长不随辐射方向变化，但特征谱辐射强度随辐射方向变化剧烈，在90°左右方向上特征谱辐射明显强于本体辐射，这是因为这些方向上喷焰辐射投影面积大，强辐射区穿过的喷焰气体距离短。此外，在该时刻的本体红外辐射主要集中在短波红外波段，主要是因为此时飞行器本体大部分区域温度已到上千度。而在飞行初始阶段，本体红外辐射应主要集中在中长波红外波段，这是因为飞行器在初始飞行阶段表面温度还不高造成的。

图7 飞行器在典型辐射方向上的光谱辐射强度(H=30 km、Ma=6、α=8°、t=800 s)Fig.7 Spectral intensity at different θ (H=30 km,Ma=6, α=8°, t=800 s)

图8所示为飞行器在30 km高度飞行800 s时典型特征光谱(波长1.52 μm、2.68 μm、4.39 μm)的辐射强度与波长1.6 μm处的光谱辐射强度(基本为本体光谱辐射强度)随辐射方向的变化曲线。从图8中可以看出，在特征谱辐射波长位置，喷焰气体特征谱辐射明显强于本体辐射，如果用特征辐射波长做窄带测量，特别是在飞行器正上方附近方向，将具有最高的信噪比。

图8 飞行器在不同方向的特征光谱辐射强度(H=30 km、Ma=6、α=8°、t=800 s)Fig.8 Typical spectral intensity at different θ (H=30 km,Ma=6, α=8°, t=800 s)

飞行器在40 km高度飞行时的辐射出射度与光谱辐射强度与在30 km高度飞行时类似，只是相同时刻辐射出射度和光谱辐射强度变小。表3给出了飞行器在70°辐射方向上、两个飞行高度下在3～5 μm和8～12 μm波段本体辐射强度和总辐射强度。从表3中可以看出：飞行器喷焰辐射随飞行时间变化不大，而本体辐射随飞行时间增加而增加；3～5 μm波段的辐射远强于8～12 μm波段的辐射；在3～5 μm和8～12 μm波段，喷焰辐射占比在13%～17%左右，与在特征辐射波长的窄带内辐射相比差异十分明显，窄带特征辐射波长内喷焰辐射占比可达50%～90%。