摘要
为研究火箭发动机横向不稳定燃烧特性,采用详细化学反应机理(GRI Mech 3.0)建表的小火焰生成流型,对模型火箭发动机中出现的横向不稳定燃烧进行数值模拟。通过与实验数据对比验证了模型的准确性;采用动态模态分解对压力场进行分析,研究了流场的动态特性;结合瑞利因子定量分析了不稳定燃烧的驱动特性。结果表明,数值模型能够有效捕捉横向不稳定燃烧,其主频与实验值相差不到1%;燃烧室横向压力振荡与喷嘴氧管纵向压力振荡相耦合,引起推进剂质量流量振荡;不稳定燃烧的驱动源主要位于燃烧室两侧,最边缘喷嘴对维持不稳定燃烧的贡献最大;推进剂与燃烧室侧壁面的相互作用极大增强了释热脉动,周期性释热为压力振荡提供能量,形成了不稳定燃烧极限环。
Abstract
To study the transverse combustion instability characteristics of the rocket combustor, numerical simulations of transverse combustion instability in a model rocket combustor were conducted based on the detailed chemical reaction mechanism (GRI Mech 3.0) and the flamelet-generated manifolds method. Accuracy of the numerical model was verified by comparing it with the experimental data. Pressure field was analyzed by the dynamic mode decomposition method, and the dynamic characteristics of the flow fields were investigated. Driving characteristics of combustion instability were quantitatively estimated by Rayleigh index. Result shows that the transverse combustion instability that occurred in the experiment can be effectively captured by the numerical model. Dominant frequency identified by the numerical study differ from the experimental value by less than 1%. Transverse pressure oscillations in the combustion chamber are coupled with that the longitudinal mode in the oxidizer post, leading to the pulsated propellant mass flow rate. Driving regions of combustion instability are mainly located on both sides of the combustion chamber, and the most marginal injectors played a critical role in keeping combustion instability. Heat release pulsations which periodically provide the energy source for the pressure oscillations are highly enhanced by the interactions between the propellant and the sidewall of the combustion chamber. And the combustion instability limit-cycle is formed.
高频不稳定燃烧是在一个受限空间内热声耦合的振荡过程,它的出现往往伴随着尖锐的啸叫以及压力和释热的脉动。液体火箭发动机不稳定燃烧是一个从20世纪40年代持续至今的问题,其产生机理至今仍未认识清楚[1],已经成为限制大推力火箭发动机发展的重要因素之一。从历史上看,几乎每一种火箭发动机在研发过程中都会遇到不同程度的燃烧不稳定困扰[2]。由于火箭发动机具有极高的功率密度,燃烧室的声学阻尼较低,只要将总释热量的很小一部分转移到声场中即可产生快速增长的压力振荡,诱发产生不稳定燃烧[3],从而使燃烧室壁面的热负荷急剧增加,液膜冷却失效,进而损坏发动机,产生严重后果。
瑞利准则[4]是理解不稳定燃烧的基础,由瑞利准则可知不稳定燃烧的产生与压力振荡和释热振荡的耦合密切相关。在后来的研究中,人们采用瑞利因子来衡量不稳定燃烧被驱动或抑制。多年来,学者们发展了多种瑞利因子[5-7]。本文采用的瑞利因子是由Harvazinski等[5]提出的,定义为:
其中,p和q为瞬态压力和释热率,“”表示取平均值,t为时间。
瑞利因子是压力振荡和释热振荡在单个体网格中的积分。积分值为正,则表示不稳定燃烧被驱动,反之则代表不稳定燃烧被抑制;积分值越大说明对不稳定燃烧的驱动作用越显著。
对全尺寸发动机的热式车研究是洞察不稳定燃烧最直接的方法,但采用全尺寸发动机研究不稳定燃烧不仅周期长、费用高,而且获得的数据有限。全尺寸发动机燃烧室中的高温高压环境导致很难对流场进行诊断,对揭示不稳定燃烧机理的贡献较小[8]。鉴于此,缩比模型发动机成为世界各国研究燃烧不稳定性的有效手段。美国普渡大学采用单喷嘴模型发动机[9-10],结合数值仿真[11],详细研究了纵向不稳定燃烧;采用二维矩形模型发动机系统研究了火箭发动机的横向不稳定燃烧[12-15]。德国宇航中心采用多喷嘴模型发动机研究了低温推进剂的切向不稳定燃烧[16-17],发现燃烧室切向压力振荡与喷注器纵向压力振荡相耦合。在国内,Bai等[18] 和Cao等[19]采用单喷嘴模型发动机研究了喷雾自激振荡与不稳定燃烧的关系,发现一定条件下自激振荡能够引起不稳定燃烧。Guo等[8,20-22]研究了模型发动机的纵向与横向不稳定燃烧,发现喷注耦合是诱发不稳定燃烧的重要因素之一。近年来,随着计算能力和数值模型的发展,计算流体力学(computational fluid dynamics,CFD)数值计算已经成为研究不稳定燃烧的重要手段之一,尤其是模型实验与CFD结合的方法得到了充分的发展[22-23]。CFD数值计算能够准确捕捉到火箭发动机的不稳定燃烧,获得丰富的流场信息,得到详细的流场动态演变,对研究不稳定燃烧机理意义重大。深入研究不稳定燃烧机理,深化对不稳定燃烧动态特性的认识,对抑制不稳定燃烧具有现实工程意义。
采用详细化学反应机理(GRI Mech 3.0[24],53组分325步),结合火焰生成流型(flamelet-generated manifolds,FGM)方法对模型火箭发动机开展了数值研究。成功复现了实验中出现的不稳定燃烧问题,分析了模型发动机横向不稳定燃烧的动态特性,初步揭示了不稳定燃烧维持机理。
1 数值计算模型
1.1 横向不稳定燃烧模型发动机
图1所示为矩形模型发动机示意图。它主要由氧化剂腔、喷嘴、燃烧室和喷管组成。它是Orth等[12]开发的一种矩形截面,用于研究横向不稳定燃烧的模型发动机。通过精心设计,令矩形燃烧室的横向模态与全尺寸燃烧室的横向模态相匹配,这样就能使模型发动机一定程度上模拟全尺寸火箭发动机中的横向不稳定燃烧[25]。如图1所示,9个同轴剪切喷嘴呈横向均匀分布在矩形燃烧室头部。在氧化剂腔的上游(图中未显示)还包含预燃室,通过氢气在富氧环境中燃烧产生高温(636 K)氧化剂。氧化剂包含质量分数为3.5%的水蒸气和约96.5%的氧气。氧化剂通过包含声学截止入口的中心氧管供应进入燃烧室。气态甲烷通过同轴环缝,经燃料腔进入燃烧室。喷嘴包含缩进段,在缩进室内燃料与氧化剂相互作用,然后进入燃烧室完成燃烧过程。模型发动机详细的结构介绍以及几何尺寸参见文献[12],在此不再赘述。本文的仿真是在Orth等[12]实验工作的基础上进行的,其目的在于揭示模型发动机中横向不稳定燃烧维持机理。因此,本文数值仿真了实验中出现横向不稳定燃烧最剧烈的工况之一(=213.21 g/s,=757.53 g/s)。

