高马赫数下超声速燃烧的自点火查表方法1)

张锦成 王振国 孙明波 汪洪波 王亚男 刘朝阳

(国防科技大学空天科学学院,长沙 410072)

引言

燃烧室是高超声速冲压发动机的核心部件,其内部超声速气流中燃料高效稳定的燃烧一直是备受关注的课题[1].高马赫数下,燃料的燃烧驻留时间极短难以充分预混,而燃料与氧化剂的掺混对点火及火焰稳定的影响非常重要[2-3].化学反应和流动的强耦合性是超声速反应流场的另一个关键问题.许多实验和数值工作都以建立更准确精细的湍流化学反应模型为目标.本文借鉴火焰面的降维建库查表思想,初步发展了具备描述高马赫数自点火效应的查表方法.进一步与大涡模拟求解器耦合形成高马赫数反应流求解器.有利于进一步理解自点火在火焰传播和稳定的作用并为高马赫数工况冲压发动机设计提供高效的仿真平台.

高马赫数下超声速燃烧的自点火效应非常普遍,并且对燃烧室内火焰稳定有重要作用[4-6].Cheng等[7]开展了一系列同轴射流实验,观察到了稳定在下游的抬举火焰.Boivin 等[8]应用氢氧三步简化反应机理进行仿真,对比发现必须通过针对自点火效应的反应速率修正,可以得到与详细机理较吻合的结果;捕捉流场自点火火焰基和描述湍流扩散火焰,并提出基于化学爆炸模型分析(CEMA)的自点火识别方法,结果表明上游“指状”自点火核是流场最初反应阶段.Moule 等[9]在文献[7]实验的基础上进行数值仿真,对火焰稳定机制进行分析.指出自点火效应和可压缩效应共同控制化学反应过程,自点火过程通过生成水来开启链式反应,下游激波后的高温反应区主导释热,使抬举火焰稳定.Bouheraoua 等[10]进行相同工况的数值研究,发现流场中激波和点火的相互作用:入口激波提供了自点火发生的速度和混合条件,反应释热会在同一位置诱导出间歇的激波,其导致的压力波动驱动自点火火核位置和强度产生不稳定.Zhao 等[11]通过采用理想搅拌器模型考虑湍流影响结合化学反应建库方法捕捉到了剪切涡中自点火火核,同时能够对实验的温度场提供吻合很好的预测.

Ben-Yakar 等[12]开展了飞行Ma=10 的横向射流实验,高速纹影和OH-PLIF 结果显示OH 基主要分布在射流迎风剪切层和近壁区域.Won 等[13]对Ben-Yakar 实验构型和工况进行数值仿真,所用的DES (detached-eddy simulation)耦合有限速率模型,成功重现了高焓来流中横向射流形成的三维非定常反应流场.计算结果说明源于上游回流区反向涡对(CVPs)相互作用持续地产生分离涡.Liu 等[14]也进行了相同条件的数值工作,采用高精度格式捕捉到了大尺度涡结构的产生和输运;进行了化学反应特性分析,发现生成OH 基的链反应是吸热的,并且发生在贫燃一侧而释热明显的链反应消耗OH 基同时生成水发生在富燃一侧,总结了上游自点火维持下游扩散火焰的机制.Oliver 等[15]和Goldfeld[16]最新的工作揭示了自点火在横向射流火焰时空演化中的作用.低释热的自点火发生在贫燃一侧,主要分布在上游回流区,释热的扩散导致自点火开始级联.发生在富燃侧高释热的自点火,导致的膨胀扩展产生褶皱的火焰锋面.

数值研究中,体现化学反应和湍流的交互,再现反应流场多时空尺度的演化非常困难.现有模型往往建立在一些假设上,只适用于一些条件下的仿真,许多针对性的模型改进由此开展.实际燃烧室的仿真模拟计算量大,所用燃料反应机理复杂,需要更加高效且适用于详细化学反应机理的模型方法.化学反应的查表方法可以降低组分输运带来的计算量,同时可以使用详细机理.

