三种典型固体推进剂排气羽流特性参数的模拟计算

李　猛，孙　美，赵凤起，轩春雷，仪建华，杨燕京，孙志华，王长健，许　毅

(西安近代化学研究所燃烧与爆炸技术重点实验室，陕西 西安 710065)

引　言

固体推进剂排气羽流是一种超音速排出喷管的气相和凝聚相粒子混合物，并伴随有燃烧产物的二次燃烧、湍流、电子激发和电离等，在喷管出口处下游常形成羽毛状的发光火焰流场[1-3]。排气羽流的气动力结构不仅导致速度和压力的非连续性，还影响气流与周围空气的混合过程[3-4]。如机载导弹发射时排气羽流所形成的压力扰动使得进气道增压，产生的温度与压力畸变引起压力机失速或停车；舰载导弹对排气羽流的防护要求较高，由于舰面范围有限，其上层建筑和雷达外接设备等都可能受到羽流的危害；排气羽流对飞行中的火箭导弹的空气动力特性和飞行稳定性等都有较大影响。流速和温度是排气羽流重要的特性参数，国外学者对其开展了数值模拟及实验验证，并将模拟结果和实验数据进行了对比分析[5-7]；国内近年在固体推进剂排气羽流数值模拟方面也开展了一些研究[8-9]，建立了排气羽流特性测试表征手段[10-12]，但对于排气羽流场数值模拟的准确性和可靠性未进行验证，使得模拟计算与工程实际严重脱节。

本研究采用基于最小自由能原理的能量计算星程序[13-15]和流体仿真商业程序ANSYS-Fluent对双基(DB)推进剂、复合改性双基(CMDB)推进剂及硝酸酯增塑聚醚(NEPE)推进剂进行了流场模拟计算，获得了3种推进剂装药排气羽流特性参数，采用红外热像仪和TDLAS流速测试系统对其进行了实验验证，将模拟结果与实验结果进行对比，为数值模拟技术用于工程实践及数值模型修正提供参考。

1　数值模拟

1.1　热力学计算

根据推进剂燃烧产物组分没有自发的变化趋势时体系自由能最小的原理建立自由能函数方程，联立质量和能量守恒方程，获得给定压强下燃烧产物组分、绝热燃烧温度等热力学参数[13-15]。

在燃烧室压强为7MPa、喷管出口压强为0.1MPa的假设条件下，应用能量计算星程序计算获得了双基(DB)推进剂、复合改性双基(CMDB)推进剂及硝酸酯增塑聚醚(NEPE)推进剂3种推进剂配方的燃烧温度、燃烧产物组分及含量等热力学参数，并将获得的热力学参数作为羽流场计算的初始条件。

1.2　计算区域及网格结构

对于推进剂装药排气羽流来说，流场中各处温度、速度等特性分布都具有轴对称性。为简化前处理、提高计算速度及满足计算精度，采用二维轴对称计算区域对羽流场进行计算，图1为排气羽流计算区域示意图。

边界条件：喷管入口处设置为压力入口边界，其总压为恒定值；喷管壁面为壁面边界条件，假设流体和壁面之间没有热量传递，固体壁面静止不动，且为无滑移、绝热壁面；喷管及羽流轴线为轴对称边界条件；其余边界为压力出口边界，出口处的气体参数按照试验过程中的大气条件给出。由推进剂羽流特性参数实验可知，其采用Ф50mm标准固体火箭发动机，喷管几何型面均为收敛-扩张的拉瓦尔型，不同的是喷管喉部半径及出口半径和扩张段长度。在Gambit中编写并固化了Journal文件对计算区域进行建模和对称结构化网格划分，在壁面处和喷管出口后缘速度剪切层进行加密，网格数量约为两万。因此只用在Journal文件中改变喷喉半径、出口半径和扩张段长度，即可实现3种推进剂装药羽流特性计算的网格自动划分。图2为采用上述方法获得的双基推进剂装药羽流计算网格。

1.3　排气羽流场计算

