火箭燃气射流流场的无网格方法模拟①

卓长飞，封 锋，武晓松

(南京理工大学机械工程学院，南京 210094)

0 引言

研究火箭燃气射流在工程技术领域有着重要意义，如机载火箭导弹发射时燃气射流所形成的压力扰动往往使得进气道内增压，它产生的温度与压力畸变都是引起压气机失速或发动机停车的重要因素。燃气射流的吸入与药柱燃烧后产物的再燃烧也对飞机发动机正常工作有影响。此外，发动机燃气射流对飞行火箭外部流分离和底压也有影响，可以改变它们的空气动力特性［1］。

近年来，火箭燃气射流动力学问题的数值模拟水平不断提高，使得火箭燃气射流动力学获得迅速发展。然而，在数值研究火箭燃气射流流场方面，均是建立在网格算法基础上［2－4］。实际上，计算流体力学的数值模拟方法可划分为2种:网格算法和无网格算法。网格算法分为非结构网格算法和结构网格算法，网格算法发展较早，技术比较成熟。相比之下，新兴的无网格算法在非结构网格算法打破网格结构化思维的基础上，进一步抛弃了网格化的思维，彻底地打破了结构化和网格化的思想约束，放弃了网格而不需要像网格方法那样形成网格单元。无网格算法的区域离散只涉及离散点，在点云结构上离散控制方程，从而不存在网格质量和拓扑结构限制等问题。因此，无网格算法比网格算法具有更大的几何灵活性，在处理复杂外形和复杂流动过程方面具有独特的优势。在过去的20多年中，国外计算流体力学研究者对无网格算法的格式精度和计算效率进行了大量研究［5－7］。而国内对无网格算法研究较少，但无网格算法仍在不断发展之中，如将基于网格算法的数值方法，如Roe格式、AUSM+－UP格式等迎风格式推广到无网格算法［8－9］、无网格的隐式算法研究［8］、无网格算法的精度分析［10、无网格在移动边界流场中的应用研究［11］，并逐渐将无网格算法用于工程问题［12－13］。

火箭燃气射流含有富燃气体，能与外流中的空气发生化学反应，整个流场是典型的化学非平衡流问题。为了研究采用无网格算法进行火箭燃气射流流场的可行性，本文详细给出在无网格条件下求解带化学反应的多组分Euler方程具体过程，并验证了采用无网格算法进行化学非平衡流数值模拟的可行性。最后，对不同工况火箭燃气射流流场进行了数值研究。结果表明，无网格方法模拟火箭燃气射流流场是可行的。

1 数值方法

1．1 控制方程

在笛卡尔坐标系下，微分守恒形式的二维轴对称化学非平衡流Euler方程为

式中 Q为守恒变量;F、G分别为轴向和径向无粘通量;H为轴对称源项;S为化学反应源项。

这里仅列出主要的变量:

式中 ρ为密度;u、v分别为轴向、径向速度;p为压力;E 为单位体积的总能;ci(i=1，．．．，N－1)为 i组分质量分数;ωi为i组分生成速率;N为总组分数。

1．2 点云的概念

无网格算法只需将计算域离散成一系列点，并按一定规律组织形成点云，计算域内每个点都有自己的点云结构。流动计算过程中，控制方程的离散即建立在点云基础上。计算域内某点的点云结构是由该点和周围的一些相邻点所构成，该点称为中心点，而周围相邻点称为卫星点。某一点云结构见图1，点i的点云结构由中心点i和6个卫星点组成。

图1 离散域内节点i的点云结构Fig．1 The structure for clouds of point i in discrete domin

1．3 无网格空间导数逼近原理

对控制方程进行离散，首先利用曲线拟合方法，确定每个点云结构内物理变量的空间导数。本文采用一阶线性函数进行曲线拟合，函数为

每个点云的中心点及其卫星点都满足该点云结构的解函数:

式中 k表示i的点云结构中卫星点总个数。

则存在如下矩阵关系:

解该矛盾方程组，得到该点云结构内a0、a1、a2的最小二乘解。则点云结构中i处的空间导数可表示为

式中 下标j表示i点的点云结构中卫星点。

由此可确定该点云结构内中心点与每个卫星点之间的传播系数 bij、cij。

多组分化学非平衡流Euler方程在点云i上半离散形式为

由式(6)可得到多组分化学非平衡流Euler方程的对流项空间导数为

其中，Eij表示中心点i与卫星点j之间的通量，其值由i和j的守恒量得到，如图2所示。求解该通量采用高精度高分辨率的AUSMPW+迎风格式［14］。在求解通量之前，首先对点i和j的守恒变量进行重构到中点的左、右间断处:

