气氧/气甲烷针栓发动机燃烧过程数值仿真研究

2019-01-18李清廉常一冰

火箭推进 2018年6期

周康，李清廉，成鹏，常一冰

(国防科技大学,湖南长沙 410073)

0 引言

相对传统喷注器，针栓喷注器结构简单、可靠性高且燃烧较为稳定，在应用上具有独特的优势。1961年，首个针栓喷注器成功应用于MIRA 500火箭发动机上，推力变化为111.2～ 2 224 N。此后，针栓喷注器又先后应用到阿波罗登月下降级发动机、姿轨控发动机TR 312以及Space X的梅林1 D发动机[1-2]上，足见其应用十分广泛。

目前针栓喷注器雾化研究大部分集中于对喷雾锥角的研究。对于径向环缝式针栓喷嘴，Son等[3]发现喷雾锥角随径向/轴向动量比的增大而增大，并总结出了与气液动量比TMR和We(韦伯数)有关的气液针栓喷嘴半锥角的指数关系式，数值仿真工作[4]也得到了同样的规律。方昕昕等[5]发现气液比的增加会减小喷雾锥角，径向环缝夹角、无量纲跳跃距离、液氧环缝厚度的增加均会增加喷雾锥角。成鹏[6]也研究了液液针栓喷嘴的喷雾锥角随工况的变化，并通过理论分析提出了一个与径向/轴向动量比TMR有关的半喷雾锥角的预测公式。数值模拟的结果还表明轴向液膜与径向液膜相互撞击在径向环缝出口位置形成局部高压区，该局部高压区促进了喷雾锥角增大。

但是国内外在针栓喷注器的燃烧性能与燃烧过程基础研究方面还比较有限。针栓发动机还存在燃烧效率较低的问题。国外普渡大学针对667.5 N的针栓小推力发动机，研究了总动量比、孔径比以及针栓长度与针栓直径之比对燃烧性能的影响，但并没有得出燃烧性能与这些参数之间明确的相互作用规律[7-8]。2016年日本东京大学的Sakaki等[9-10]以LOX/Ethanol为推进剂，进行了平面针栓(开窗观测)和轴对称针栓的对比试验，发现两种情况下推力室燃烧效率均对TMR呈现反相关的关系，故而Sakaki认为可以用便于使用光学观测等手段的平面针栓来进行观察针栓燃烧过程的基础研究。仿真方面，2017年国防科技大学方昕昕的进一步研究表明：跳跃距离Ls与针栓直径之比为1时，特征速度效率最高(约0.96)。方昕昕[5]和韩国航天大学的Min Son等[11-12]研究表明了径向环缝尺寸影响燃烧效率。方昕昕[5]的研究还表明燃烧室特征长度越长燃烧效率越高，这一点与常规喷注器发动机的结论是一致的。总的来说，对燃烧效率的基础研究还不够深入，需进一步明确结构工况参数对燃烧效率的影响规律。

在针栓发动机的燃烧过程基础性研究上，仿真方面国内学者以及韩国学者一般采用Lagrange随机轨道模型来模拟推进剂液滴的运动轨迹[13-16]，并以此作为雾化、蒸发乃至最终燃烧计算的初始条件。西北工业大学的李进贤[15]及北京航空航天大学的张连博等[14]对NTO/MMH针栓发动机燃烧室流场结构和推力室性能的影响因素进行了2D仿真研究，结果表明针栓深入燃烧室的长度对内流场结构有一定影响。2017年国防科技大学方昕昕的研究表明：跳跃距离过小时，径向喷出的液氧射流将在甲烷气流的撞击下加速破碎，并与燃烧室壁面发生碰撞并反弹，然后向燃烧室中心移动，造成回弹位置附近和燃烧室中心富氧而燃烧不充分，出现相对低温区。径向射流环缝大小控制着液氧流速大小，开口愈大，则液氧流速愈小，对维持燃烧稳定性起着非常重要的针栓尖部下方的回流区愈小，使得其对推进剂雾化混合的促进作用减弱。Min Son等[11-12]的研究也给出了类似结论，研究还表明相同结构下，喷注质量流率减小，火焰张角随之减小，针栓尖部下方回流区变小。总而言之，现有资料缺乏燃烧流场深入的研究，应从基本流场分布和火焰结构方面来解释结构工况参数对燃烧过程的作用规律。