为获得较高的计算精度，采用二阶迎风格式的有限体积法对羽流流场的气凝两相多组分控制方程进行离散，并采用密度基求解器来模拟高马赫数流场，湍流模型选用RNGk—ε模型[16]，然后采用有限反应速率模型和DPM模型分别研究排气羽流多组分的复燃和凝聚相颗粒的运动规律。

由于高速运动的气流对凝聚相Al2O3颗粒的剪切作用，排气羽流中颗粒的直径通常小于10～15μm，且是多种直径组成的颗粒群，考虑到工程实际情况及保证相应计算精度的同时尽可能简化，在数值计算过程中认为铝颗粒直径分布在1～10μm之间，使用Rosin-Rammler分布规律获得羽流中凝聚相颗粒的粒度分布情况，即：

(1)

式中：F为分布函数；d为粒径尺寸；D43为平均粒径；n为指数系数。

对于某些包含原子及自由基的反应，需要考虑第三种物质以接受过剩的反应能量，从而使反应得以稳定进行。为更真实地反映排气羽流二次燃烧所带来的温度场变化，在分析上述3种固体推进剂的燃烧产物组分特点及含量的基础上，对DB推进剂和CMDB推进剂，主要采用C /H /O 组分系统9种组分、10个基元反应的有限速率化学反应模型来模拟羽流中的二次燃烧现象；由于NEPE推进剂的燃烧产物包含Cl元素，故采用C /H /O /Cl组分系统11种组分，14个方程的有限速率化学反应模型来研究羽流二次燃烧。所用反应机理数据如表1所示。

表1　3种典型推进剂排气羽流中的化学反应

为节省计算资源，首先计算单组分的气相流场，待计算收敛且结果符合理论解后，在现有流场结果的基础上开启组分反应模型，设置各组分在各个边界处的含量分布和组分之间的反应模型，从而进行组分二次燃烧的模拟。最后再开启DPM模型，设置颗粒的材料、速度、流量等属性，并加载UDF来控制颗粒的运动、二次燃烧及气相、颗粒相之间的耦合来完成最终的模拟。

2　实　验

2.1　推进剂及装药

选取DB、CMDB、NEPE 3种固体推进剂，其中DB推进剂和CMDB推进剂采用螺旋压伸工艺制备，NEPE推进剂采用淤浆浇铸工艺制备。推进剂配方见表2。

表2　3种典型推进剂的配方

根据推进剂装药排气羽流检测需要，选定Ф50mm标准固体火箭发动机，制备3种推进剂单孔管状装药，其中外径为45mm，内径为8mm；单端包覆装药。同时按GJB770A方法706.1靶线法要求制备药条试样，并测试燃速，由3种推进剂在7MPa下的燃速、推进剂密度以及计算获得的特征速度，根据内弹道理论设计出3种推进剂装药地面静止发动机实验所需要的喷管尺寸，见表3。

表3　实验喷管构型尺寸

注：R为喷喉直径；L为喷喉长度；r为入口半径；θ1、θ2分别为收敛角和扩张角；η为扩张比

2.2　实验装置

羽流特性参数红外辐射温度采用德国IRCAM327型红外热像仪测试[17-18]。羽流流速采用可调谐半导体激光吸收光谱技术(TDLAS)测量[12]，即利用半导体激光器的窄线宽和可调谐特性对气体分子的吸收谱线进行快速扫描和测量，实时反演出气流速度。

2.3　实验原理及方法

2.3.1羽流温度的测量

利用红外热像仪对羽流流场进行实时温度测量，根据羽流发射率计算羽流温度场的红外辐射温度平均值和最高值的分布[17-18]。以火箭发动机羽流场作为红外辐射源，在一定位置上布设红外热像仪和比色计对辐射温度进行测量，利用比色温度对被测目标发射率不依赖的特性，计算获得了发动机羽流场发射率与时间的关系，对红外热像仪测得的温度分布进行实时修正，使得测量结果的准确度大幅度提高。接通点火电源前触发红外热像仪采集数据，在发动机工作完后停止采集数据。红外热像仪放置位置如图3所示。

2.3.2羽流速度的测量