现有的处理湍流化学反应的方法主要有两类:非降维方法和降维方法.典型的非降维方法为设定型PDF (assumed PDF)方法[17-18].该方法通过已知的高阶矩信息设定概率密度函数,可以给出反应速率在包含脉动量的状态下的值.通过设定温度高斯分布,组分联合β分布[19-20]可以得到与实验结果较为一致的仿真结果.但是当使用详细化学反应机理时,设定型PDF 方法需要求解组分输运方程,导致巨大的计算负荷.内在低维流形(ILDM)是最早也是最简单的化学反应降维建库方法,由Maas 和Pope[21]提出.由于ILDM 只采用了反应平衡后的状态进行建库,在反应预热的低温区域没有被覆盖,只能通过线性外推预测,会导致较大的偏差.扩展的火焰低维流形(FPI)方法[22]作为ILDM 扩展和改进被提出,能够对低温区反应得到更好的预测结果.FPI 使用了不同混合分数的一维层流预混火焰来描述低温区.自点火发生时反应速率非常低,特征时间大于局部流动特征时间,以上方法不能被准确描述非平衡过程.针对自点火现象,点火动力学查表(TKI)是一个被成功应用于低速压燃发动机中的方法[23].TKI 采用一定初始温度和压力下的等压预混反应器,在不同进度标量,生成燃料消耗速率和释热率数据库[24].但是其在高超声速反应流还没有进一步的应用.经典的稳态扩散燃烧火焰面模型引入薄反应层假设(Da≫ 1),建立火焰面方程.垂直于火焰面可以用混合分数标识,火焰面上可以用标量耗散率表示不同的状态,从而进行降维.稳态扩散火焰面模型[25-26]不包含非稳态燃烧分支,不能描述反应发生的瞬态过程.因此,火焰面进度变量模型被提出,增加一个定义的进度变量来标识反应进行的程度[27-28].但是自点火发生时,火焰面薄反应层假设不再成立,限制了火焰面方法在自点火效应明显的高马赫数反应流中的应用.

本文在火焰面进度变量模型的思路基础上,提出在时间维度上降维以考虑自点火效应,建立一种超声速燃烧的自点火查表方法.能够与大涡模拟求解器耦合,适用于高马赫数反应流的仿真应用.介绍了自点火数据库的生成和映射过程,阐述流动求解器的耦合方式和参数更新方式,然后针对无流动的零维化学反应进行验证.最后将耦合自点火查表方法的大涡模拟应用在两种典型的高焓超声速模型冲压发动机中,证明该方法的有效性和适用性.

1 模型与方法

1.1 可压缩大涡模拟方法

本文的模型实现和流场仿真工作以大涡模拟(LES)求解器为平台.该求解器基于有限差分法求解超声速反应流场.为保证数值方法的相容性和稳定性,方程中不同数学物理性质的项采用不同的数值格式进行处理.对于空间离散项,采用五阶WENO格式重构对流项,而二阶中心格式离散黏性项.对于时间离散项,使用隐式LU-SGS 算法进行时间积分.有关方程和数值方法的详细说明,请参见文献[4,29].

1.2 自点火查表方法

借鉴于火焰面方法对火焰空间结构的降维,可以定义时间序列上一组热力学状态向量Vn.通过解如下的化学反应动力学常微分方程求得随时间变化的热力学参数Vn(如图1 所示)