纵观国内外关于针栓发动机的燃烧性能和燃烧过程基础研究还是比较有限的。跳跃距离和动量比等对燃烧过程包括流场结构、温度分布和火焰结构等的影响机理缺乏更深层次的研究，因此还需要继续开展仿真研究来明确燃烧性能和燃烧过程与针栓喷注器关键设计参数之间的规律。

1 物理模型与计算方法

1.1 物理模型与研究问题

本次研究的对象如图1所示的二维发动机结构，发动机总流量为152 g/s，其中甲烷流量为36.2 g/s。发动机喉部半径为9.55 mm，圆筒段半径为22.5 mm，圆筒段长度为48.6 mm。In1(甲烷喷注口)喷进燃烧室的是甲烷，In2(氧气喷注口)喷进燃烧室的是氧气，箭头方向表示燃料和氧化剂的流动方向。工作过程中，两股气流在燃烧室内混合燃烧，发生化学反应，产生高温高压气体，再经过喷管加速之后，以高速气流排出，在这次仿真中暂时不考虑推力室壁面传热对内流场和温度场的影响。

研究采用商业软件FLUENT进行数值仿真，研究的问题在于气氧气甲烷针栓发动机的燃烧特性，考察推进剂发生化学反应后燃烧室的内流场结构和温度场分布，研究参数主要有不同针栓喷孔直径、位置以及环缝宽度，研究的孔的位置有两个，分别为位置A和B，如图1所示。对应仿真算例如表1所示，对比算例1和2，可分析跳跃距离对内流场结构和温度场分布的影响，对比算例1，3，4和5可分析动量比对内流场结构和温度场分布的影响。

图1 仿真模型的结构Fig. 1 Configuration of simulation model

算例序号孔的位置孔径/mm环缝宽度/mm跳跃距离/mm径向与轴向动量比备注Case 1B位置0.51.311.40.66Case 2A位置0.51.370.66动量比不变Case 3B位置1.11.311.40.47变径向动量Case 4B位置0.51.511.40.81Case 5B位置0.5211.41.15变轴向动量

1.2 控制方程

处理含有化学反应的传热和流动问题需要求解质量、动量、能量守恒和化学组分输运的联立方程组，为更清楚地表达研究所涉及的控制方程及其物理意义，这里采用计算机程序计算时所采用的通用方程[17]：

(1)

式中：ρ为密度；v为速度向量;φ为待求变量;г为变量φ的扩散系数;S为源项。由该式可看出通用微分方程一般包含4项，即非稳态项、对流项、扩散项和源项，本问题研究燃烧室内定常流场，所以非稳态项直接忽略。

动量方程包括x和y两个方向的方程，源项中的两项分别代表单位质量流体的压力差和体积力。

对于化学反应，这里采用最简单的一步总包反应：

CH4+2O2→CO2+2H2O

(2)

反应共有4种组分，由于质量分数之和为1，故只需求解3个组分输运方程，能量方程是以温度T为变量来表达。

另由于流动是可压缩流动，故需要补充状态方程，假定流体满足理想气体状态方程：

p=ρRT

(3)

式中：T为温度；R为气体常数。

1.3 计算模型

这里主要说明湍流模型、化学反应模型以及湍流和化学反应相互作用的模型。

湍流模型采用Standardk-εmodel，其可适用于航空发动机和火箭发动机燃烧室内的流动计算[18]。

化学反应模型选择组分输运模型(Species Transport)。本问题研究的是气气燃烧，所以选择了组分输运模型，该模型应用广泛，能够满足本问题的仿真要求[19]，且计算成本也较小。

