本次研究的对象如图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

表1 仿真算例
Tab.1 The simulation cases

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

(1)

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

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

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

CH₄+2O₂→CO₂+2H₂O(2)

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

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

p=ρRT(3)

式中:T为温度; R为气体常数。

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

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

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

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

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

表2 边界条件
Tab.2 Boundary condition

其中表示换热量; p表示压力; 表示质量流率。燃气混合物各组分的比热采用Piecewise polynomial计算,热传导率和粘性选择Kinetic-theory计算。

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

以甲烷喷孔位于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]。

图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

对针栓喷孔处于不同位置算例进行燃烧仿真,得到图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结果总结
Tab.3 Summary of Case 1 and Case 2

不同针栓喷孔直径的温度和混合比仿真结果如图 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

表4 Case 1与Case 3结果总结
Tab.4 Summary of Case 1 and Case 3

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

对不同环缝宽度(H_F)进行燃烧仿真,得到图 13的温度分布,可得出Case 4 和5的火焰锥角要大于Case 1。采用以上相同的数据处理方法,得到如图 14所示曲线,从曲线变化可知,Case 4和5燃烧室带状高温区域要比Case 1的宽,靠近壁面和对称轴处的温度Case 1要高于Case 4 和Case 5。通过计算Case 4和5平均温度分别为2 090 K和2 043 K,平均室压Case 4为0.906 MPa,Case 5为0.905 MPa。

图 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结果总结
Tab.5 Summary of Case 1,Case 4 and Case 5

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

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

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

图 17 理想混合区模型
Fig.17 Model of ideal mixing zone