固发羽流流场及辐射特性随飞行状态的变化规律

王志博 宋兆龙 张 彪 李 健 许传龙

(1东南大学能源与环境学院， 南京 210096)(2东南大学火电机组振动国家工程研究中心， 南京 210096)

固体火箭发动机在地面搭载实验或者高空低压飞行时，发动机尾喷管喷射出的固体燃料燃烧产物的火焰流呈现羽毛状，其称为羽流.固发羽流主要来源于推进剂的燃烧，少部分来源于喷管材料的烧蚀、衬层和绝热层的裂解，这种高温高压的羽流会向空间中辐射出强烈的红外特征信号.掌握固发羽流的流场和辐射特性，可为制导反导系统探测、识别及反识别提供方案，同时还可为飞行器跟踪、拦截提供依据.

近年来，许多学者对羽流进行了大量的实验和模拟仿真研究.在羽流复燃流场计算方面，牛青林等[1]针对飞行状态对固体火箭发动机尾焰的影响，根据热力计算结果对某固体火箭发动机的化学反应流和单一组分流进行了仿真计算.吴睿等[2]采用常用的3种k-ε湍流模型和2阶/3阶迎风格式，对单UDMH/NTO火箭尾焰随飞行高度和来流速度的影响进行了数值仿真，得出尾焰在不同飞行高度、不同来流马赫数下流场变化规律以及复燃反应对组分分布的影响规律.郭东升等[3]研究了不同飞行马赫数下Al2O3凝相粒子的固体火箭发动机的流场结果，并与纯气相流场结构进行了比较，Al2O3凝相粒子的加入会使得波系膨胀加强，长度增加.目前关于羽流流场及辐射的计算，多是在一定飞行高度及飞行速度上的研究，对于实际火箭升空过程中流场及辐射特性的研究较少.

在羽流辐射特性计算方面，Rankin等[4]采用高速中红外相机对湍流非预混射流火焰的红外辐射强度进行了实验测量，并验证了大涡模拟(LES)与窄带辐射模型.李霞等[5]建立了喷焰复燃数值模拟方法和视线光辐射传输数值计算方法，结合单线组谱带模型(SLG)对不同高度下的喷焰进行了模拟及辐射特性计算.Gu等[6]采用灰气体加权和模型(WSGGM)并结合有限体积法研究了羽流流场及辐射特性，发现了辐射热流随飞行高度的变化规律.Cai等[7]采用有限体积法并利用Elsasser窄带模型计算了火箭羽流中燃烧气体的辐射特性，并利用Mie散射理论计算了颗粒的辐射特性.王扶辉等[8]基于HITEMP数据并结合全光谱K分布(FSK)模型对平板内的混合气体介质辐射情况进行了计算，并与其他方法进行比较，验证了FSK模型的准确性.目前FSK模型已经成功应用在多种工况，用于计算气体的光谱辐射特性，相对于其他模型而言其计算速度更快.

本文研究了固体火箭在不同飞行高度和来流马赫数下羽流流场及辐射特性的变化规律，并进一步计算了实际火箭升空过程中随着时间推移羽流流场及辐射特性的变化.利用FLUENT软件对某固体火箭发动机羽流进行计算，得到羽流流场分布，继而利用Mie理论和HITEMP数据库结合FSK模型计算了粒子和气体的辐射参数，使用广义源项有限体积法计算羽流红外辐射特性.

1 计算模型及方法

1.1 羽流复燃流场计算方法

羽流中含有大量未完全燃烧的成分，在进入大气后，由于卷吸效应，会产生复燃反应，因此本文选用含有组分运输项和化学反应项的二维轴对称N-S方程，表示如下：

(1)

式中，t为时间；x和y分别为喷焰的轴向和径向坐标；Q为守恒变矢量；Ec、Fc分别表示坐标系2个方向上的对流通量；Ev、Fv分别为2个方向上的黏性通量；S为源项.

羽流中的复燃反应主要发生在近场与空气的混合处，该处气体为超声速区域，因此选用有限速率化学反应模型.根据Arrhenius公式来计算化学反应源项，其中第r个反应的化学方程式为

(2)

式中，N为系统中化学物质数目；v′τ,r、v″τ,r分别为第r个反应中组分τ的反应物和生成物的当量系数；Mτ为组分τ的相对分子质量；kf,r、kb,r分别为反应r的正、逆反应速率常数.

反应r的正向反应速率通过Arrhenius公式计算：

(3)