对于湍流和化学反应的相互作用模型，选用涡耗散模型(Eddy-Dissipation Model)，只要湍流出现(k/ε>0)，燃烧即可进行，不需要点火源来启动燃烧。这通常对于非预混燃烧是可行的。在涡耗散模型中，每个反应都有同样的湍流速率，因而模型能用于单步或是双步整体反应，并适用于高雷诺数燃烧过程[20]。该研究中采用一步总包反应，且经计算本问题的雷诺数达105量级，所以综合来看选择涡耗散模型是较为合适的。

1.4 边界条件与物性参数

仿真模型边界条件设置如表2所示。

表2 边界条件

1.5 算法使用

综合考虑各种选项的精度、计算开销、数值稳定性等，本文采用SIMPLE算法，空间离散上除湍动能和湍流耗散率采用一阶迎风各式外，其余各变量均采用二阶迎风各式，压力采用二阶格式，梯度采用最小二乘法计算。在本次研究中，氧气与甲烷的混合比为3.2，在该混合比下，根据蔡震宇[21]研究的不同室压下气氧/气甲烷理论燃烧温度，最高反应温度不超过3 470 K，超出的温度无研究意义，因此考虑算法中加入限制器，并且研究关注的是温度分布规律，而不是最高温度，这是合理化的处理。

2 网格无关性及计算方法验证

基准算例的网格无关性验证采用3套网格，网格数量级分别为6万、10万和20万。仿真主要关注燃烧室内的温度场和内流场，因而网格无关性验证的参数选择温度和压力。压力和温度的计算都是基于Area-weighted Average算法求值，原理为将所有网格面上的物理量φi与对应网格面的面积Ai乘积之和除以总面积A，即可求的对应物理量的平均值，公式如式(4)所示：

(4)

在轴向上选取某点处(该次选取x=63 mm)径向直线，导出数据进行网格无关性验证，结果如图2所示。对于不同数量网格计算的最大压力误差为1.04%，在温度方面，粗网格与密网格误差为3.01%，中度密网格与密网格基本符合，最大误差仅为0.85%，因此综合考虑精确度与计算成本使用10万数量级的中度密网格。

图2 网格无关性验证结果Fig. 2 Result of grid independent verification

对于计算方法的验证，采用高玉闪[22]的同轴剪切式气氧/气甲烷单喷嘴物理模型，计算模型和本研究的保持一致，边界条件设置和试验参数设置一致。由于研究关注的燃烧特性与温度和压力相关，而试验结果给出的是压力结果，所以此次模型验证指标为室压。通过仿真得出的室压结果为2.75 MPa，试验结果为2.56 MPa，误差为7.4%。论文关注温度分布规律的差别，在误差允许范围内(小于10%)，用此计算方法来研究气氧/气甲烷针栓发动机的燃烧特性可行的。

3 结果讨论与分析

3.1 基准算例结果与分析

以甲烷喷孔位于B位置，喷孔直径d=0.5 mm为基准算例，对其仿真结果进行分析，分析压力和组分分布云图。

此前，定义初始燃烧反应带状区域氧气浓度梯度为0.001(即氧气浓度几乎处于不变，如公式(5)所示)时的界面线与轴线方向形成的锐角为火焰锥角。燃烧效率使用特征速度效率表征。

(5)

燃烧仿真得出的基准算例燃烧室压力分布如图3所示，其平均室压为0.922 MPa。甲烷和氧气的组分分布云图如图4所示，甲烷和氧气的混合主要在于两组分的分界面之间，形成了一个混合层，化学反应也主要发生在这个混和层内。

图3 压力分布云图Fig. 3 Contour of pressure distribution

图4 组分分布云图Fig. 4 Contour of species distribution

图5的流线分布图展示了针栓喷注器存在两组明显的回流区，并且角回流区和针栓头部的回流区会影响其燃烧特性[8]。

