矩形喷管外尾焰红外辐射特性的数值计算

冯云松，李晓霞，路远，金伟

(解放军电子工程学院 安徽省红外与低温等离子体重点实验室，安徽合肥230037)

0 引言

飞机发动机的排气系统是飞机的主要红外辐射源，采用矩形喷管是一种有效抑制排气系统红外辐射的手段。由于矩形喷管在实现隐身性上的优势，轰炸机、战斗机以及无人机都广泛采用了矩形喷管[1]。因此，对矩形喷管外尾焰红外辐射特性研究颇有必要。

目前，国外在飞机排气系统红外辐射特性研究方面的理论、技术和方法都进行了研究，如Varney[2]，Jim[3]对飞机排气系统的红外辐射特性进行了较详细的数学建模和数值计算。国内的学者从20 世纪90 年代开始已采用多种方法对火箭、导弹和飞机等发动机尾焰的红外辐射特性进行了计算。如董士奎等[4]用贴体坐标系下的离散坐标法(DOM)研究了不同工况下尾喷焰的红外光谱辐射特性;樊士伟等[5]用有限体积法(FVM)计算了固体火箭发动机羽流的红外特性;帅永等[6]用反向蒙特卡罗法(BMCM)模拟了高温含粒子自由流的红外辐射特性;阮立明等[7]提出了用源项六流法(SSFM)求解导弹尾喷焰的红外辐射特性。综合比较针对尾焰红外辐射的多种算法，其中有限体积法具有原理简单、计算精度高和可以求解任意方向的辐射特性等优点，因此，本文针对矩形喷管外尾焰，采用有限体积法数值计算其红外辐射特性。

1 喷管外尾焰流场的计算

1.1 喷管几何模型

矩形喷管的几何模型如图1 所示，喷管入口取在XOY 面上，Z 轴正向为射流方向。矩形喷管的宽w=0.9 m，高h =0.45 m，长l =0.5 m[1]，喷管出口的当量直径D*=(4wh/π)1/2=0.718 m.

图1 喷管模型Fig.1 Nozzle model

1.2 数学模型

为了简化尾焰的流动模拟和分析，假设:

1)喷管燃气流动为定常流动;

2)忽略燃气在膨胀过程中组分的变化，并认为燃气的定压比热是常量;

3)假设喷管中燃气膨胀流动过程是理想等熵流动过程，忽略燃气对喷管壁的传热;

4)喷管内流动为纯气相流动，燃气为理想气体，服从理想气体状态方程;

5)不考虑辐射作用，忽略重力的影响。

根据以上假设，包括连续方程、动量方程、能量方程的通用形式[8]为

式中:ρ 为气体密度;v 为气体速度矢量;φ 为通用求解变量;Γφ为有效扩散系数;Sφ为广义源项。

湍流采用RNGκ－ε 模型，湍流动能κ 方程和湍流耗散率ε 方程[8]分别为

式中:ueff为流体在i 方向上的有效速度;Gκ为由层流速度梯度而产生的湍流动能;Gb为由浮力而产生的湍流动能;YM为在可压缩湍流中，过渡的扩散产生的波动;C1z、C2z、C3z为常量;ακ和αε分别为湍流动能κ 和耗散率ε 对应的湍流Prandtl 数。

1.3 尾焰流场数值计算

流场计算区域在Z 方向的长度为14D*，在X和Y 方向的长度均为5D*.图2 给出了流场计算区域的整体示意图，喷管进口位于ABCD 平面上，出口位于EFGHIJKL 组合面上，区域MNOPQRST 为整个计算区域边界。计算中采用六面体结构化网格，流场计算区域网格数目均在13 万左右。

假定喷管入口气体为完全燃烧的燃气，基于燃烧理论，确定主流进口处N2、CO2和H2O 蒸气质量百分比分别为0.706，0.209 和0.085;外场边界和进口中引射的气体均为空气，成分为N2和O2，其质量百分比为0.756 和0.244[9].

设置飞机的飞行高度为3 km，飞行马赫数为1.5，发动机工作状态为非加力状态[1]。入口边界条件:喷管入口为压力入口，给定的总温度为800 K，总压力为0.16 MPa，入口流动角θ = 0°[9－10];计算IJKLMNOP 区域的边界条件:设定该长方体的几个面为压力出口，温度为280 K，总压力为71 kPa;壁面边界条件:采用无滑移固壁边界条件，壁面设定为流、固耦合面，温度场计算时不考虑壁面之间的辐射传热。

图2 流场计算区域Fig.2 Computational area of flow field