式中，Ar为指数前因子；T为温度，K；βr为温度指数；Er为反应活化能，kJ/mol；R为通用气体常量，J/(mol·K).

为求解以上控制方程，采用雷诺平均方法建立计算模型，并且利用基于涡黏性假设的带旋流修正的k-ε湍流模型中的雷诺应力和雷诺输运项.为了模拟羽流中的复燃反应，参考文献[9]中的13组元30个反应化学动力模型，如表1所示.

表1 羽流中主要复燃反应模型

1.2 羽流红外辐射特性计算方法

固发羽流红外辐射特性主要与其流场参数(温度、压力)和成分有关，其成分包括凝聚相微粒(Al2O3等)和气体分子(H2O、CO2、CO、HCl等).要计算羽流红外辐射特性，需要确定流场介质的吸收系数、散射系数等物性参数.

关于气体辐射物性参数的计算，本文采用FSK模型计算.FSK模型是指整个光谱上普朗克函数不变，将吸收系数进行重新排列，并根据吸收系数出现的几率进行统计，获得吸收系数在该光谱范围内出现的几率分布函数f(κ)；若进一步统计小于κ的吸收系数出现的几率，则可获得吸收系数的累积分布函数g(κ).由于累积分布函数g(κ)是单调上升的光滑函数，可采用简单高效的高斯数值积分求解[10].

本文以HITEMP 2010光谱数据库为基础，将1.5～5.5 μm谱带均匀分为40个波数区间.同时，考虑的压力范围为0～0.4 MPa，间隔0.05 MPa进行取值；温度范围为300～3 000 K，间隔300 K进行取值.计算每种气体在不同光谱波数区间、温度和压力下的吸收系数，得到12个Guass积分点对应的吸收系数矩阵.根据吸收系数矩阵，任意状态下每种气体在某个高斯点的吸收系数函数形式为

(4)

ygi(T)=a1T6+a2T5+a3T4+a4T3+a5T2+a6T+a7

(5)

(6)

式中，z为吸收气体摩尔分数，取值范围为0～1；gi为第i个高斯点，i=1,2, …, 12；κgi,gas为某种气体在第i个高斯点处的吸收系数；P为实际状态下的大气压，Pa；Pe为标准大气压，Pa；ap、bq为多项式拟合系数，p=1,2,…,7，q=1,2,3.由于CO2、CO、HCl谱线的自加宽半宽和空气加宽半宽相当，摩尔分数的变化对吸收系数温度拟合影响很小，因此3种气体在任意摩尔分数条件下的吸收系数是通过理想状态即摩尔分数为1.0时的吸收系数乘以实际状态下摩尔分数数值所得.而对于H2O谱线而言，其自加宽半宽是空气加宽的5倍，因此对摩尔分数0.01、0.1、1的水蒸气进行拟合，再乘以其摩尔分数.混合气体的K分布吸收系数κgi可由单一气体的K分布吸收系数直接相加得到[11].

对于大多数固体火箭发动机而言，羽流中含有部分固体颗粒且粒子为稀疏粒子群，粒子之间独立散射，故将其看作球形粒子，即可利用Mie散射理论计算其消光、散射和吸收系数，如下式所示：

(7)

式中，下标x′=e,s,a;κη,e、κη,s、κη,a分别为波长η时颗粒的消光、散射和吸收系数，m-1；Qx′,n为对应因子；Dn、Nn分别表示第n种粒径粒子的直径和数密度.

在均匀温度、浓度环境中，带有普朗克权重的全光谱K分布函数可以表示为

(8)

式中，Ib,η表示光谱辐射强度，W/(m2·μm·sr)；Ib表示黑体辐射强度，W/(m2·μm·sr)；δ(κ-κη)为狄拉克函数；κη表示波长为η下的吸收系数.

累积K分布函数g(κ) 表示为

(9)

考虑谱带模型的吸收、发射、散射性非灰介质的辐射传递方程[12]，两边同乘δ(κ-κη)，在全光谱区间上积分后再分别除以几率分布函数f(κ)，方程转变为如下形式：

(10)

式中，Igi为累计分布函数对应的辐射强度，W/(m2·μm·sr)；l为微元段的长度；Ωm表示散射出射方向，Ωm′表示辐射入射方向；Iη(l,Ωm)为Ωm方向、l微元段的光谱辐射强度，W/(m2·μm·sr)；Ib，Δη(l,Ωm)为l微元段的黑体光谱辐射强度，W/(m2·μm·sr)；Δη为总波数间隔；Φ(Ωm,Ωm′)为散射项函数.