式中 φ表示抑制数值振荡的限制器;r表示矢径。

守恒变量的梯度则由式(6)计算得:

图2 中心点和卫星点通量重构Fig．2 The flux reconstruction between center point and satellite point

1．4 化学反应动力学模型与刚性问题

化学反应动力学模型是决定成功模拟化学非平衡流的关键之一。常规的双基固体推进剂产生的燃气主要是由 CO、H2、CO2、H2O、N2、HCl等组成，因此采用CO—H2—O2化学反应系统。本文采用8组分(CO、H2、O2、CO2、H2O、H、OH、O)12 个基元反应的 CO—H2—O2系统化学反应模型［15］，基元反应表达式和系数见表1。表1中，M代表第三碰撞体;A为指前因子，其中n为基元反应级数;b为温度指数;E为活化能。

表1 CO—H2—O2基元反应模型Table 1 Detailed reaction model for CO—H2—O2

化学非平衡流控制方程组分为流动部分和化学反应部分，两者相互耦合，并会产生刚性问题。本文采用时间算子分裂算法处理这种耦合过程，即首先计算流体流动效应，得到物理变量的过渡值。然后，继续计算，将化学反应的贡献叠加到物理变量过渡值，最终得到下一时刻体现整体效应的物理量值。其中，在计算化学反应对流场的贡献时，需把求解流动偏微分方程时采用的时间步长进一步细分，作为求解化学反应刚性常微分方程的步长。具体做法是先冻结化学反应，求解得到流场参数;然后，将化学反应看做等容吸热或放热过程，保持内能、速度参数不变，计算各组分的质量变化率;最后，迭代求解温度。

1．5 物理模型与边界条件

本文取燃气射流场轴向计算长度为弹体直径的40倍，径向计算长度弹体直径的10倍。火箭燃气射流流场计算区域中喷管出口附近空间布点如图3所示，布点总数为28 049个。在本文计算中，不考虑火箭发动机燃烧室和喷管内流动，直接在喷管出口给定燃气超声速射流条件，并保持喷管出口燃气射流的马赫数Maj=2．5、静温Tj=1 500 K、燃气组分与质量分数(如表2所示)不变。

表2 燃气主要组分与质量分数Table 2 The main mass fraction of rocket gas

图3 射流流场空间布点示意图(喷管出口附近)Fig．3 Clouds of points for rocket gas efflux field

2 数值验证

2．1 NACA0012翼型跨声速绕流

首先，针对NACA0012翼型的跨声速流场进行计算。图4给出了NACA0012翼型的空间布点局部放大图。来流条件为马赫数0．8，攻角1．25°，计算区域布点总数为6 327个。该算例为二维平面问题，计算时关闭轴对称源项，且整个流场气体为理想完全气体，不考虑化学反应。计算得到的翼型上下表面压力系数图5所示。由图5可看出，跨声速流场存在复杂的波系，上下表面的压力系数与试验值吻合较好，在激波处没有出现数值振荡，这说明本文发展的无网格方法数值模拟精度较高，具有模拟复杂波系的能力。

图4 NACA0012翼型空间布点示意图(表面局部)Fig．4 Clouds of points for NACA0012 airfoils(surface)

图5 NACA0012翼型上下表面压力系数分布Fig．5 Pressure coefficient distributions of NACA0012 airfoils

2．2 楔体诱导斜爆轰流场

计算工况:混合气体为2H2+O2+3．76N2，来流马赫数7．5，静温 293 K，静压4 000 Pa，尖劈角度为 25°，H2－O2化学动力学模型采用7组分8反应模型［16］。计算区域布点如图6所示，布点总数为15 074个。

计算得到的温度场如图7所示。高速可燃气体通过尖劈会发生偏转，形成斜激波，波后混合气体受斜激波压缩的作用，压力和温度有所提高，但不能直接点燃混合气体。压缩后的混合气体经过一定的诱导期后形成斜爆轰波。整个流场是由斜激波、横波、斜爆轰波、滑移线组成，这些波系相交于三波点。本文采用无网格算法的数值结果与文献［17］数值计算得到的温度云图完全吻合，各种波系均清晰可见，说明无网格算法能适用于具有复杂化学非平衡流场的数值模拟。

图6 楔体诱导斜爆轰流场的空间布点示意图Fig．6 Clouds of points for wedge oblique detonation

图7 斜爆轰流场温度云图与流场结构Fig．7 Temperature and fields structure of oblique detonation fields

3 计算结果与分析