利用TDLAS进行流场速度的非接触式测量，其实质是利用多普勒效应，利用频移来实现流场速度测量[10-12]。将TDLAS流速测量装置放置于固体火箭发动机喷口后方，TDLAS羽流流速测试系统测量点(激光交叉点)距发动机喷口30cm，以喷口轴心为原点，则TDLAS的激光发射点与接收点的坐标分别为(0.05，0.25)、(0.55，0.25)、(0.05，-0.25)、(0.55，-0.25)。在接通点火电源使发动机点火时，触发TDLAS流速测量系统采集数据，在发动机工作完毕后停止采集数据。系统示意图如图4所示。

3　实验与数值模拟结果对比分析

3.1　排气羽流温度结果的对比

采用红外热像仪在实验室条件下对3种推进剂装药的羽流红外辐射温度进行测试，结果见图5。

由图5可以看出，CMDB推进剂羽流红外辐射面积最大，NEPE推进剂羽流红外辐射面积次之，且辐射温度最高，DB推进剂羽流红外辐射面积最小，红外辐射温度最低。

表4列出了3种推进剂羽流红外辐射温度测试结果的最大值和平均值、能量星程序计算的燃烧室温度及模拟计算的羽流温度。由表4可知，3种推进剂红外辐射温度最大值和平均值的大小顺序一致，NEPE推进剂最高，CMDB推进剂次之，DB推进剂最小。通过与能量星程序计算获得的燃烧室和完全膨胀后羽流的温度值进行对比，发现3种推进剂羽流温度值之间大小顺序一致，能量星程序计算结果可以对羽流温度进行定性判断。

表4　3种典型推进剂羽流红外辐射温度数据

注：t1和t2分别为红外辐射温度的测试值和计算值。

图6为3种推进剂装药排气羽流数值模拟温度场云图。

由图6的温度云图可以看出，由于出口压强高于环境压强，燃气出喷管即开始膨胀解压，温度降低，进而截断斜激波造成温度间断升高，激波遇到混合层边界反射成膨胀波，从而导致新一轮的膨胀开始，如此往复，马赫盘和温度间断不断出现，直到射流与外界环境压力相同。同时，二次燃烧主要发生在排气羽流混合层内，并且随着排气羽流距离的增加而逐渐与主流混合，由于混合层二次燃烧的大量放热，混合层的温度比核心区温度要高几百摄氏度。由于DB推进剂的氧平衡较低，使得DB推进剂二次燃烧最剧烈，CMDB推进剂次之，NEPE推进剂二次燃烧相对最弱。图7为3种推进剂装药发动机轴线上的温度分布曲线。

由图7可看出，NEPE推进剂羽流温度最高，CMDB推进剂次之，DB推进剂最小，流场计算获得的羽流温度与红外热像仪测试结果一致。由于排气羽流的膨胀压缩波系，使得排气羽流温度出现振荡，但经过马赫盘后，燃气的温度突然升高。由于NEPE推进剂的能量较高，燃烧温度高，喷管喷出燃气的温度也相应较高，同时由于温度的升高，在混合层处化学反应速率加快，复燃产生的热量使混合层局部温度更高，从而使此区域等温线分布更密。CMDB推进剂次之，DB推进剂能量最小，推进剂燃烧温度也最小。

综上所述，羽流温度作为排气羽流特性重要参数之一，本实验测试采用的红外热像仪是非接触式测量，由于辐射率等因素的影响，测试值比真实值偏小；所采用的羽流场数值模型是基于能量守恒定律的，相对于外界是绝热状态，与实际排气羽流温度有差别，计算值比真实值要偏高。

3.2　排气羽流流速结果对比

采用TDLAS对3种推进剂装药的羽流流速进行测试，结果见图8。

由图8可知，羽流流速测试曲线与推进剂装药工作过程p—t曲线一致，经分析和数据处理，获得3种推进剂装药羽流流速平均值分别为879.0、1103.0、1457.1m/s，总体来看，由于NEPE推进剂能量高、CMDB推进剂次之，DB推进剂能量最小，由一维等熵流动假设可知，羽流流速与推进剂能量性能大小顺序一致，本研究测试结果可靠。在羽流场计算数据中，选取TDLAS的激光发射点、接收点及激光交叉点坐标处的流速数据，经数学处理获得羽流流速模拟值为1002.38、1279.01、1243.16m/s，模拟值与计算值之间偏差为14.04%、15.9%、-14.7%。