(11)

总辐射强度采用12点Gauss积分来计算：

(12)

式中，ωi为Gauss积分点对应的权值.

2 计算实例及模型验证

2.1 计算模型

考虑到固发羽流流场具有轴对称的特性，故羽流流场模型选择二维轴对称模型对流场进行数值计算.本文研究对象为BEM-2发动机，喷管的几何尺寸及喷口参数参考文献[13].其中，喷口尺寸为15 mm，计算域大小为X轴方向0.6 m，Y轴方向2.5 m，网格数目为58 254.图1和图2分别为流场模型简图和模型网格分布示意图，边界条件如表2所示.图2中，由于喷管壁面以及流场轴线处流场

图1 流场模型简图

特性变化相对比较剧烈，因此选择网格局部加密.

图2 模型网格分布示意图

表2 边界条件参数

辐射计算域是将羽流流场旋转为三维结构，即得到圆柱体的羽流三维模型，建立长方体，沿x、z、y方向节点数分别为30、30、200，并分割为立方体.以立方体中心到轴线距离为基准，匹配圆柱体中各点流场参数，该点参数代表立方体中流场参数.继而，采用广义源项有限体积法求解辐射传输方程，并获得射线路径上的辐射强度，如图3所示，通过谱带积分即可在探测器像素上获得总的辐射强度.

图3 羽流辐射传输示意图

2.2 模型验证

为验证本文计算模型的正确性，利用本文建立的辐射传输模型，同时参考文献[13]中流场边界条件以及化学反应模型，计算了发动机尾焰流场及光谱红外辐射模型，并与文献[13]中实验及仿真数据进行了对比，如图4所示.

图4 本文光谱辐射强度计算值与文献[13]测量值、计算值对比曲线

固体火箭发动机羽流中存在多种气体组分，气体对不同波长具有选择性吸收的特性，并且辐射强度的数值与气体组分的浓度有关.从图4可以看出，在1.5～5.5 μm的波长范围内，光谱辐射强度随波长的整体变化规律基本相同，仅在部分点存在差异，本文计算结果与文献[13]结果均在2.5～3.0 μm及4.2～4.7 μm的波长范围内出现了较强的光谱辐射强度，而在其他波长范围内较弱.由于本文与文献[13]中采用的计算模型有一定差异，且羽流复燃模型涉及范围比较大，以上误差在可接受的范围之内，因此本文模型是正确可靠的.

对于很多固体火箭发动机而言，在羽流中存在大量Al2O3颗粒.因此，在计算光谱辐射强度时，需加入Al2O3固体粒子，本文假设Al2O3固体粒子的直径为10 μm.图5表示含有Al2O3颗粒与未含有Al2O3颗粒(仅考虑气体辐射)羽流光谱辐射强度对比曲线.从图中可以看出，在4.1～4.5 μm波段下考虑Al2O3粒子的影响时误差相对较小，与文献[14]在4.5 μm波段下加入Al2O3颗粒造成羽流辐射强度的误差数值仅为1%的结论相符，侧面证明本文粒子散射辐射模型的正确性.在其他多数波段下，考虑Al2O3颗粒对羽流光谱辐射强度影响较大，如在1.6～1.7 μm波长范围内，未考虑Al2O3辐射时羽流光谱辐射强度相对于考虑Al2O3时下降了99.34%.因此，在研究羽流光谱辐射强度时，应考虑Al2O3颗粒的吸收与散射.

图5 含有Al2O3颗粒与未含有Al2O3颗粒羽流光谱辐射强度对比

3 结果与分析

3.1 飞行参数介绍

考虑到本文不同飞行状态对于流场的影响较大，为了保证羽流充分发展，在此计算域的尺寸为1 m×8 m，网格数目为47 588.基于本文上述流场仿真模型以及化学反应动力学模型，分别对不同状态下的羽流流场进行模拟仿真.不同飞行高度H下的环境参数参考美国标准大气USSA—1976数据，部分高度大气条件如表3所示.考虑到固体火箭在飞行过程中，随着时间的推移，飞行高度与来流马赫数随之变化，参考文献[15]中火箭飞行轨迹，选取典型轨迹点的数据，其高度大气条件如表4所示.