计算中，设外界气流静止，尾喷管壁面绝热，忽略尾喷管厚度。流场计算使用Fluent 6.3 软件进行，采用Couple 隐式算法，壁面附近采用标准壁面函数进行修正，收敛精度为10－4.当Fluent 对尾焰流场的迭代计算收敛时，获得尾焰的静压力p 和静温度T 云图，如图3 所示。

图3 Y=0 m 处XOZ 切面的尾焰静压力和静温度分布Fig.3 Distribution of static pressure and static temperature of exhaust plume in XOZ face in Y=0 m

2 气体的辐射传递方程

结合介质的发射、吸收和散射的相互关系，文献[11 －12]建立了求解光谱辐射亮度的辐射传递方程，其表达式如下:

式中:λ 为窄带l 的中心波长;r 为位置向量;s 为方向向量为气体平均吸收系数;σsλ(s)为气体散射系数;Lλ(s，r)为位置r 处微元体在s 方向上的光谱辐射亮度;Lbλ(s)为微元体在s 方向上的黑体光谱辐射亮度;s'为散射方向向量;Φλ(s'，s)为散射相函数;Ωi为空间立体角。

航空发动机在非加力状态下工作，且燃烧较充分时，其尾焰几乎不含固体颗粒，因此，可不考虑散射对辐射的影响。由此，辐射传递方程简化为

根据(5)式，若要求解辐射亮度，必须首先获得传输介质的平均吸收系数为了求解平均吸收系数αλ(s)，先采用洛伦兹线型的统计窄带模型(6)式计算某个窄带l 内的平均透过率[13]

式中:ΔC、k 和d 分别为碰撞增宽的半宽、平均吸收因子和谱线间距;u 为光学厚度。

由于尾焰中含有多种吸收气体，可以先求出每种吸收气体的吸收系数，然后对所有吸收系数求和，即可得尾焰的吸收系数。

3 尾焰红外辐射的计算

3.1 确定核心计算区域

根据对矩形喷管外尾焰流场核心区域的分析，当尾焰的密度和压强在喷口附近取得最大值后，随着喷口径向和轴向距离的增加，尾焰对红外辐射吸收和发射作用的影响越小。可以用指定的CO2百分比含量ηCO2作为设定尾焰辐射计算区域的阈值，并对区域进行简单的填补式修正，得到长方体形尾焰[10]。具体的方法为:在X 轴方向上找到尾焰流场的CO2百分比含量大于等于ηCO2时的最大坐标Xmax，同理在Y 轴方向上和Z 轴方向上找到最大坐标Ymax和Zmax.把Zmax，2Ymax，2Xmax分别作为辐射计算区域的长宽高，从而得到一个长方体形的核心辐射计算区域，如图4 所示。对计算区域进行网格化，NX=8，NY=12，NZ=20，如图5 所示，则每个网格(微控制体)的编码为(nX，nY，nZ).

图4 尾焰核心区域示意图Fig.4 Sketch of plume core area

图5 计算区域网格化示意图Fig.5 Gridding of computational area

3.2 光谱亮度计算

FVM 的基本思想是保证微控制体在每个立体角内的辐射能量守恒[11]。为此，需要对计算区域和4π 空间分别进行空间离散和角度离散。前文已经对尾焰辐射区域进行了空间离散，得到的单个微控制体P，其体积为VP.

对4π 空间进行角度离散，离散为互不重叠的立体角Ωmn，如图6 所示。Ωmn中心对应的天顶角和方位角分别为θm和φn.

对(5)式在控制体P 和立体角Ωmn积分，得到[15]

图6 立体角离散示意图Fig.6 Discrete solid angle

式中:re为辐射亮度沿Ωmn中心方向的单位矢量，在直角坐标系中的计算表达式为

经过分析与计算，尾焰计算区域中的每个控制体都会得到类似于(8)式的方程，从而形成以控制体中心节点沿离散立体角中心方向上的光谱辐射亮度为未知数的方程组，联合迭代求解该方程组，便可求得任一控制体在任何离散立体角中心方向上的光谱辐射亮度。

3.3 红外辐射强度计算

在得到尾焰某一微控制体的方向光谱辐射亮度Lλ(θm，φn)后，还必须利用(10)式和(11)式计算尾焰的某一微控制体外表面在某一方向的光谱辐射强度ΔIλ(θm，φn)和波段辐射强度ΔIλ1－λ2(θm，φn).通过对尾焰所有控制体外表面同一方向的辐射强度求和，如(12)式，进一步得到尾焰在该方向的波段总辐射强度Iλ1－λ2(θm，φn).

式中:φ 为探测方向与微控制体外表面法线的夹角;ΔA 为微控制体外表面的面积。

