1 引言
目前，以超燃冲压发动机为特征的吸气式高超声速飞行器技术已成为世界各军事强国的重点发展项目。以X51-A为代表的乘波体飞行器，将超燃冲压发动机与飞行器机身相耦合设计，使其可利用飞行器前体作为外部压缩面为进气道提供压缩空气，并利用后体作为喷管型面产生附加推力[1]。自2010年5月26日美国波音公司“驭波者”的第一次试飞成功以来[2]，因其快速高效的全球打击潜力，军方和学者对X-51A高超声速飞行器的发展和关注度日趋提高。在飞行器目标特性捕获、识别和预(告)警技术以及攻防对抗中，对目标光辐射信号的研究显得尤为重要[3]，有必要开展相关的数值模拟和分析研究。
高超声速飞行器通常会以超过5马赫的速度在大气层内飞行，由于对周围空气强烈的压缩和摩擦作用，会在头部周围形成高温激波层，该激波层内伴随着空气组分的离解、电离、复合和光化学反应等复杂的物理化学和辐射跃迁过程，从而产生气体光辐射效应。同时，高速飞行时飞行器表面受到气动加热作用，壁面温度会升高，且典型部件的温度差异明显，目标的红外辐射信号也会受到显著影响。此外，冲压发动机工作工程中，高温燃烧产物(完全及未完全燃烧产物)与周围大气剧烈相互作用，产生湍流卷吸、大气掺混以及复燃等作用，喷焰中的多原子气体(H2O、CO2、CO)及颗粒物(炭黑)会产生强烈的红外辐射，此过程涉及复杂的物理化学过程。然而，在被动巡航阶段，对类X-51A飞行器而言，目标辐射信号主要来源于飞行器激波绕流场气体和本体表面部分。
由于相关X51-A的公开文献较少，国内的初探类研究均参考赖特-帕特森空军基地气动所Hank[4]等的公开文献，文中给出了该飞行器的概略尺寸、典型飞行高度、飞行时间和状态以及推进剂的种类等基本信息。中国科学技术大学的程迪[5]等以及中国空气动力研究与发展中心的周正[6]等分析总结了大量资料数据，建立了一套可行的外形反设计重构方法，实现了X-51A飞行器外形的反向建模和外形确认。北京航空航天大学的Chen J[7]等以X-51A为对象，对“乘波体”高超声速飞行器的气动特性进行了研究。
针对类X-51A高超声速飞行器辐射特性研究的需求，本文将围绕飞行器壁面温度预估和红外辐射特性建模工作开展，包含目标乘波体的外形、结构和冲压发动机几何参数，飞行器关机巡航状态下的流场参数和壁面温度预测，红外辐射特性计算模型的建立，以及飞行器在典型工况下不同探测方向的红外辐射特性计算分析。
2 计算模型
2.1 流动控制方程
在连续流假设下，描述三维流动的N-S方程有如下的形式：
(1)
式中: ρ、ui、p分别表示气体的密度、速度和压强；总能量表示为
黏性张量表示为
热流项表示为
采用雷诺平均方法对上述方程进行求解计算，采用剪切应力输运(SST) k-ω湍流模型对上述方程进行封闭处理。计算方程组采用有限体积法进行离散，黏性通量项采用二阶中心差分格式，无黏通量采用二阶迎风格式，并采用时间推进法获得计算模型的定常结果，详细的计算过程参考文献[8]。
2.2 壁面温度特性计算模型
高温绕流场对飞行器壁面的加热热流主要由对流热流、组分扩散热流和气体辐射热流3部分组成。在不考虑辐射加热作用下，壁面气动热流有如下形式：
(2)
其中，n表示壁面的法向方向坐标。等式右边第一项为对流热流项，第二项为组分扩散热流项。
本文忽略了烧蚀效应和壁面催化效应，因此在垂直壁面方向上的热流平衡方程为：
(3)
其中：ε为表面材料发射率，本文取ε=0.8； δ为斯蒂芬-玻耳兹曼常数，其值为5.67×10-8 W/(m2·K4)。qc表示以热传导形式从壁面到机体结构内部的热流。在本文中壁面采用辐射平衡边界，即：
(4)
2.3 红外辐射特性计算模型
在流场参数的基础上，建立流场辐射传输计算模型。辐射传输示意图见图1。由辐射传输理论，考虑介质内的吸收、发射、散射时辐射传输方程为[9]：
(5)
其中：Ιλ(s)和Ιbλ(s)表示光谱辐射强度和黑体光谱辐射强度；κaλ和κsλ分别为介质的吸收系数和散射系数。
在不考虑粒子散射效应下，引入光学厚度τλ = κeλs，辐射传输方程简化为：
(6)
上述辐射传输方程可采用视线踪迹法(LOS)[10]求解，其基本思想是将射线在流场中的传输简化为一维多层介质的辐射传输问题，即沿探测方向平行线的LOS与流场相交时，LOS会一直延续下去，直到遇到飞行器表面或离开流场。如果将射线沿着L方向通过流场的路径分为N层，每一层的介质认为是均匀、等温的，则考虑每层介质的吸收、发射，逐层递推最终可得到LOS射线的辐射出射强度。如图1所示，沿L方向的逐层递推公式可表达为：
图1 LOS辐射传输模型示意图
Fig.1 Radiative transfer model of LOS method
(7)
式中：
为第i层的方向光谱辐射强度；
为第i层的黑体光谱辐射强度；
为第i层的介质透过率。
通过以上方法可求得以任意方向穿过流场的射线的光谱辐射强度，对于在λ1～λ2谱段内的辐射强度有
IΔλ=
Iλdλ
(8)
3 计算条件
3.1 几何模型
本文以类X-51A气动外形的飞行器为研究对象，利用反向建模技术[5]获得目标乘波体的外形结构尺寸，冲压发动机几何参数参考文献[6]中的数据。图2为X-51A飞行器的外形结构[11]，其主要建模参数在表1中给出。
图2 类X-51A飞行器巡航体外形结构示意图
Fig.2 Configuration of cruising body of X-51A type waverider
表1 类X-51飞行器主要建模参数
Table 1 Main modeling parameters of X-51 typewaverider
主要参数参考值进气道参数参考值机体总长/mm4 267喉道高度/mm44最大宽度/mm580内流道宽/mm225最大高度/mm566总收缩比4.9翼展/mm927内收缩比2.4
3.2 流场计算域和边界条件
依据模型的对称性，在不考虑侧滑角情况下，为降低计算量，选取1/2计算域进行网格划分，网格数约为224万。模型边界层处区域采用六面体网格，远壁面区为四面体非结构网格。为了保证模型表面摩阻的计算误差和激波层的捕获，依据文献[12]给出的边界层网格建议，本模型的法向壁面网格选取为y+≈2。图3所示为部分计算域及模型网格分布情况。
图3 计算域及模型局部网格示意图
Fig.3 Computational domain and local model surface grid
3.3 流场计算条件与方法
本文利用ANSYS-FLUENT商业软件计算飞行器的流场分布，采用了SST方程k-ω湍流模型[13]，AUSM+格式离散对流项和中心差分离散黏性项的数值格式，化学反应模型含空气5组分(O、N、NO、O2和N2)的有限速率化学反应模型[14]，壁面温度通过辐射平衡条件来表征[15]，即高温气体对壁面的气动热流和壁面向外的辐射热流保持平衡。边界条件类型包含压力远场、压力出口以及辐射无滑移壁面边界。
依据文献公开的飞行时序[4]，X-51A的典型巡航设计飞行条件为：飞行马赫数为6.01，飞行高度21.3 km。该飞行高度下标准大气的静压为3 981 Pa，静温为217 K。
4 计算结果与分析
4.1 流场参数分布
类X-51A乘波体在以近6马赫数的飞行速度关机巡航状态下，流场参数发生复杂的波系交互作用。在图4中，表示了计算域对称面上的马赫数分布。可以看出，在飞行器周围呈现出较高的马赫数，在内流道内马赫数较低，且在平直喉道段保持马赫数为1，在喷管扩张段马赫数骤然升高，且呈现出两型马赫波节。这是由于在经由平直通道的气体经历喷管扩张段和喷口与大气2个膨胀过程引起的。图5为对称面上的温度分布云图，仿真结果表明其绕流场及尾迹流场内温度水平低于2 000 K的水平，高温区域主要集中在飞行器头部、鸭嘴型面下部以及进气道前段。在这一温度水平下，空气中氮氧成分尚未达到解离的温度水平，可认为绕流场尚无NO等红外强发射组分产生。
图4 对称面上的马赫数分布云图
Fig.4 Contour of Mach number on symmetry plane
图5 对称面上的流场温度分布云图
Fig.5 Contour of temperature distribution of flow field
on symmetry plane
4.2 壁面温度分布
飞行器壁面高温部件会发生强烈的红外光谱辐射。对壁面而言，其发射光谱符合普朗克黑体定律，发射连续的光谱辐射。在图6中，给出了类X-51A飞行器的表面温度分布。可以看出，高温区域主要集中在飞行器头部、机身前体侧部和翼缘部分，最高温度高达1 430 K。机身大部分区域的温度在600～800 K的水平。
图6 蒙皮温度分布云图
Fig.6 Contour of skin temperature distribution
4.3 红外辐射亮度
为了探究不同观测角度下飞行器的本征红外辐射特征，本文考虑了4种观测角度，分别为俯视、正侧视、仰视和前视观测。图7分别表示了4种观测角度下类X51-A飞行器红外波段(1～14 μm)的光谱辐射亮度分布。图7中可以看出，不同观测角度下，红外辐射亮度分布差异明显，其差异主要来自机身蒙皮温度的大小和分布不同。由于高温区域主要集中在机身头部、进气道前段以及尾翼部分，因此红外辐射亮度较高的区域也集中该高温部件所在区域。
图7 不同观测角度下的本征红外辐射亮度分布图
Fig.7 Intrinsic infrared radiance distribution under
different observation angles
机体的本征红外辐射亮度在各观测角度下的最大值差别并不显著，均在110～140 W/(sr·m2)范围内。观测角度不同时红外辐射亮度区域范围和分布差异明显。例如，正侧视观测角度下，红外辐射亮度的特征点为机身前体侧部、进气道边缘和翼缘尖端；前视观测则呈现出两个条状高亮度特征，这是机体头部前缘和进气道前段的高温区域引起的；俯视观测和仰视观测在前端辐射亮度呈现出明显的分布差异，在机身后端则差异较弱。这意味着，不同观测角度下红外辐射特性差异明显，可用于判别目标的特征属性。
4.4 红外光谱辐射强度
图8为4种典型观测角度下类X-51A飞行器的辐射光谱强度曲线，波长为1～14 μm。图8中可以看出，4个观测角度下光谱辐射峰值强度的大小依次为仰视>俯视>侧视>前视，这是由于高温区域的温度水平、温度分布和发射率大小共同决定的。飞行器的下表面具有较大的高温区域，高温意味着较强的红外发射能力。此外，虽然飞行器背部高温区域较腹部低，但在俯视观测中可以观察到飞行器较大的区域，因此其辐射强度峰值较仰视观测略偏低。
图8 不同观测角度下的光谱辐射强度曲线
Fig.8 Spectral radiation intensity under different
observation angles
此外，前视观测下的辐射强度峰值较低，约为2 000 W/(sr·μm)，较其他3种观测角度偏低3倍以上，且光谱峰值约在2.2 μm处，较其他工况的峰值波段(约在3 μm处)而言向短波方向偏移，这是本体壁面发射红外光谱的峰值波长随温度的升高而变短的缘故。
5 结论
1) 高温区域主要集中在飞行器头部、机身前体侧部和翼缘部分，最高温度高达1 430 K。
2) 机体的本征红外辐射亮度在各观测角度下的最大值差别并不显著，观测角度不同红外辐射亮度区域范围和分布区差异明显，可用于判别目标的特征属性。
3) 典型飞行工况下，类X-51A的红外光谱辐射强度峰值随观测角度不同而不同，谱辐射峰值强度的大小依次为仰视>俯视>侧视>前视观测角度。
[1] 邓艳丹，黄生洪，杨基明，等.一种X-51A相似飞行器模型的气动特性初探[J].空气动力学学报，2013，31(3).
Deng Y D，Huang S H，Yang J M，et al.A preliminary study on aerodynamic characteristics of an X-51A similar aircraft model[J].Chinese Journal of Aerodynamics，2013，31(3).
[2] Caitlin H.USAF successfully tests X-51A waverider[J].Jane’s DefenceWeekly，2010，47(22)：130-139.
[3] Rao A G，Mahulikar S P.Effect of atmospheric transmission and radiance on aircraft infrared signatures[J].Journal of Aircraft，2005，42(4)：1046-1054.
[4] Hank J M，Murphy J S，Mutzman R C.The X-51A scramjet engine flight demonstration program[J].AIAA Paper，2008，2540：2008.
[5] 程迪，黄生洪.X-51A 高超声速乘波体气动外形的反向建模和数值模拟[C]//第十四届全国激波与激波管学术会议论文集(上册).中国力学学会.黄山：2010.
Cheng D，Huang S H.Reverse Modeling and Numerical Simulation of the Aerodynamic Shape of X-51A Hypersonic Wave Wider[C]//The 14th National Conference on Shock Waves and Shock Tubers(Volume 1).2010.
[6] 周正，贺旭照，吴颖川.X-51A飞行器模型的建立及性能初步研究[C]//第十五届全国激波与激波管学术会议论文集(上册).中国力学学会.杭州：2012.
Zhou Z，He X Z，Wu Y C.The Establishment of X-51A Aircraft Model and Preliminary Study on Its Performance[C]//The 15th National Conference on Shock Waves and Shock Tubers.2012.
[7] Jiang CHEN，Songping WU.Numerical Simulation Study of Aerodynamic Characteristic on Hypersonic Aircraft[C]//2012 National Conference on Information Technology and Computer Science.Atlantis Press，2012.
[8] Niu Q，Yuan Z，Chen B，et al.Infrared radiation characteristics of a hypersonic vehicle under time-varying angles of attack[J].Chinese Journal of Aeronautics，2019，32(4).
[9] 谈和平，夏新林，刘林华.红外辐射特性与传输的数值计算的数值计算-计算热辐射学[D].哈尔滨：哈尔滨工业大学出版社，2006.
Tan H P，Xia X L，Liu L H.Numerical calculation of numerical calculation of infrared radiation characteristics and transmission computational thermal radiation.Harbin Institute of Technology Press，2006.
[10] 穆磊，马宇，贺志宏，等.高速钝锥体非平衡绕流场光辐射特性模拟[J].工程热物理学报，2012，33(11)：1958-19652.
Mu L,Ma Y,He Z H.Radiation characteristics simulation of non-equilibrium flow field around the hypersonic blunted cone[J].Journal of Engineering Thermophysics,2012,33(11):1958-1962.
[11] Mu L，Ma Y，He Z H，et al.Simulation of optical radiation characteristics of high-speed non-equilibrium flow around a blunt cone[J].Acta thermophysicsSinica，2012，33(11)：1958-1965.
[12] Mutzman R，Murphy S.X-51development：a chief engineer's perspective[C]//17th AIAA international space planes and hypersonic systems and technologies conference.2011，13.
[13] Thivet F，Besbes O，Knight D D.Effect of grid resolution accuracy of skin friction and heat transfer in turbulent boundary layers[C]//AIAA，Aerospace Sciences Meeting and Exhibit，38th，Reno，NV.2000.
[14] Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA journal，1994，32(8)：1598-1605.
[15] Park C.Review of chemical-kinetic problems of future NASA missions.I-Earth entries[J].Journal of Thermophysics and Heat Transfer，1993，7(3)：385-398.
[16] 聂涛，刘伟强.高超声速飞行器前缘流固耦合计算方法研究[J].物理学报，2012，61(18)：184401-1-5.
Nie T，Liu W Q.Computational method of fluid-structure coupling at leading edge of hypersonic vehicle[J].Acta Physica Sinica，2012，61(18)：184401-1-5.