图9为3种推进剂装药排气羽流速度场云图。

由图9可知，气流在喷喉处被加速到音速，随后流通截面增大，气流继续加速到超音速，由于喷管出口处的压强小于外界环境压强，处于高度过膨胀状态，随着气流经过喷管出口，伴随着压强升高、激波产生，激波形状是从喷口边缘发出的收缩截锥形，截锥顶部为一盘形正激波即马赫盘，经过马赫盘后，速度迅速降低。此后，羽流再次经过一个膨胀波系，速度升高，直到下一个马赫盘出现，如此反复，羽流核心区和边界层最终混合均匀。

图10为3种推进剂在发动机轴线上的羽流流速。

由图10可看出，NEPE推进剂羽流流速最高，CMDB次之，DB推进剂最小。流场计算获得的羽流流速与TDLAS羽流流速测试结果一致。单就一种推进剂羽流流速轴线分布情况来看，羽流流速数据经过多次振荡后达到平衡状态，随着远离喷口速度逐渐衰减到零，其中NEPE推进剂羽流流速振荡次数最多，其次为CMDB推进剂，DB推进剂振荡次数最少，这与排气羽流马赫盘的形成分析一致。

3种推进剂羽流流速模拟结果与实验测试结果大小顺序一致，分析其原因主要是：(1)TDLAS实验技术选取的是1392nm附近的H2O吸收谱线，并且在均匀气相风洞中进行标定，而本研究的推进剂排气羽流是包含气、凝两相多组分、含化学反应的流场，其速度呈现明显的滞后，根据文献[11-12]分析，测量值与理论值偏差可能是流场的不均匀导致；(2)TDLAS实验测试羽流流速是在测量光学路径上的平均气流流速，由于考虑到实验装置的安全性，TDLAS装置架设位置不可能完全贴合排气羽流场，造成排气羽流场外缘很多静止的空气也被考虑到羽流流速中，导致测试结果偏低。

4　结　论

(1)通过日志文件驱动Gambit和ANSYS-Fluent程序读入Journal文件来自动生成3种推进剂羽流对称结构化网格和开展气凝两相化学反应流场模拟，采用C/H/O 9种组分、10个基元反应模拟DB、CMDB推进剂羽流二次燃烧，采用C/H/O/Cl 11种组分、14个方程来模拟NEPE推进剂羽流二次燃烧，只需输入燃气成分及含量即可获得羽流温度和速度分布。

(2)采用红外热像仪和可调谐半导体激光吸收光谱技术(TDLAS)对排气羽流特性参数羽流温度和羽流速度进行了实验验证，将模拟结果与测试结果进行对比，对产生偏差的原因进行了分析，分析结果可对模型修正提供一些途径，为实验系统的设置及实验方法的改进提供参考。

参考文献:

[1]Thorn L B, Wharton W W. Terminology and assessment methods of solid propellant rocket exhaust signature, AGARD-AR-287[R]. Loughton: AGARD, 1993.

[2]Chastenet J C, Kessel P A , Lawrence R, et al. Rocket motor plume technology, AGARD-LS-188[R]. Loughton: AGARD, 1993.

[3]李猛,王宏.火箭发动机排气羽流特征信号预示技术研究综述[J].飞航导弹,2010(3):54-57.

LI Meng, WANG Hong. A review studies on the technology for prediction of rocket motor plume signature[J]. Aerodynamic Missile Journal,2010(3):54-57.

[4]Guery J F, Chang I S, Shimada T, et al. Solid propulsion for space applications: an updated road map[J]. Acta Astronautics, 2010, 66: 201-219.

[5]Troyes J, Dubois I, Borie V, et al. Multi-phase reactive numerical simulations of a model solid rocket motor exhaust jet[C]∥42th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. New York: AIAA, 2006:1-13.