3．1 欠膨胀状态下射流流场

计算工况为来流马赫数2．0，来流静压为0．4 atm，火箭喷管出口的燃气射流压力分别为 1．0、2．5、5．0 atm。当喷管出口的燃气射流压力大于外流压力时为欠膨胀状态，且随着燃气射流压力的增大，欠膨胀程度增加。不同喷管出口燃气射流压力下燃气射流流场的马赫数云图如图8所示，轴线马赫数、压力、温度如图9所示。

图8 欠膨胀状态下射流流场马赫数云图Fig．8 Mach contours of field flow at under-expanded state

由图8看出，相交射流激波、反射激波、射流边界等主要流动特征非常清晰，且由图9中轴线参数分布可看出，在轴线上相交射流激波第1次相交形成的强间断处，各参数均发生突变，具有较大的梯度变化率，这完全符合激波物理性质。这些均说明本文采用无网格方法模拟得到的激波结构非常清晰且合理，具有较强的捕捉强间断流场的能力，能较好地模拟火箭燃气射流流场。

欠膨胀状态下，燃气流出喷管后在喷管唇口附近发散成一束扇形膨胀波族，形成Prantl－Meyer流，气体将继续膨胀加速，流动参数中的马赫数升高，压力、温度降低。同时，燃气流出喷管后，会受到外流的压缩作用，而形成相交射流激波，上、下相交射流激波相交于轴线上，并发生发射形成反射激波，反射激波又在混合层发生反射。这样，在射流轴线上重复发生着膨胀－压缩过程，且强度不断减弱，这个现象可由轴线上各参数变化看出，各流动参数在一定范围内呈振荡衰减型，最后趋于稳定变化。

由图8、图9可看出，相交射流激波在轴线上的交点，随着欠膨胀程度的增加而远离喷管出口。这是由于欠膨胀程度的增加，导致燃气射流对外流的作用增强，即相交射流激波强度增强，射流激波张角增大，从而交点更靠后。由流动参数轴线分布图还可看出，随着欠膨胀程度的增加，燃气流出喷管后膨胀程度也大大提高，在相交射流激波前的马赫数更高，压力、温度更低。

射流场轴线上离开喷管一定距离后，不同喷管出口燃气射流压力对应下的轴线压力也基本一致，且均约为外流的静压。这主要是由于燃气射流与外流之间混合层的相互作用，通过剪切层的传输，导致燃气射流轴线上的静压与外流的静压保持一致，而马赫数、温度仍保持较高水平。可见，火箭燃气射流对尾流场的影响范围非常大。

图9 欠膨胀状态下轴线上主要参数分布Fig．9 Distribution of main parameter along the axis at under-expanded state

3．2 过膨胀状态下射流流场

计算工况为来流马赫数2．0，来流静压为1．0个大气压，火箭喷管出口的燃气射流压力分别为0．8、0．6、0．4 atm。当喷管出口的燃气射流压力小于外流压力时属于过膨胀状态，且随着燃气射流压力的降低，过膨胀程度将增加。不同喷管出口燃气射流压力下射流流场的马赫数云图如图10所示，轴线马赫数、压力、温度如图11所示。

图10 不同过膨胀状态下射流流场马赫数云图Fig．10 Mach contours of field flow at over-expanded state

图11 过膨胀状态下轴线上主要参数分布Fig．11 Distribution of main parameter along the axis at over-expanded state

在过膨胀状态下，主要流动特征仍然非常清晰可见。从流场结构来看，过膨胀和欠膨胀状态相比的重要区别之一是过膨胀状态下相交射流激波比较平直，而欠膨胀状态下相交射流激波呈曲线。

随着过膨胀程度的增大，相交射流激波在轴线上的交点不断靠近喷管。这是由于随着燃气射流压力降低，外流对燃气射流的“挤压”作用更明显，从而使相交射流激波在轴线上的交点更靠前。随着过膨胀程度的增大，射流激波会在轴线以及射流与外流的混合层处发生的反射次数增多，轴线上的参数也明显表现为呈振荡衰减型，最后趋于稳定。值得一提的是在燃气射流压力为0．8 atm时，由轴线参数分布看出，在相交激波交点前，燃气射流会先进行一定程度的膨胀加速，马赫数升高，压力、温度降低。这可能是由于喷管出口燃气射流与外流压力很接近，处于过膨胀和欠膨胀之间，但不是完全膨胀状态，因而产生既含欠膨胀特征，也含过膨胀特征的现象。

3．3 来流马赫数对射流流场的影响