田志芳一个人在地窝子里，看看头顶，看看脚下，一屁股坐在土台上，叹口气，心想姆妈拼死阻拦都没拦住她和哥哥，现在怪谁呢，自己跳起脚要支边。她垂下头，把手中的沙枣花捧起来瞧，带沙点的叶根处，确实有细小的花苞，同样泛出密密麻麻的沙尘，形如青色的米粒，一粒一粒挤在一起，好似家乡中秋的桂花。猜想，沙枣花开了，是不是真有桂花那样的千里香？

表3 不同高度下环境参数

表4 不同轨迹点下环境参数

3.2 飞行高度对羽流流场的影响

固体火箭在飞行过程中，飞行高度变化，环境压力和温度也会发生变化，对羽流的复燃产生很大影响，进而对流场光谱辐射强度产生影响.基于以上辐射计算模型，分别对飞行高度为0、5、15、20 km，来流马赫数Ma=1.0下的羽流流场进行模拟仿真.图6表示来流马赫数为1.0时不同飞行高度下的温度云图.从图中可以看出，随着飞行高度的增加，环境压强逐渐下降，羽流的膨胀更加充分和明显；同时，核心区域的位置逐渐后移，核心区域最高温度相对呈现下降趋势，羽流压缩气体重复出现几次膨胀—压缩—膨胀，从而温度达到峰值.

图7与图8分别表示不同飞行高度下的光谱辐射强度随波长变化曲线以及4.3～4.4 μm波长范围内的光谱辐射强度图像.从图中可以看出，随着飞行高度的增加，羽流的光谱辐射强度逐渐增加，尤其是H=20 km时尤其明显.若以H=20 km时的光谱辐射强度为标准，在4.3～4.4 μm波长范围内，羽流的光谱辐射强度下降幅度最大，H=0，5，15 km时光谱辐射强度下降分别为96.05%、86.45%、41.16%.分析光谱辐射强度曲线以及辐射强度图像可以得出，随着飞行高度的增加，虽然核心区域温度峰值减小，但是羽流核心温度区域逐渐增大，导致羽流的核心区域光谱辐射强度增强.飞行高度为0时，即使核心区域温度峰值较高，但

(a) H=0

(b) H=5 km

(c) H=15 km

(d) H=20 km

图7 Ma=1.0时不同飞行高度下的光谱辐射强度曲线

是由于核心区域面积太小，导致光谱辐射强度下降很多.而4种飞行高度下辐射强度的峰值分别为13 643.0、29 139.4、35 882.0、30 599.7 W/(m2·sr)，因此随着飞行高度的增加，局部的辐射强度呈现先增大后减小的规律.

本文中辐射计算模型为完全对称模型，因此以7.5°为间隔，选取与流场x轴正向不同夹角θ的探测角度进行探测，即流场轴线呈0°、7.5°、15.0°、22.5°、30.0°、37.5°、45.0°.图9表示H=15 km，Ma=1时不同探测角度下的羽流光谱辐射强度曲线.从图中可以看出，随着探测角度的增加，羽流光谱辐射强度呈现下降的趋势，若以θ=0°下的光谱辐射强度为标准，在波长范围为4.2～4.3 μm时，

(a) H=0

(b) H=5 km

(c) H=15 km

(d) H=20 km

探测角度7.5°、15.0°、22.5°、30.0°、37.5°、45.0°下的光谱辐射强度分别下降0.73%、2.95%、6.73%、11.77%、16.54%、18.63%.

图9 H=15 km，Ma=1.0时不同探测角度下的羽流光谱辐射强度曲线

3.3 来流马赫数对羽流流场的影响

固体火箭在飞行过程中，飞行速度会随着飞行状态的改变而改变，但由于其与来流具有相对运动，因此研究飞行速度等效于研究来流马赫数.图10表示来流马赫数分别为0.5、1.0、1.5、2.0下的温度云图.从图中可以发现，随着来流马赫数的增加，环境与羽流之间动压增加，导致流场尺寸逐渐压缩，羽流中的高温区域也逐渐减少.随着来流马赫数的增加，来流中的氧气速度也会加快，与羽流掺混的强度减小，羽流复燃反应强度也随之减弱.同时，由于来流空气的温度远小于羽流温度，自由来流会带走羽流中的热量.因此，随着来流马赫数的增加，热量的交换也会加剧，羽流中温度峰值随之减小，同时羽流周围动压增加，并且峰值的位置随之后移.

(a) Ma=0.5

(b) Ma=1.0

(c) Ma=1.5

(d) Ma=2.0

图10 H=15 km时不同来流马赫数的温度云图