联创汽车电子有限公司（DIAS）由上海汽车集团股份有限公司和中联汽车电子有限公司联合投资，主要从事汽车电子产品的研发、制造、销售和技术服务，致力于成为国内领先、国际一流的汽车电子系统制造商；DIAS发展迅速，已成为国内主要汽车制造公司（如上汽、上海通用和奇瑞）的供应商。

图5 流线分布云图Fig. 5 Contour of streamline distribution

定量分析温度场(如图6所示)的变化，在轴向方向上，截取x=50 mm和x=100 mm两条径向直线上温度变化(如图7所示)，温度随Y向距离增加呈现出先增加后减少的趋势。这是因为反应主要集中在甲烷与氧气混合的带状区，因此反应释热的高温区在氧化剂和燃料之间混合的带状区。基准算例的火焰锥角为54°，如图6所示，燃烧效率的结果为87.8%，平均排气速度为2 122 m/s。

图6 温度分布云图Fig. 6 Contour of temperature distribution

图7 不同位置温度随径向距离分布Fig. 7 Distribution of temperature with different radial distances

3.2 针栓喷孔位置对温度场分布的影响

对针栓喷孔处于不同位置算例进行燃烧仿真，得到图8中结果。定义温度处于3 400 K与3 500 K之间的带状区域为高温区，对应的温度变化带状区域的宽度为带宽。对温度仿真结果进行定量分析，在Case 1和Case 2的温度云图中，取x= 50 mm如和x=100 mm径向直线上的温度进行分析，其结果见图9所示。A位置(Case 2)燃烧室带状高温区域要比B位置(Case 1)的宽，但火焰锥角没有随着喷孔位置改变而改变，通过计算Case 1燃烧室内平均温度为2 179.5 K，Case 2为2 398.5 K。

图8 不同孔位置的温度云图Fig. 8 Contour of temperature with different positions of pintle orifice

图9 不同孔位置的温度随径向距离分布Fig. 9 Distribution of temperature with different positions of pintle orifice

对于以上结果，分析原因在于不同的喷注位置(Case 1和Case 2)，实质上改变了跳跃距离，Case 1的跳跃距离为11.4 mm，Case 2的跳跃距离为7 mm。跳跃距离减小，导致Case 2中氧气与甲烷混合界面带扩大，燃烧反应释热带相对Case 1宽，稳态状态下，靠近壁面的温度和对称轴的温度比Case 1温度高。同时跳跃距离减小，使得氧气甲烷在燃烧室内滞留时间比Case 1长，等效增加了特征长度，混合燃烧更加充分，导致Case 2 的平均温度要高于Case 1，室压同样如此，Case 2平均室压为0.97 MPa，Case 1为0.922 MPa，这表明Case 2的燃烧效率高于Case 1，即跳跃距离减小，燃烧效率增大，为92.3%，平均排气速度为2 235 m/s。但是跳跃距离没有改变火焰锥角。对两个算例燃烧结果总结如表3所示。

表3 Case 1与Case 2结果总结

3.3 针栓喷孔直径对温度场分布的影响

不同针栓喷孔直径的温度和混合比仿真结果如图10和图11所示，从温度云图可看出，火焰锥角随喷孔直径增加而减小。对温度进行分析，得到如图12所示变化曲线，从曲线中可知，Case 1燃烧室带状高温区与Case 3的带状高温区宽度接近，接近壁面处Case 3温度高于Case 1，而靠近对称轴附近刚好相反。并通过计算得出Case 3平均温度为2 236 K，此时Case 3平均室压为0.925 MPa，与Case 1相近。