1.2 湍流燃烧模型
本文采用非稳态雷诺时均(Reynolds-averaged Navier-Stokes equations,RANS)方法,对于湍流的处理采用了由Menter[27]提出的k-ω 剪应力输运(shear stress transport,SST)模型,它结合了k-ω和k-ε方法的优点。此模型在近壁面处采用k-ω模型,能够准确解析壁面附近的剪应力传输过程;在自由剪切层则采用k-ε方法。k-ω SST模型由k和ω两个方程组成,详细的理论可见文献[27]。
燃烧模型采用了基于部分预混燃烧的FGM模型[28]。在本研究中,燃料和氧化剂在进入燃烧室之前已经在缩进室进行了剧烈的相互作用,可以认为推进剂在进入燃烧室之前已经处于部分预混状态。因此,本文采用了基于部分预混燃烧的FGM模型。大量的研究[8,20,22]表明,此方法能够很好地描述湍流与火焰的相互作用,以及捕捉火箭发动机中的不稳定燃烧过程。
本文采用详细的化学反应机理(GRI-Mech 3.0)生成预设概率密度函数(probability density function,PDF)表,它将物理空间的层流扩散火焰方程映射到混合分数空间,通过求解混合分数方程,结合查表方法能够大大降低计算量,特别适合采用详细的化学反应机理。组分的质量分数方程为:
(1)
温度方程为:
(2)
其中,Yi为第i种组分的质量分数,T、ρ 和 f 分别是温度、密度和混合分数,cp是混合平均的比热容,cp,i、Hi和 Si 分别为第i种组分的比热容、比焓和反应速率。其中标量耗散率χ定义为:
(3)
其中,erfc(·)为互补误差函数。
假设混合分数f和反应进程变量c是相互独立的,则f和c联合PDF可以表示为:
(4)
其中,取决于平均混合分数(Favre平均)及其二阶矩′;同理,取决于平均反应进程变量及其二阶矩′。
基于上述模型,通过求解、′、和′的输运方程,结合查表就能确定物理场的标量值,在考虑详细化学反应的基础上大大缩减了计算量。FGM模型的详细理论可参见文献[29]。
1.3 计算域及边界条件
计算域及监测点分布如图2所示。如图2(a)所示,计算域包含了声学截止下游的氧管、燃料喷嘴、燃烧室和喷管。氧化剂和燃料入口为恒定质量流量入口。与实验相同,氧化剂温度设置为636 K,包含质量分数为3.5%的水蒸气。出口被设置为压力出口,初始压力为101.325 kPa,壁面为绝热无滑移壁面。网格采用六面体结构,数量约为353万个。图2(b)所示为计算域的监测点分布,其中为了便于与实验直接对比,probe_R1的位置与实验中高频压力监测点位置一致。probe_L1、probe_C1、probe_R1在一条直线上,其中probe_C1在燃烧室的中间,根据实验中出现的不稳定燃烧模态,probe_C1大约在压力波节的位置,probe_L1和probe_R1大约在压力波腹位置。这样,结合上述三点的压力相位关系与特征频率,可以初步判断燃烧室中出现的不稳定燃烧模态。监测点probe_ox位于最右侧喷嘴中心,位于缩进室内部,距离喷嘴出口5 mm。