其中 ωi由基元反应模型和化学反应动力学公式[30]计算.图1 展示了初始温度1050 K,当量比为1,压力100 kPA 的H2-O2反应曲线.通过在曲线上采集不同时刻的点,记录诸如组分质量分数,温度,反应速率,释热率等信息组成一组原始数据.计算不同初始状态下的化学反应曲线,并进行数据采集则可以形成原始数据库.可以通过初始状态和不同初始状态下的反应时间索引到每个时间序列上的所有数据点.考虑影响化学反应的主要因素,初始条件选定为混合分数Z,初始温度T0,压力p.初始条件设置的范围和取值如表1,覆盖了冲压发动机燃烧室内涉及的大部分参数范围,同时为了控制数据量对混合分数设置了不均匀的取值.

图1 化学反应动力学计算出的温度和H2 曲线Fig.1 Temperature and H2 curves calculated by chemical reaction kinetics

表1 初始参数的范围和取值Table 1 Range and value of control parameters

值得注意的是,化学反应计算需要使用的参数随时间变化是高度非线性的.图2 展示了原始数据库中产物H2O 反应速率和反应释热率随时间变化的曲线,释热率和产物反应源项是化学反应与流动耦合的重要参数;反应延迟时间,对流场中火焰分布位置和结构有重要影响.建立自点火模型的重点是:准确预测延迟时间,反应速率和释热率的极值;尽量还原反应速率和释热率随反应进行的分布细节.

图2 H2O 的生成率和反应释热率Fig.2 Generation rate of H2O and reaction heat release rate

在流场仿真计算中,反应经历时间和全局的物理时间是不一致的,通过其他方法也不容易求得,这给自点火查表耦合流动求解器带来困难.通过组分构造一个随时间单调变化的新变量,一方面可以代替时间进行数据库索引,另一方面代替所有组分在流场中输运,这是数据库方法实现反应流动计算的关键.基于以上考虑在生成数据库时要考虑将时间映射到一个新的变量,能够通过组分反应反应进度.类似火焰面/进度变量模型,定义自点火进度变量

由于YC,max只取决于初始条件,这样的归一化是可行的.由于水的积累是单调的,所以归一化过后的变量仍然保持对时间的单调性.进一步需要对进度变量的采样点在点火延迟之前加密,以保证对自点火行为的精确捕捉.加密方式如图3 所示,局部加密还可以减少数据库的规模.

图3 进度变量采样点分布Fig.3 Distribution of value points for progress variables

图4 展示了OH 基质量分数和温度随进度变量的变化,与图1 对比可以看出进行映射的效果:反应过程在进度变量维度上变化的更均匀,这样降低了化学反应在时间维度上的刚性.图5 展示H2反应速率和进度变量反应速率曲线,进一步证明了选取进度变量的优势,反应速率曲线更光滑,能保证插值计算的精度.

图4 OH 基质量分数和温度随进度变量的变化Fig.4 Variation of OH mass fraction and temperature with progress variables

图5 H2 反应速率和进度变量生成率变化曲线Fig.5 H2 generation rate and progress variable generation rate curve

1.3 耦合查表与验证

本节首先介绍数据库的查表插值方法以及与大涡模拟求解器耦合方法.之后以零维的预混反应工况,验证数据库生成方式和查表插值方式.

生成的数据库是一个四维 Φ (T0,p,Z,YC) 正交的大型数组,采用分段线性插值的方式进行查表输出.如图6 所示输入当前的 (T0,p,Z,YC),通过查找并进行线性插值,得到进度变量速率,释热率hrr,组分Yi,然后进行温度组分更新得到下一时间步的参数.

图6 数据库查表插值示意图Fig.6 Schematic diagram of lookup and interpolation from the auto-ignition database

在与大涡模拟耦合求解流场时,首先要在流场中输运混合分数和进度变量用来代替组分的输运,并模拟自点火火核的传播和火焰扩散.输运方程与火焰面/进度变量方法中构造的方程一致

其次定义数据库查询的输入变量,根据求解方程组1 的定压和守恒条件,从大涡模拟求解器输入数据库的 (T0,p,Z,YC) 可以如下定义:温度根据当地的混合分数定义纯混合温度T0=Tox(1−Z)+TfZ,其中Tox,Tf分别为氧化剂测温度和燃料一侧温度,Z为混合分数.p,Z为当地的压力和混合分数,YC由输运方程求得.