对此进行分析，喷孔直径增大，径向动量变小，使得动量比减小，Case 1的动量比为0.66，Case 3为0.47。动量比减小，火焰锥角减小，如图 10所示，所以实质上火焰锥角随动量比减小而变小。仿真混合比为3.2，动量比减小，径向气流穿透轴向气流能力减弱，相对的轴向气流增强，扩散到壁面处氧气增加，从图11可知，3和4 之间的局部混合比区域面积Case 3大于Case 1，即动量比减小导致靠近壁面处氧气与甲烷混合反应更加充分，而总体性能燃烧效率与Case 1接近，结果为88.1%，平均排气速度为2 152 m/s。两者燃烧性能总结如表4所示。

图10 不同孔直径的温度云图Fig. 10 Contour of temperature with different pintle orifice diameters

此外，从图11可看出在带状高温区域下方的混合情况很差，这是因为在二维情况下，氧气穿透甲烷气膜进入到带状高温区域下方的量很少，导致混合变差。改善的方法在于减小径向气膜喷射速度，减小径向动量，增加轴向气膜的穿透能力，使得进入到带状高温区域下方氧气的量增加，提高氧气与甲烷的掺混，从而改善其燃烧性能，Case 1和Case 3的仿真结果验证了该方法的可行性。

图11 不同孔直径的混合比Fig. 11 Contour of mixture ratio with different pintle orifice diameters

算例燃烧效率/%排气速度/(m·s-1)火焰锥角/(°)Case 187.82 12254Case 388.12 15235

图12 不同位置不同孔径温度随径向距离分布Fig. 12 Distribution of temperature with different positions of pintle orifice diameter

3.4 针栓环缝宽度对温度场分布的影响

图13 不同环缝宽度的温度分布Fig. 13 Contour of temperature with different annular gaps

图14 不同位置不同环缝宽度温度随径向距离分布Fig. 14 Distribution of temperature with different annular gaps

图15 不同环缝宽度流线图Fig. 15 Streamline with different annular gaps

分析以上结果，环缝宽度变大，使得轴向动量减小，动量比增加，混合带状区域向燃烧室头部倾斜，导致燃烧时Case 4和Case 5火焰锥角变大。动量比增加，径向气流撞击壁面，使得扩散到靠近壁面处氧气减少，氧气与甲烷混合变差，反应释热减少，温度降低，如图14曲线分布，并且动量比增加，靠近针栓头部的回流区轴向尺寸减小，如图15所示，对靠近喉部气体卷吸作用减弱，气体与高温区交换热量减小，轴向温度降低。环缝宽度增加，导致动量比增加，氧气与甲烷混合燃烧变差，使得燃烧室内平均温度降低，平均室压降低，燃烧效率降低，Case 4和Case 5的燃烧效率分别为86.3%和86.1%，平均排气速度分别为2 055 m/s和2 036 m/s。以上燃烧结果可总结如表5所示。

表5 Case 1，Case 4 和Case 5结果总结

3.5 火焰锥角对燃烧效率的影响

在这次仿真中发现喷孔直径、位置以及环缝宽度均会对燃烧效率产生影响，而其中火焰锥角也会发生相应变化，因此燃烧效率与火焰锥角其中规律可以通过曲线图来说明，如图16所示。燃烧效率随着火焰锥角增加而减小。

图16 不同火焰锥角下的燃烧效率Fig. 16 Combustion efficiency with different flame angles

对此现象可以从混合角度解释，本文理论混合比为3.2，因此局部混合比处于3到4之间可认为是理想混合范围，因而可从3.3节中图11可以简化建立一个理想混合区模型，如图17所示。理想混合主要有两个，一个靠近燃烧室头部，另一个靠近燃烧室壁。在前面分析中可知，火焰锥角的变化是由动量比变化引起的，动量比增加，火焰锥角增加，而此时径向射流由于动量相对轴向增加，使得射流撞击壁面，而使得靠近壁面氧气减小，而甲烷浓度没有变化，因而混合变差，靠近壁面的理想混合区域被压缩变小，燃烧不充分，燃烧效率也会相应降低。因而提高燃烧效率可以通过减小火焰锥角实现，但火焰锥角不能无限减小，推测存在一个变化范围，这需要通过进一步研究来确定。