式中:λ1－λ2为某一波段;A 为尾焰核心计算区域所有外表面。

4 计算结果及分析

4.1 红外辐射光谱特性

随机取核心计算区域表面的一编码为(8，5，12)的微控制体，其外侧表面中心沿θ=5π/28，φ=π/24 rad 和θ=13π/28 rad，φ =11π/24 rad 方向上的光谱辐射亮度曲线如图7 所示。在图7 的两条光谱辐射亮度曲线上，都存在几处较为明显的辐射峰，辐射峰值波长分别在2.7、4.3、5.97 μm 和6.55 μm附近，CO2在2.7 μm、4.3μm 处附近产生了2 个强吸收带，H2O在2.7、5.97 μm 和6.55 μm 处附近产生了3 个强吸收带，其中，左起第一处辐射峰2.7 μm的出现是H2O 与CO2共同作用的结果。在图7(a)中，辐射亮度的最大值为6 ×10－2W/(cm2·μm·sr);而在图7(b)中，辐射亮度的最大值为9 ×10－2W/(cm2·μm·sr)，后者约是前者1.5 倍，产生这种现象的原因是:随着辐射方向的变化，辐射传输经过尾焰区域的温度和光学程长是变化的，导致光谱辐射亮度发生变化。

图7 控制体(8，5，12)在不同方向的光谱辐射亮度Fig.7 Spectral radiant luminance of control volume(8，5，12)in different directions

4.2 红外辐射强度分布

对整个核心计算区域在宽边对称面XOZ 和窄边对称面YOZ 面内的探测范围均为－13π/28 ～13π/28 rad、探测间隔为π/14 rad 的各个方向利用(12)式，计算得到尾焰在3 ～5 μm 波段红外辐射强度分布，如图8 所示。图8(a)和图8(b)为尾焰分别在XOZ 与YOZ 平面内的红外辐射强度分布。在XOZ 平面内，尾焰在θ =π/28 rad 方向的红外辐射强度IXOZ仅为267 W/sr，随着θ 的增大，尾焰红外辐射强度迅速增大，当在θ =13π/28 rad 方向时，红外辐射强度IXOZ增大到1 965 W/sr，红外辐射强度增加了7 倍;在YOZ 平面内，尾焰在θ=π/28 rad 方向的红外辐射强度IYOZ仅为142 W/sr，随着θ 的增大，尾焰红外辐射强度迅速增大，当在θ =13π/28 rad 方向时，红外辐射强度IYOZ增大到69 W/sr，红外辐射增加了约5 倍;YOZ 平面内最大红外辐射强度IYOZ，max仅是XOZ 平面内的最大红外辐射强度IXOZ，max的35%，这主要因为矩形喷管的宽高比AR 为2∶1，导致尾焰为扁平状，所以造成了宽边对称面内的IXOZ，max大于窄边对称面内的IYOZ，max.

图8 尾焰3 ～5 μm 波段总辐射强度分布Fig.8 Total intensity distribution in 3 ～5 μm of exhaust plume

5 结论

本文在只考虑燃气中H2O、CO2的光谱吸收与发射影响的情况下，针对矩形喷管外尾焰流场进行了数值模拟计算，并采用FVM 和气体辐射的窄谱带模型，计算了矩形喷管尾焰的红外辐射的光谱特性与在3 ～5 μm 波段的总强度分布。最终的模拟计算结果与文献[9 －10]所得实验数据相吻合，证实了模型建立的合理性与数值计算结果的有效性。并由结果分析，初步得出以下结论:

1)在大气窗口3 ～5 μm 内尾焰辐射出现了2.7 μm和4.3 μm 两个辐射峰，且4.3 μm 处的光谱辐射亮度较2.7 μm 处的要大，因此一般采用中心工作波长为4.3 μm 的红外探测器探测飞机的尾焰辐射。

2)矩形喷管的尾焰为扁平状，且宽边对称面内的红外辐射强度要大于窄边对称面内的红外辐射强度，因此红外探测器在尾焰宽边对称面内将更容易探测到飞机的中波辐射。

3)本文的计算方法与结果可为实现飞机红外隐身和红外探测隐身飞机提供理论依据。

在实际的飞机尾焰中，必然有固体粒子存在(如未完全燃烧的碳粒子)，也就存在固体粒子对红外辐射的散射，但是本文在求解辐射传递方程时，忽略了散射项，这必将造成一定的计算误差，因此，本文数值计算结果精确度有进一步提高的空间。

References)

[1]魏钢.F－22“猛禽”战斗机[M].北京:航空工业出版社，2008:5 －8，80 －81.WEI Gang.F－22“Accipiter”fighter plane[M].Beijing:Aviation Industry Press，2008:5 －8，80 －81.(in Chinese)