图2计算域及监测点分布
Fig.2Computational domain and distributions of probes
1.4 离散方法
本文的数值计算基于三维FGM模型,采用商业CFD求解器ANSYS Fluent V19.2。本文中湍流模型采用k-ω SST。标量通量采用二阶迎风格式离散,采用压力-隐式分裂算子(pressure-implicit with splitting of operators,PISO)算法处理速度压力耦合。动量方程的离散方式为二阶迎风格式。为取得更准确的结果,采用双时间步进行求解。
1.5 网格无关性验证
采用六面体结构网格对计算域进行空间离散,网格如图2(a)所示。为了保证较高的时间分辨率,仿真采用的时间步长为5×10-7 s。在进行数值研究之前,首先进行了网格无关性验证,以确保仿真结果的合理性。采用三种不同数量的网格Mesh1(265万)、Mesh2(353万)和Mesh3(465万)进行了网格无关性验证。在数值计算过程中,三种网格均能捕捉到自发的横向不稳定燃烧。压力值的均方根(root mean square,RMS)-时间轨迹如图3所示,可以看出当网格量达到353万后压力均方根对网格量不再敏感。网格无关性验证总结如表1所示。可以看出,Mesh1在平均压力、压力均方根和振荡主导频率上和Mesh2、Mesh3均存在明显的差异。对比Mesh2和Mesh3发现压力均方根和振荡主导频率相差很小,说明其对网格已不再敏感。因此,本文后续的结果均是基于Mesh2。