然后通过查表插值后获得的进度变量速率,释热率,组分,对流场进行更新.释热率,组分用于更新温度,采用牛顿下山法进行迭代求解.进度变量速率用于更新进度变量,并添加到进度变量方程的源项上.

化学反应过程中遵循质量和原子守恒定律,如下定义的混合分数Z不随时间变化

图7 混合分数Z 守恒的验证Fig.7 Verification of mixture fraction Z conservation

对初始温度T=1200 K,压力p=100 kPa,混合分数Z=0.0283 (当量比为1)工况进行计算,对比自点火数据库插值方法和化学反应速率积分方法得到的结构来验证组分求解的准确性.图8 为插值和积分方式得到的H2和H2O 质量分数在反应中的变化,图9 为不同初始温度下两种方式温度随时间和进度变量的变化.图中红色带方块标记的是查表进行插值得到的曲线,蓝色圆圈标记的是反应速率积分的结果.积分结果可以当作验证数据库插值结果的基准.其中,作为进度变量的H2O 的质量分数吻合情况最好,说明进度变量的选取和进度变量速率的插值计算都是合适的.H2在快速反应阶段吻合较好,在反应稳定阶段略偏低.总体上看自点火数据库方法基本上可以实现化学反应的计算,可以成为有限速率方法的高效率替代方法.

图8 组分H2 和H2O 质量分数在反应中的变化Fig.8 Variation of H2 and H2O mass fraction in the reaction

图9 不同初始温度下温度随时间变化Fig.9 Temperature variation with time at different initial temperatures

本小节主要验证自点火数据库方法的三个方面:混合分数守恒性,主要组分和温度更新正确性,为该方法在超声速反应流仿真中的验证提供基础.

2 算例验证

当冲压发动机燃烧室工作在高马赫数工况下,入口来流空气的总焓非常高.很多研究表明此时氢气燃料能够立即与热空气发生反应,无需借助任何稳燃装置能够形成稳定的火焰结构.在简单的燃烧室构型产生的复杂射流燃烧流场,包含自点火与火焰之间的复杂的演化与竞争.本节基于 Gamba 的横向射流燃烧实验和Burrows-Kurkov 台阶射流实验,针对高焓环境下的氢气自点火主导的燃烧过程开展基于LES 的自点火查表方法的应用验证.

2.1 高焓超声速来流中横向射流的自点火

Gamba 和Mungal[31]通过 OH 自发辐射和纹影技术研究了不同动量比(J=0.3～ 5.0)的声速射流喷注进入超声速来流诱导的混合燃烧特性.发现动量比J直接控制着下游的反应模式和分布区域.实验条件下自由来流马赫数 2.4,其滞止温度可达到3000 K,相应的静压和静温分别为 40 kPa 和 1400 K.一个直径 2 mm 的圆形喷注器安装在距离平板前缘64 mm 位置.氢气射流以声速垂直喷注到超声速流场中,其静压为 1074.5 kPa、静温为 250 K.计算域构型和参数如图10,展向宽度为80 mm.

图10 Gamba 实验构型计算域示意图 (单位:mm)Fig.10 Schematic diagram of the computational domain of the Gamba’s experiments (unit:mm)

燃料射流沿垂直壁面方向喷注到高焓超声速气流中,会形成非常复杂的反应流场.射流与来流之间的强相互干扰诱导出一系列激波和湍流涡等结构,在此过程中氢气与热空气快速掺混并发生剧烈反应.