[2]Varney G E.Infrared signature measurement techniques and simulation methods for aircraft survivability，AIAA－1979－1186[R].Las Vegas:AIAA，SAE and ASME Joint Propulsion Conference，1979.

[3]Jim C.F/A－22 IR signature flight test model validation[J].Aircraft Survivablity，2003，29(4):9212 －9216.

[4]董士奎，于建国，李东辉，等.贴体坐标系下离散坐标法计算尾喷焰辐射特[J].上海理工大学学报，2003，25(2):159－162.DONG Shi－kui，YU Jian－guo，LI Dong－hui，et al.Numerical modeling of infrared radiation properties of exhaust plume by the discrete ordinates method in body－fitted coordinates[J].Journal of University of Shanghai for Science and Technology，2003，25(2):159 －162.(in Chinese)

[5]樊士伟，张小英，朱定强，等.用FVM 法计算固体火箭羽流的红外特性[J].宇航学报，2005，26(6):793 －797.FAN Shi－wei，ZHANG Xiao－ying，ZHU Ding－qiang，et al.Calculation of the infrared characteristics of the solid rocket plume with FVM method[J].Journal of Astronautics，2005，26(6):793 －797.(in Chinese)

[6]帅永，董士奎，刘林华.高温含粒子自由流红外辐射特性的反向蒙特卡罗法模拟[J].红外与毫米波学报，2005，24(2):100－104.SHUAI Yong，DONG Shi－kui，LIU Lin－hua.Simulation of infrared radiation characteristics of high temperature free－stream flow including particles by using backward Monte－Carlo method[J].Journal Infrared Millimeter and Waves，2005，24(2):100 －104.(in Chinese)

[7]阮立明，齐宏，王圣刚，等.导弹尾喷焰目标红外特性的数值仿真[J].红外与激光工程，2008，37(6):959 －962.RUAN Li－ming，QI Hong，WANG Sheng－gang，et al.Numerical simulation of the infrared characteristic of missile exhaust plume[J].Infrared and Laser Engineering，2008，37(6):959 －962.(in Chinese)

[8]李进良，李承曦.精通Fluent 6.3 流场分析[M].北京:化学工业出版社，2009:55 －58.LI Jin－liang，LI Cheng－xi.Mastery of flow field analysis in Fluent 6.3[M].Beijing:Chemical Industry Press，2009:55 －58.(in Chinese)

[9]范仁钰，张靖周，单勇.遮挡板结构对二元喷管红外辐射特性的影响[J].航空动力学报，2011，26(2):333 －348.FAN Ren－yu，ZHANG Jing－zhou，SHAN Yong.Effects of sheltering baffles on the infrared radiation characteristics of two－dimensional nozzles[J].Journal of Aerospace Power，2011，26(2):333－348.(in Chinese)

[10]侯晓春，季鹤鸣，刘庆国，等.高性能航空燃气轮机燃烧技术[M].北京:国防工业出版社，2002:375 －378.HOU Xiao－chun，JI He－ming，LIU Qing－guo，et al.Combustion technology for high performance aviation gas turbine[M].Beijing:National Defense Industry Press，2002:375 － 378.(in Chinese)

[11]Siegel R，Howell J R.Thermal radiation heat transfer[M].Washington DC:Hemisphere and McGraw－Hill，1981.

[12]王伟臣，李世鹏，张峤，等.工作条件对固体发动机羽流温度场的影响[J].兵工学报，2011，32(12):1493 －1497.WANG Wei－chen，LI Shi－peng，ZHANG Qiao，et al.Influence of operating conditions on temperature distributions of exhaust plume of solid rocket motor[J].Acta Armamentarii，2011，32(12):1493 －1497.(in Chinese)

[13]Ludwing C B，Malkus W，Reardon J E，et al.Hand－book of infrared radiation from combustion[R].US:NASA－SP－3080，1973:222 －230.

[14]Liu F，Gülder Ö L，Smallwood G J.Three－dimensional non－grey gas radiative transfer analyses using the statistical narrow－band model[J].International Journal of Heat Mass Transfer，1998，37(9):759 －768.

[15]谈和平，夏新林，刘林华，等.红外辐射特性与传输的数值计算[M].哈尔滨:哈尔滨工业大学出版社，2006:194 －198.TAN He－ping，XIA Xin－lin，LIU Lin－hua，et al.Numerical calculation of transmission and characteristic of infrared radiation[M].Harbin:Harbin Institute of Technology Press，2006:194－198.(in Chinese)