计算工况为来流马赫数分别为 1．5、2．0、2．5 Ma，来流静压为0．4 atm，火箭喷管出口的燃气射流压力为1 atm。不同来流马赫数下射流流场的马赫数云图如图12所示，轴线马赫数、压力、温度如图13所示。

由计算结果可看出，超声速来流条件下，仅改变不同来流马赫数，对射流流场总体结构基本无影响，相交射流激波相交点位置基本相同，轴线上参数分布也基本一致。值得注意的是离开喷管一定距离后，由于射流与外流的充分混合和发展，轴线上压力会以不同程度降低，逐渐与外流保持一致。

图12 不同来流马赫数下射流流场马赫数云图Fig．12 Mach number contours of field flow variation with Mach number of incoming flow

4 结论

(1)发展的无网格计算方法能有效地进行化学非平衡流的数值模拟，能准确反映流场的流动特性，具有布点灵活、求解效率高、精确模拟化学非平衡流中如激波等复杂流动现象的特点。采用无网格方法模拟火箭燃气射流流场，为数值研究燃气射流流场开辟了一种新途径。

(2)在欠膨胀状态下，相交射流激波在轴线上的交点随着欠膨胀程度的增加而远离喷管出口，且相交射流激波呈曲线;在过膨胀状态下，相交射流激波在轴线上的交点，随着过膨胀程度的增加而不断靠近喷管出口，且相交射流激波呈直线。

(3)在欠膨胀或过膨胀状态下，射流流场轴线上离开喷管一定距离后，不同喷管出口燃气射流压力对应下轴线压力也基本一致，且均约为外流的静压。

(4)超声速来流条件下，仅改变不同来流马赫数，对射流流场总体结构基本无影响，相交射流激波相交点位置基本相同，轴线上参数分布也基本一致。

图13 不同来流马赫数下轴线主要参数分布Fig．13 Distribution of main parameter along the axis variation with Mach number of incoming flow

［1］ 张福祥．火箭燃气射流动力学［M］．哈尔滨:哈尔滨工业大学出版社，2004．

［2］ 孙振华，徐东来，何国强．飞行参数对导弹发动机羽流的影响［J］．固体火箭技术，2005，28(3):188－191．

［3］ 聂赟，乐贵高，马大为．矢通量分裂法在喷管超声速喷流计算中的应用［J］．固体火箭技术，2013，36(1):56－60．

［4］ 常见虎，周长省，李军．高低空环境下火箭发动机射流场的数值分析［J］．系统仿真学报，2007，19(16):3672－3675．

［5］ Balasubramanyam S，Raghurama Rao S V．A grid free upwind relaxation scheme for inviscid compressible flows［J］．International Journal for Numerical Methods in Fluids，2006，51(2):159－196．

［6］ Luo H，Baum J D，Löhner R．A hybrid Cartesian grid and gridless method for compressible flows［J］．Journal of Computational Physics，2006，214(2):618－632．

［7］ Jaisankar S，Shivashankar K，Rao S V R．A Grid－free central scheme for inviscid compressible flows［R］．AIAA 2007－3946．

［8］ 陈红全．隐式无网格算法及其应用研究［J］．空气动力学报，2002，20(2)．

［9］ 孙迎丹，王刚，叶正寅．AUSM+－UP格式在无网格算法中的推广［J］．空气动力学报．2005，23(4)．

［10］ 胡世祥，李磊，佘春东，等．二维Euler方程无网格算法的精度分析［J］．计算力学学报，2005，22(2):232－236．

［11］ 周星，许厚谦．计算含运动边界非定常流动的无网格算法［J］．力学与实践．2010，32(3):16－25．

［12］ 马新建，谭俊杰，任登凤．三维无网格及其在超声速弹丸流场模拟中的应用［J］．弹道学报．2010，22(3):54－57．

［13］ 周星，许厚谦．无网格算法在膛口流场模拟中的应用［J］．弹道学报．2011，23(1):24－26．

［14］ Kim K H．Accurate computation of hypersonic flows using AUSMPW+scheme and shock－induced grid technique［R］．AIAA 98－2442．

［15］ Charles J．Nietubica．Navier－stokes computations for a reacting，M864 base bleed projectile［R］．AIAA 93－0504．

［16］ 刘君，周松柏，徐春光．超声速流动中燃烧现象的数值模拟方法及应用［M］．长沙:国防科大出版社，2008．

［17］ Viguier C，Silva L F F，Desbordes D，et al．Onset of oblique detonation wave:comparison between experimental and numerical results for hydrogen－air mixtures［C］//Symposium(International)on Combustion．Elsevier，1996，26(2):3023－3031．