[6]Wang X Y ，Yuko J R.Thermal analysis on plume heating of the main engine on the crew exploration vehicle service module, NASA/TM-215049[R].Washington D C:NASA, 2007.

[7]Gueyffier D, Fromentin-Denoziere B, Simon J, et al. Numerical simulation of ionized rocket plumes[C]∥44th AIAA Plasmadynamics and Laser Conference. New York: AIAA, 2013:1-13.

[8]李猛, 张晓宏, 孙美, 等. 改性双基推进剂两相化学反应羽流特性研究[J]. 弹箭与制导学报, 2012, 32(1): 123-126.

LI Meng, ZHANG Xiao-hong, SUN Mei, et al. The study on plume characteristic of two phase chemical reaction of modified double-base propellant[J]. Journal of Projectiles, Rockets, Missiles Guidance, 2012, 32(1): 123-126.

[9]李峥, 向红军, 张小英. 复合推进剂固体火箭发动机喷流流场数值模拟[J]. 固体火箭技术, 2014, 37(1): 37-42.

LI Zheng, XIANG Hong-jun, ZHANG Xiao-ying. Numerical simulation of composite solid propellant rocket motor exhaust plume[J]. Journal of Solid Rocket Technology, 2014, 37(1): 37-42.

[10] 杨燕京, 赵凤起, 孙美, 等. 基于TDLAS技术的双基系推进剂装药排气羽流流速特性研究[J]. 火炸药学报, 2016，39(3)：70-74.

YANG Yan-jing, ZHAO Feng-qi, SUN Mei, et al. Investigations on the exhaust plume velocity characteristics of double base and CMDB propellants charge based on TDLAS technique[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2016，39(3)：70-74.

[11] 杨斌, 黄斌, 刘佩进, 等. 利用TDLAS技术评估火箭基组合循环发动机试验性能[J]. 宇航学报, 2015，36(7)：840-848.

YANG Bin, HUANG Bin, LIU Pei-jin, et al. Test performance evaluation of rocket based combined cycle engine using TDLAS technique[J]. Journal of Astronautics, 2015，36(7)：840-848.

[12] 贾良权, 刘文清, 阚瑞峰, 等. 波长调制-TDLAS 技术测量风洞中氧气流速方法研究[J]. 中国激光, 2015 (7): 304-311.

JIA Liang-quan, LIU Wen-qing, KAN Rui-feng, et al. Study on oxygen velocity measurement in wind tunnel by wavelength modulation-TDLAS technology[J]. Chinese Journal of Lasers, 2015 (7): 304-311.

[13] McBride B J, Gordon S. Computer program for calculation of chemical equilibrium compositions and applications: II users manual and program description, NASA RP-1311[R]. Washington D C:NASA, 1996.

[14] 李猛, 赵凤起, 徐司雨, 等. 三种能量计算程序在推进剂配方设计中的比较[J]. 火炸药学报, 2013, 36(3): 73-77.

LI Meng, ZHAO Feng-qi, XU Si-yu, et al. Comparison of three kinds of energy calculation programs in formulation design of solid propellants[J]. Chinese Journal of Explosives & Propellants(Huozhayao Xuebao), 2013, 36(3): 73-77.

[15] 田德余,赵凤起,刘剑洪.含能材料及相关物手册[M].北京:国防工业出版社,2011.

[16] Yakhot V , Orszag S A. Renormalization group analysis of turbulence: I. basic theory[J]. Journal of Scientific Computing, 1986,1(1):1-51.

[17] 仪建华. 火箭发动机羽流红外辐射温度动态测量方法：CN，201610846494.4[P]. 2016-09-23.

[18] 张波, 吕事桂. 基于双响应波段工作的红外热像仪测温原理与误差分析[J]. 宇航计测技术, 2010, 30(6): 11-15.

ZHANG Bo, Lü Shi-gui. Temperature measurement of the thermal infrared imager based on two response working wavebands and its error analysis[J]. Journal of Astronautic Metrology and Measurement, 2010, 30(6): 11-15.