图11 展示了实验OH-PLIF 和计算得到的中心截面上OH 分布云图.实验和仿真均表明化学反应主要集中在两个区域:一个靠近壁面,此处静温高点火延迟低;并且此处流动速度较低,燃料驻留时间长混合充分;因此可以维持燃烧稳定.另一个反应区位于射流迎风剪切层,形成OH 基浓度很高的反应锋面.尽管当地速度很高,但是在湍流涡的作用下自点火火核未被吹灭,并能够输运到流场下游.说明自点火查表方法能够准确捕捉到流场中的自点火现象,并能够反映自点火火核的传播.

图11 中心截面上OH 分布云图Fig.11 OH distribution on the central section

2.2 台阶壁面燃烧室中射流诱导的抬举火焰

Burrows-Kurkov 的实验[32]采用了台阶壁面的燃烧室,燃料氢气在台阶侧壁以来流同向的速度注入由一定比例的H2O,N2,O2混合的污染空气来流.台阶底壁略有扩展,具体的几何构型见图12.马赫数2.4 的来流以1237.9 K 静温,96 kPa 静压在台阶上游矩形入口进入,根据一些基于此实验的仿真研究入口处含有1 cm 厚度的边界层;声速的261.7 K的低温射流以略高于主流的压力喷注.来流和射流的参数见表2.

表2 Burrows-Kurkov 实验射流和来流参数Table 2 Jet and inflow parameters in Burrows-Kurkov experiment

图12 Burrows-Kurkov 实验构型计算域示意图 (单位:mm)Fig.12 Schematic diagram of the computational domain of the Burrows-Kurkov experimental configuration (unit:mm)

燃烧室内燃料和氧化剂剪切形成的涡的存在提供了混合的场所,自点火一般发生在剪切涡贫燃一侧,局部温度升高会降低点火延迟时间;自点火会发生级联形成火焰基,向下游传播形成稳定的湍流火焰.由于存在点火延迟时间,火焰被抬举维持在距射流出口一定位置.

图13 为自点火数据库查表法预测的中心截面上温度和OH 基的分布,高温区在上游成破碎的指状,而在下游高温区连片状.这里由于上游的自点火火核的扩散和输运形成了稳定的火焰.氢氧基作为很重要的组分用来标记化学反应的锋面,氢氧基出现的位置较反应升温明显的火焰区域提前,说明上游存在低温自点火反应.自点火查表方法再现了自点火发展形成火焰的过程,证明其适用于火焰自点火共存的复杂燃烧流场,扩大了建库类方法的应用范围.但是细节上,本方法预测的火焰抬举位置略低,可能是由于还没有针对湍流对火焰的作用进行修正.

图13 中心截面上温度和OH 基的分布Fig.13 Distributions of temperature and mass fration of OH at the central slice

图14 通过出口位置(x=356 mm)总温剖面定量对比了实验数据和仿真结果.可以看到在近壁的反应区,总温的上升与实验数据吻合很好,远离反应区的远壁端总温也基本与实验一致.但是得到的反应区中心位置的总温最大值略比实验所得高.这是由于数据库方法还是层流状态,使得预测的反映温度偏高.定量的对比证实了自点火数据库查表方法的可靠性.

图14 出口位置(x=356 mm)总温剖面Fig.14 Total temperature profile at the outlet (x=356 mm)

3 结论

超声速高焓来流条件下,燃料室可能发生自点火行为对火焰稳定有重要意义,并且给超声速燃烧问题建模和仿真研究带来了新的挑战.本文提出一种超声速燃烧的自点火查表方法.首先介绍基于化学反应动力学建立自点火数据库的思路并就零维预混火焰验证了数据库的建立和插值方法.然后对两个三维射流诱导的燃烧流场进行仿真,对基于大涡模拟的自点火数据库查表方法进行验证,可以得到如下结论.

(1)选择H2O 作为进度变量,降低了求解化学反应时的刚性.并且通过进度变量成功地耦合了流动求解和反应求解,可以正确地处理流动与化学反应的相互作用.

(2)算例验证可以证明本文提出的方法初步具备对超声速条件下包含自点火模式的反应流场进行仿真的能力.