图3三种网格下的压力均方根-时间轨迹
Fig.3Pressure RMS-time traces for three grids
表1网格无关性验证结果
Tab.1 Results of grid independence study

2 结果与讨论
2.1 不稳定燃烧的频谱特性
图4所示为监测点probe_L1、probe_C1和probe_R1的压力时间轨迹。从图中可以看出,在燃烧室两侧(probe_L1、probe_R1)压力振荡幅值较高,但相位几乎完全相反,说明燃烧室的两侧为压力振荡的波腹。在燃烧室中心处(probe_C1)压力振荡幅值明显降低,频率增加(几乎是两侧振荡频率的两倍),说明燃烧室中心为压力振荡的波节。结合图4,在数值计算过程中得到的平均室压为1.17 MPa,而实验中的平均室压为1.14 MPa[12]。数值计算得到的平均室压与实验值相差仅为2.63%,说明数值计算模型具有适当的合理性。数值计算得出的平均室压比实验值略高,这主要是因为实验中燃烧不充分及其过程中的热量损失。
主导频率是不稳定燃烧的重要动态特性之一,因此功率谱密度(power spectral density,PSD)分析能够提取变量振荡主导频率,通过对比实验与数值仿真得到的主导频率能够进一步验证数值模型的合理性。图5(a)为实验过程中燃烧室右侧压力监测点的PSD图,可以看出燃烧室的压力振荡主频为2 625 Hz,对应为一阶横向不稳定燃烧。由于不稳定燃烧的非线性特征,又伴随着一系列的高次倍频谐波。图5(b)为相同位置(probe_R1)处数值计算得到的压力振荡PSD图。对比图5(a)和图5(b)可以看出数值计算得到的主导频率与实验值相差不足1%,这比之前的数值仿真更加精确[26]。同时,数值仿真过程中也捕捉到了多个高次倍频谐波,进一步说明了数值计算模型的合理性,展现了考虑详细化学反应机理的模型优势。

图4probe_L1、probe_C1和 probe_R1的压力时间轨迹
Fig.4Pressure time traces of probe_L1, probe_C1, and probe_R1
图6所示为监测点probe_L1和probe_C1的压力时间轨迹和色谱图。结合上述分析可以看出,数值计算捕捉到了燃烧室中的一阶横向模式(first width mode,1W)不稳定燃烧,在压力振荡极限环形成后持续了整个计算过程,并伴随着高次倍频谐波二阶横向模式(second width mode,2W)和三阶横向模式(third width mode,3W)。在燃烧室中心处数值计算捕捉的振荡主频为5 394 Hz,对比实验值5 310 Hz [26],两者相差1.6%,进一步验证了数值计算的准确性。对比图6(a)和图6(b)可以看出,燃烧室中心处(probe_C1)的压力振荡主导频率是燃烧室边缘的两倍,呈二阶横向模态,其振荡幅值相对较小。这种特性与燃烧室中压力波的传播特性有关,燃烧室两侧为压力振荡波腹,而中心则为波节。压力波在燃烧室中横向传播,导致一个周期内压力波锋面经过燃烧室中心两次,因此燃烧室中心的主导模态为二阶横向模态。