图11表示H=15 km时不同来流马赫数下的光谱辐射强度随波长变化的曲线.由图中可以看出，随着来流马赫数的增加，羽流光谱辐射强度逐渐减小.以Ma=0.5时的光谱辐射强度为标准，羽流的光谱辐射强度在波段4.2～4.3 μm处下降的幅度最大，Ma=1.0，1.5，2.0时光谱辐射强度下降分别为46.56%、73.35%、72.79%.分析光谱辐射强度曲线可以得出，随着来流马赫数的增加，核心区域温度峰值减小，并且羽流核心温度区域逐渐减小，区域减小会减弱羽流的光谱辐射强度.

图11 H=15 km时不同来流马赫数下的光谱辐射强度曲线

图12 H=15 km，Ma=0.5时不同探测角度下的羽流光谱辐射强度曲线

3.4 飞行轨迹对羽流流场的影响

在实际固体火箭发射过程中，飞行高度及来流马赫数会随着时间推移而增大.因此，研究固体火箭随时间变化对羽流流场和辐射特性的影响就显得尤为重要.图13表示羽流不同时间飞行轨迹下

(a) t=13.0 s

(b) t=22.4 s

(c) t=29.3 s

(d) t=34.8 s

(e) t=41.1 s

图13 不同轨迹点下的温度图像

的温度分布云图.从图中可以看出，随着时间推移，飞行高度以及来流马赫数随之增加，从喷口处至第一个马赫盘处，温度逐渐下降达到最小值，并且最小值逐渐减小.在第一个马赫盘结束后，轴线温度逐渐上升，但由于来流马赫数过大，冷空气卷吸更多的热量以及复燃反应强度的降低.再者随着时间推移，羽流中的高温区域逐渐减少，并且羽流流场区域范围增大，这是由于环境压强降低，尾焰流场中的气体向外膨胀.

图14表示不同轨迹点下的羽流光谱辐射强度曲线.从图中可以看出，随着时间的增加，羽流的光谱辐射强度呈逐渐减小的趋势，尤其在t=41.1 s，H=50 km时，在整个波长范围内光谱选择性就不明显，这是由于此时复燃反应强度降低，在流场中温度相对较低.可见，在固体火箭发动机发射过程中，随着时间的增加，飞行高度以及来流马赫数随之增加，会对羽流的流场产生影响，从而使得羽流的光谱辐射强度降低.

图14 不同轨迹点下的羽流光谱辐射强度曲线

图15为t=13.0 s，H=5 km时不同探测角度下的羽流光谱辐射强度曲线.从图中可以看出，随着探测角度的增加，羽流光谱辐射强度呈现递减的趋势.若以θ=0°下的光谱辐射强度为标准，在波长范围为4.3～4.4 μm时，探测角度7.5°、15.0°、22.5°、30.0°、37.5°、45.0°下的光谱辐射强度下降

图15 t=13.0 s, H=5 km时不同探测角度下的羽流光谱辐射强度曲线

均是最大，分别为0.78%、6.41%、14.58%、23.17%、33.77%、41.89%.

4 结论

1) 当考虑Al2O3固体粒子时，相较于仅考虑气体影响的羽流光谱辐射强度，在多数波长范围内影响较大，如在1.6～1.7 μm波长范围内，未考虑Al2O3颗粒辐射时羽流光谱辐射强度相对于考虑Al2O3颗粒辐射时下降了99.34%.

2) 在来流马赫数一定时，随着飞行高度的上升，温度曲线的峰值逐渐减小，峰值位置后移.同时羽流光谱辐射强度逐渐增加，在不同波段下的增强效果有差异.若以Ma=0.5时的光谱辐射强度为标准，在波长为4.2～4.3 μm时，羽流的光谱辐射强度下降的幅度最大，Ma=1.0,1.5,2.0时光谱辐射强度下降分别为46.56%、73.35%、72.79%.

3) 在飞行高度一定时，随着来流马赫数的增加，温度峰值减小，曲线峰值的位置后移.同时羽流光谱辐射强度逐渐减小，若以H=20 km时的光谱辐射强度为标准，在波长为4.3～4.4 μm时，羽流的光谱辐射强度下降幅度最大，H=0,5,15 km时光谱辐射强度下降分别为96.05%、86.45%、41.16%.

4) 在火箭升空过程中，随着时间推移，温度峰值逐渐减小趋于平缓，同时光谱辐射强度急剧减小，到达一定时间后，复燃反应强度极低，光谱辐射强度也趋于平缓.