固体火箭发动机羽流场波系结构及特性数值计算

2018-08-31崔立堃

兵器装备工程学报 2018年12期

崔立堃，叶伟

(陕西理工大学机械工程学院，陕西汉中 723001)

固体火箭发动机羽流流场的结构分布特性对固体火箭的总体设计、火箭发射装置设计、结构材料的选择等方面均会产生很大的影响[1-2]，特别是其红外辐射对发动机故障监控和军事侦察有着非常重要的作用。为此，国内外许多学者对此进行了研究：H.B.Ebrahimi等[3]运用GPACT模型计算分析了Titan II 液体火箭发动机红外波段辐射场；K.Berer等[4]比较了运用逐线模型和分子波带模型计算弹道导弹羽流红外辐射的结果；Kim等使用改进的离散坐标法[5]和有限体积法[6]计算了固体发动机羽流辐射和喷管探照灯(searchlight)发射的底部加热效应；张小英等[7]研究了某型固体火箭发动机在地面试车和不同高度工况下发动机红外波段2～6 μm的辐射情况；金云飞等[8]基于羽流辐射的监测和光谱特性，对发动机的工作状态、羽流特性等进行了在轨观测试验，试验结果良好；陈杰等[9]采用差分求解N-S方程与直接模拟蒙特卡洛耦合的方法对火星环绕器姿轨控发动机羽流力、热以及污染效应进行了仿真研究。这些研究主要集中在羽流红外辐射的计算分析方法和特性探究方面，而对羽流场的流场结构及其特性的研究较少。而羽流场的流场结构及特性对包括羽流红外辐射等其他羽流效应有着重要的意义。因此，本文基于三维N-S方程，采用Spalart-Allmaras湍流模型，通过大型CFD软件FLUENT对某固体火箭发动机羽流场进行了数值分析，并对不同飞行参数条件下、纯气相作用下的发动机羽流场进行了综合分析比较，拟为深入研究固体火箭发动机羽流效应提供一定的指导。

1 数值计算模型

可压缩流动的强守恒型N-S方程在笛卡尔坐标系中可以表示为：

(1)

式(1)中各系数的具体含义及表达式见文献[10]。

在动量方程中，涉及到湍流粘性系数μ1，它必须通过湍流模型求得，本文采用工程上广泛使用的Spalart-Allmaras模型，该模型是一个相对简单的模型，只求解一个有关涡粘性的输运方程，计算量相对较小。此模型是专门应用于涉及到束缚壁面流动的航空领域，对于受到反压梯度作用的边界层，能给出很好的模拟结果。湍流粘性系数μ1按下式计算：

(2)

式(2)中，各系数的具体含义及表达式见文献[11]。

1.1 几何模型及网格划分

本文以某舰载导弹固体火箭发动机为例，选定发动机几何尺寸。整个弹体结构简图及其网格划分如图1所示。发动机喷管主要几何尺寸参数见表1。

表1 喷管尺寸参数

高质量的网格划分能够提高计算精度、加快计算速度，加大计算稳定性。本文中计算域形状不规则，因此为了提高网格质量将其分割为10个部分，

分别对其进行结构化网格划分。由于粘性计算和激波捕捉的需要对喷管壁面附近和尾流场轴线附近进行了网格加密。总的网格数为359 577个。网格划分情况如图2所示。

模拟介质为某推进剂燃烧产物的成分，推进剂按照铝18.2%、过氯酸铵67.2%、丁羟11%、葵二酸二辛酯3.1%、甲苯二异氰酸酯0.5%配比，喷管扩张比是9.5，燃烧组分的计算采用布林克莱法，详细的计算公式和方法及程序参见文献[11]。表2中列出了计算后得到的主要组分及其摩尔数。

表2 燃气中各组分含量

流场计算中取所占比重最大的前七种组份，分别为：CO、HCL、H2、N2、H2O、CO2、H，这七种物质占气态产物总量的98.461%。燃烧产物中Al2O3虽然质量占比比较大，但它是凝相，在仿真发动机羽流场波系结构时可忽略。

1.2 边界条件

本文数值计算中涉及到的边界条件有压力入口边界、壁面边界、对称边界和压力远场边界。其中，喷管入口采用压力入口边界条件，总压5 MPa，总温3 400 K，流动方向与入口边界垂直；喷管壁面采用绝热无滑移固体壁面边界条件；羽流场径向截面采用对称边界条件。外边界采用压力远场边界条件，压力与环境压力一致，采用标准大气相应高度下的参数(见表3)。

很多人心里都不舒服，说来说去的很长一段时间，之后又都闭了口。试问一下，别人除过自己的工作，去掏厕所把自己弄得臭烘烘连老婆孩子都不愿近身了？去牵着驴顶着毒日头口干舌燥满山沟地给它找吃的去了？还是在假期本该休息时一个人寡兮兮地在学校修理那些被淘气孩子损坏的课桌？还是给同事和村里人的家具义务上油漆，给学校充大工修理房子，当农工苦呛呛地去垦荒田了？

表3 大气压强和高度的关系

2 计算结果分析

2.1 相同飞行高度和攻角，不同飞行马赫数羽流场数值分析

选取飞行高度为0 km，攻角为0°，分别对飞行马赫数为0.6、0.9、1.2、1.5这4种情况下的发动机羽流场进行数值计算。其羽流场马赫数等值线图如图3所示。