图6probe_L1和 probe_C1的压力时间轨迹和色谱图
Fig.6Pressure time traces and spectrogram of probe_L1, and probe_C1
2.2 压力振荡的DMD分解
动态模态分解(dynamic mode decomposition,DMD)是一种纯数据矩阵驱动的模态分解技术,它不依赖于任何先验准则,旨在通过提取的低维数据结构重构原始数据模态。DMD能够将时间和空间模态分别提取出来,使得它特别适合分析具有特定频率的不稳定燃烧数据[30]。多年来,学者们发展了多种DMD方法,本文采用的是由Tu等[31]开发的exact DMD,其详细的理论可见文献[31],在此不再赘述。
本文采用中间截面的压力振荡数据进行了DMD分析,用振幅排序的方法按照各个模态对流场的贡献和影响进行排序。DMD提取的第一个模态(Mode_0)为静态模态,表示平均流场[32],其余模态均成对出现,其特征值为共轭复数,频率的绝对值相同,符号相反如图7(a)所示,因此每对共轭模态可以看作是单个模态,对应其正的频率值[33]。由图7(a)可以看出,DMD的前5阶模态基本落在单位圆上,结合图7(b)说明系统基本符合简谐运动特性,因此本文取DMD前五个模态进行分析。
图8展示了DMD分析的前5个模态,其中第一个模态(见图8(a))即静态流场命名为Mode_0,相当于平均压力场。图8(b)~(e)为4个共轭模态,分别为Mode_1~Mode_4,DMD分析得到的主导频率分别为2 663 Hz、5 321 Hz、7 999 Hz和10 626 Hz,与PSD分析得到的前4阶模态基本一致。图8(b)所示为DMD分解的Mode_1,其在燃烧室两侧展现了压力高振幅特性,在燃烧室中间压力振幅较低,这是典型的驻波型一阶横向不稳定燃烧,进一步确认了燃烧室中的不稳定燃烧模态。观察Mode_1~Mode_4可以看出前三阶模态在燃烧室中表现得最为明显,振荡幅值也比较高,说明Mode_1~Mode_3对燃烧室中的振荡特性贡献最大。如图8所示,对比燃烧室和氧管的模态可以发现,当燃烧室中出现一阶横向模态时氧管中伴随着一阶纵向模态,在出现二阶模态时氧管也会出现二阶纵向模态,说明在不稳定燃烧发生时氧管中也伴随着压力振荡,且氧管中压力振荡的纵向模态与燃烧室中的横向模态相耦合。

图7DMD 特征值和前五阶时间系数演化
Fig.7DMD eigenvalues and evolution of first five time coefficients

图8前5阶DMD模态(静态模态和4对共轭模态)
Fig.8The first 5 DMD modes with one static mode and four pairs of conjugate mode
2.3 不稳定燃烧的动态流场特性
如图9所示,给出了一个周期内静压、X方向速度和流线、O2质量分数和释热率的动态特性。如图9(a)t1所示,压力波向燃烧室左侧传播,同时,燃烧室中的气态混合物在压力梯度的驱动下向燃烧室左侧移动(见图9(b)t1)。在高压区域释热率较高,推进剂快速消耗,导致推进剂在燃烧室左侧集聚(见图9(c)t1)。压力波继续向左侧传播,直到与燃烧室壁面相互作用,导致局部压力迅速升高(见图9(a)t2)。如图9(c)t2所示,压力波携裹着推进剂气团撞击壁面,引起推进剂与燃烧室壁面的相互作用,增强了推进剂局部的混合效果。因此,燃烧室左侧出现了局部释热脉动并为压力振荡提供了能量,如图9(d)t2所示。燃烧室左侧的推进剂被迅速消耗,同时压力波经过壁面反射向燃烧室右侧传播(见图9(a)t3)。另外,燃烧室中的高压向上传播进入氧管,造成喷嘴短暂阻塞。推进剂在燃烧室中集聚,直到喷嘴中的压力高于燃烧室推进剂才能喷入燃烧室。因此,在燃烧室左侧的推进剂没能立刻得到补充,造成燃烧室中推进剂的不均匀分布,如图9(c)t3所示。在t4时刻压力波继续向右传播,同时释热峰也向右移动,为压力振荡提供能量。在t5时刻压力波撞击右侧壁面,经过壁面反射向左传播,准备下一个周期的循环。在维持压力振荡的过程中,压力波撞击燃烧室壁面,以及燃烧室横向压力振荡和喷嘴氧管纵向压力振荡的耦合起到了关键作用。
经过上述分析,燃烧室的横向压力波向上传播进入氧管,导致燃烧室横向压力振荡与氧管纵向压力振荡的耦合。因此,分析了喷嘴中监测点probe_ox的振荡特性,如图10所示。需要注意的是,为了实现在同一图中进行对比,各变量之间的相位关系经过了量纲处理。

图9一个周期内详细的流场演变
Fig.9Detailed evolutions of flow fields in a cycle
如图10所示,当压力波传入氧管后,氧管内压力迅速升高,速度和质量流量迅速下降。质量流率和Z方向速度甚至出现负值,说明氧管中出现了反流。氧化剂和燃料的入口条件为恒定质量流量入口,氧管压力升高,推进剂在喷嘴内积聚。当氧管内压力下降时,Z方向速度增加,推进剂加速进入燃烧室,引起推进剂质量流率周期性脉动,进而引起释热脉动。结果表明,氧管中的质量流量振荡对维持不稳定燃烧至关重要。

图10probe_ox的压力、Z方向速度和质量流量的标准化振荡特性
Fig.10Normalized oscillations of pressure, Z velocity, and mass flow rate at probe_ox
2.4 不稳定燃烧自持机理
瑞利因子是压力振荡与释热振荡的体积分,被广泛用于识别不稳定燃烧被驱动或抑制的区域。本文采用的瑞利因子计算方法与文献[5]中相同,积分区域为体网格单元,积分步长为10-5 s。图11(a)所示为喷注面板下游5 mm处燃烧室横截面的瑞利因子分布。由图可以看出,在燃烧室两侧瑞利因子为正值且较大,在燃烧室中心其值几乎为零。图11(b)为基于单个喷嘴宽度平均的瑞利因子,能够更明显地看出燃烧室两侧瑞利因子较大。仔细观察图11可以发现,位于最边上喷嘴处的瑞利因子最高,说明最边上的喷嘴对不稳定燃烧的驱动作用最大,这可能和最边上的推进剂与壁面的相互作用有关。

图11瑞利因子云图
Fig.11Contours of the Rayleigh index
综上所述,图12总结了本文中出现的横向不稳定燃烧机理。压力波在燃烧室中横向传播,同时高压向上传递进入氧管,导致氧管中的纵向压力振荡与燃烧室中的横向压力振荡相耦合。而氧管中的纵向压力振荡引起推进剂质量流量的脉动,在此过程中质量流量振荡的相位不断被燃烧室中的横向压力波调制,使得喷嘴在高压时储存推进剂,在低压时释放推进剂,从而导致释热脉动。燃烧室横向压力波携裹着未燃气团撞击壁面,使得局部压力上升,同时极大促进了推进剂混合。因此,燃烧室两侧出现巨大的脉动释热,为压力振荡提供了能量,形成不稳定燃烧极限环。

图12本文横向不稳定燃烧机理
Fig.12Mechanisms for the present combustion instability
3 结论
本文基于FGM方法对三维全尺寸模型发动机中横向不稳定燃烧进行了数值模拟。数值计算精确捕捉到了实验中出现的横向不稳定燃烧,得到的主频与实验值高度吻合,误差不足1%。通过数值研究初步得到了以下结论:
1)燃烧室中出现了一阶横向不稳定燃烧,且伴随着高次谐波,与实验现象高度吻合。与前人仿真相比[26],本文结合DMD分析发现燃烧室中的横向压力振荡与喷注器中的纵向压力振荡相耦合。燃烧室中的高压向上传递进入氧管,造成喷嘴质量流量周期性振荡,这种耦合关系对维持不稳定燃烧至关重要。
2)推进剂与燃烧室侧壁面的相互作用能够有效提高混合效率,产生巨大的释热脉动,导致燃烧室两侧的压力振幅比中间高。通过瑞利因子定量分析发现,不稳定燃烧的驱动源主要位于燃烧室的两侧,左右两侧最边上的喷嘴对驱动不稳定燃烧的贡献最大。
3)喷注器质量流量的周期性脉动,导致推进剂在燃烧室的不均匀分布,进而引起燃烧室局部释热脉动。释热脉动在推进剂与燃烧室侧壁面相互作用的过程中得到了增强,并为压力振荡提供能量,最终形成振荡极限环。