由图3可以得出，在地面状态时不同飞行马赫数下羽流流场结构分布的变化规律：当飞行马赫数为0.6时，燃气在喷管后方羽流区域形成了五个明显的膨胀-压缩波系，每个波系的长度约为0.8 m，波系的强度逐次降低；当飞行马赫数增加到0.9时，喷管后方羽流区域仍然有五个明显的膨胀-压缩波系，波系的膨胀半径略有减小，波系的长度、强度等方面变化不大；当飞行马赫数增加到1.2时，喷管后方羽流区域形成了四个明显的膨胀-压缩波系，波系的膨胀半径明显减小，波系的长度略有增加，约为0.9 m，羽流核心区的长度有所增加；当飞行马赫数增加到1.5时，流场结构发生了明显的变化，喷管后方羽流区域只有一个明显的膨胀-压缩波系，波系膨胀半径较亚声速时明显降低，第二个和第三个波系未能充分发展，羽流核心区的长度继续增加。

2.2 相同飞行马赫数和攻角，不同飞行高度羽流场数值分析

选取飞行马赫数Ma为0.9，攻角α为0°，分别对飞行高度H为0 km、5 km、15 km这3种情况下的发动机羽流场进行数值计算，结果如图4所示。

由图4可得出，在亚音速状态时，不同飞行高度下羽流流场结构分布的变化规律如下：

当飞行高度为0 km时，燃气在喷管后方羽流区域形成了五个明显的膨胀-压缩波系，每个波系的长度约为0.8 m，波系的强度逐次降低；当飞行高度为5 km时，喷管后方羽流区域形成了四个明显的膨胀-压缩波系，波系的膨胀半径有所增大，波系的长度有所增大，约为1.4 m，波系的强度逐次降低，羽流场边界和羽流核心区的长度都有所增加；当飞行高度为15 km时，喷管后方羽流区域有两个明显的膨胀-压缩波系，第三个波系未能充分发展，波系的膨胀半径和波的长度明显增大，羽流场边界和羽流核心区的长度较地面、中低空状态时明显增加。

图5、图6和图7分别为Ma=0.9、攻角α=0°，飞行高度H=0为0 km、5 km、15 km时羽流场轴线上马赫数、温度和压强的分布对比。

由图5～7可以得出，由于膨胀-压缩波系的存在，在不同飞行高度下轴线上的马赫数、温度和压强这些流场结构参数都随着膨胀-压缩波系个数、长度和强度的不同而发生着较为明显的变化。

2.3 相同飞行高度和飞行马赫数，不同攻角羽流场数值分析

图8为Ma=0.9、H=0 km，攻角α分别为0°、4°、8°、12°时的羽流场马赫数等值线图。

由图8可以得出，亚音速状态时，不同攻角下羽流流场结构分布的变化规律如下：

当攻角为0°时，燃气在喷管后方羽流区域形成了五个明显的膨胀-压缩波系，每个波系的长度约为0.8 m，波系的强度逐次降低；当攻角为4°时，喷管后方羽流区域仍然有五个明显的膨胀-压缩波系，波系的膨胀半径、长度和强度的变化不大，可以明显的看出膨胀-压缩波系形成上扬趋势，不再具有轴对称性，羽流核心区和羽流场边界偏离轴线；当攻角为8°时，喷管后方羽流区域仍然有五个明显的膨胀-压缩波系，波系的膨胀半径、长度仍有所减小，膨胀-压缩波系偏离轴线上扬角度加大，羽流核心区长度和羽流场边界减小；当攻角为12°时，喷管后方羽流区域仍然有五个膨胀-压缩波系，波系的膨胀半径、长度继续减小，膨胀-压缩波系较攻角为0°时与轴线形成明显的上扬角度，羽流核心区长度继续减小，在第四个膨胀-压缩波系后，羽流核心区和羽流场边界完全偏离轴线。

图9为Ma=0.9、H=0 km，攻角α分别为0°、4°、8°、12°，与喷管出口的距离为0 m、2 m、4 m、6 m时的羽流场截面温度云图。

由图9可以得出，亚音速状态时，不同攻角情况下，与喷管出口距离不同的羽流场截面温度分布的变化规律如下：

喷管出口截面，攻角为0°时温度为轴对称分布，随着攻角的增大，截面上温度的分布产生变化，攻角为12°时，可以明显看到温度分布不再具有轴对称性；距离喷管出口2 m处的截面，攻角为0°时温度仍为轴对称分布，随着攻角的增大，可以明显看到温度分布偏离中心，向Z轴方向移动，产生上扬；距离喷管出口4 m处的截面，攻角为0°时温度仍为轴对称分布，随着攻角的增大，可以明显看到温度分布向Z轴方向移动上扬的距离增大；距离喷管出口6 m处的截面，攻角为0°时温度仍为轴对称分布，由于羽流场膨胀-压缩波系不在形成，此时随着攻角的增大，温度分布上扬程度较羽流场过渡段更加明显，上扬程度加强，攻角为12°时，可以明显看到温度分布完全离开截面中心。通过对亚音速状态，不同攻角情况下，与喷管出口距离不同截面上温度分布的观察，我们进一步验证了攻角能够使膨胀-压缩波系、羽流核心区偏离轴线，产生上扬，随着攻角的增大，上扬的角度增大，羽流核心区后端会完全偏离轴线。