双旋流模型燃烧室高温区演化分析

2021-08-24高东硕

燃烧科学与技术 2021年4期

高东硕，金捷，沈硕，李敏，曾家，王方

(北京航空航天大学能源与动力工程学院仿真研究中心，北京 100191)

在航空发动机和燃气轮机燃烧器中，旋流火焰以良好的点火性能和稳定性得到广泛应用，且旋流火焰能在非常有限的小体积内进行高效率的能量转化[1].由于雷诺数较高，旋流一般处于湍流状态，流动发展过程中常出现多种不稳定现象如漩涡破碎和进动涡核(PVC)[2].进动涡核由剪切层的不稳定性产生，是旋流流场中的一种不稳定的大尺度拟序结构，PVC会影响附近流体使之出现进动特征.Widenhorn 等[3]对一种双旋流模型燃烧室无反应冷态流进行了数值研究与激光多普勒测速(LDA)结果作对比分析，验证了位于内剪切层的频率为1 514 Hz 的进动涡核的存在.而PVC 引起的速度脉动会导致燃烧场的不稳定性，Canepa 等[4]利用3-D 多普勒速度仪对一个双同向旋转的径向旋流器引起的不稳定流场进行了测量，表明回流区附近存在大尺度的涡脱落，涡脱落过程会影响燃烧的稳定性.

在湍流燃烧模拟中，精度高适用性强的湍流燃烧模型对于湍流燃烧的强烈非线性相互作用的解释更为合理.LES[5-6](large eddy simulation)作为现阶段描述细节最多而计算量可行的湍流模拟方法[7]，能够捕捉流场中的低频振荡运动，适合用于模拟有非定常特征的旋流流动.PDF[8](probability density function transport equation)是可以耦合详细化学机理的火焰模型，LES-TPDF 相结合能够得到比较详细的湍流与化学反应时间尺度，可以描述详细的瞬时流场和温度场，因此在有旋流场中得到广泛应用.

以往文献表明当量比的扰动与气动力学引起的涡旋形成和破碎过程会造成燃烧室的不稳定，对当量比的扰动的影响研究已经有很多，而对PVC 引起的旋涡与燃烧的不稳定关系的研究较少.本文预期通过对湍流火焰的数值模拟，分析燃烧场的数据，捕捉流场中的局部高温区，讨论进动涡核对燃烧不稳定性的影响，进而分析局部高温区演化机理.

1 计算程序与模型

1.1 程序介绍

AECSC(Aero Engine Combustor Simulation Code)是基于LES-TPDF 方法的三维亚声速两相湍流燃烧数值模拟程序，来源于帝国理工William Jones教授研究组，北京航空航天大学数值仿真中心发展了气相、两相和大规模并行版本.

AECSC 程序是FORTRAN 程序，采用模块化编程，拓展性好，计算精度高，原作者Jones 利用该程序检验并发展了许多模型，取得了良好成果.程序采用大涡模拟方法求解连续方程和动量方程，采用概率密度函数方法求解能量方程和组分方程.大涡模拟亚格子应力模型使用动力Smagorinsky-Lilly 模型，概率密度函数方程使用欧拉随机场解法.已有文献表明其能够模拟旋流流场，如Jones 等[9]在2013 年使用LES-TPDF 方法对旋流模型燃烧室进行两相模拟，与实验对比良好，火焰结构得到准确捕捉.

1.2 计算模型

本文选取的计算对象是一个双径向旋流模型燃烧室，由德国航空航天中心的Meier 等设计[10-11].实验[11]将测量区域简化为燃烧室主燃区，即在燃烧器的近场和受液体燃料主要影响的区域.为了方便实验测量，燃烧室由4 块石英玻璃构成，因此实验是在光学可达的高压单扇型燃烧室中进行的.利用LDA测量冷态流场速度，相位多普勒分析仪(PDA)用于测量雾化液滴的尺寸和速度，使用平面激光诱导荧光(PLIF)方法使火焰区域温度分布可视化.

GTMC 的几何结构如图1 所示，燃烧室截面为102 mm×102 mm 的正方形，长264 mm，实验使用的燃料为液态航空煤油Jet-A.空气通过中心旋流器和外旋流器进入燃烧室，两个旋流器的径向通道均为12 个，空气质量流量比接近1∶1.65，其设计细节如图2 所示.图3 展示了空气和航空煤油进口的设计细节，煤油由两个位置相对的燃油管道先进入一个环形的燃油通道，然后通过圆周均布的 36 个面积0.2 mm2大小的孔进入垂直槽，实验径向通道中高温空气会在垂直槽中给燃料带来显著的温升和预蒸发.燃料穿过垂直槽后以31°的倾角进入燃烧室，这会使得燃油的径向动量很小，避免了燃油从出口处脱离.

图1 GTMC的几何结构Fig.1 Geometric structure of GTMC

图2 GTMC中心旋流器和外旋流器轴向截面图Fig.2 Axial sections of GTMC center cyclone and outer cyclone

图3 空气和航空煤油进口的设计细节Fig.3 Design details of air and aviation kerosene inlets

1.3 网格信息

本文使用SolidWorks 软件创建几何模型，使用ICEM CFD 软件划分三维结构化网格，网格共划分为144 个区块，分别划分了网格数量为220 万、297 万与440 万的3 套网格，在Fluent 软件中使用3 套网格对模型燃烧室的冷态流场进行模拟，发现220 万网格的计算结果不是很好，而297 万与440 万网格的模拟结果基本一致，综合考量计算量与精确度，选择297 万的网格进行后续计算[12]，见图4.在选用的297 万网格中，轴向(流向)方向是Y 方向，图4(a)展示了X＝0 截面的网格，在燃烧室轴向方向上由于前段有着非常强烈的两股气流的掺混，这里的网格进行了加密，中后段的流动和反应都不是很强烈，网格调得稀疏一些.图4(b)展示了内外旋流器部分的网格细节，其中蓝色和绿色的方形部分为旋流器的空气进口.燃烧室的内外旋流器采用5 层O-block 进行加密.

图4 GTMC整体和旋流器部分的计算网格Fig.4 Computational grid of overall GTMC and swirler section

1.4 边界条件

实验的模型燃烧室工况见表1.本算例的燃料为航空煤油Jet-A，模拟时真实组分使用分子式C12H23代替，程序采用烷烃通用的四步七组分反应机理[13].

表1 实验工况Tab.1 Experimental conditions

2 模拟结果及高温区分析

2.1 燃烧场模拟结果

按照表1 中的冷态边界条件，使用AECSC 程序模拟冷态流场，实验测量了无反应流的速度分布，取距离喷射器出口轴向高度h＝10 mm 处的轴向、径向与切向速度与实验数据进行对比，如图5，3 个时均速度均与实验测量速度分布趋势相同，且曲线的极值点位置相近，验证了AECSC 程序对模型燃烧室冷态模拟的准确性.

图5 冷态工况下轴向、切向与径向时均速度沿径向分布Fig.5 Radial distribution of time-averaged velocity in axial，tangential and radial directions under cold condition

对于热态工况，先按表1 中边界条件模拟出没有喷油的冷态流场，以此为基础进行喷油并使用点火模型进行点火和燃烧.实验没有对进口燃油的温度，速度和d32进行测量，考虑到燃油在垂直槽中的温升，将初始液滴 d32设置为 20，液滴初始温度设为340 K，液滴初始速度为50 m/s.使用一个随机场求解PDF 模型.在点火成功，燃烧流场基本稳定后统计了3 个周期，约30 万步.在实验中使用PDA 测量了雾化液滴的尺寸和速度，图6 展示的是在距离喷射器出口轴向高度h＝10 mm 处液滴的轴向、径向与切向速度与实验数据的对比，可以看出程序模拟结果与实验相差不大，而在内回流区(IRZ)实验测得的轴向速度为负值，而程序模拟的轴向速度为零，猜测是由于在模拟时采用的是圆锥空心雾化方式，直径较小的液滴分布在中心，直径较大的液滴分布在边缘，当小液滴进入温度较高的区域时会快速蒸发，从而导致以上现象.

图6 热态工况下液滴轴向、切向与径向时均速度沿径向分布Fig.6 Radial distribution of time-averaged velocity of droplets in axial，tangential and radial directions under combustion condition

温度场可以直观地表征燃烧状态和火焰形态，图7 展示了流场在X＝0 截面瞬时温度云图，表明了沿流向方向温度的分布情况.图8 分别展示了在X＝0截面上实验测量的温度云图和AECSC 程序计算得到的时间平均温度云图，展示的均是图7 中粉色方框区域.实验中流场温度由OH 的绝对密度计算得到，而OH 的绝对密度由PLIF 和吸收同时测量，由于OH 浓度随温度迅速下降，该方法可获得的低温极限在1 500 K 左右，因此实验温度中色条覆盖范围是1 500 K 至2 200 K.图8 中程序模拟结果很好地再现了实验的V 型火焰，且高温区的温度值与实验结果很接近均在2 200 K 附近.模拟结果很好地还原了低温的内、外回流区，其中内、外回流区位置见图9，内、外回流区的大部分区域温度低于1 500 K 或在1 500 K 附近，与实验结果基本吻合，说明在内、外回流区没有发生明显的化学反应.

图7 X＝0截面瞬时温度云图Fig.7 Instantaneous temperature map at section X＝0

图8 X＝0 截面时均温度云图，包括AECSC 程序计算和实验结果Fig.8 Time-averaged temperature map at section X＝0，including AECSC program calculation and experimental results

空气气流流经内外旋流器，从旋流器出口高速喷出，形成了具有较大的轴向、径向和切向速度的高速旋流流场.两级旋流旋向相同，旋流产生的速度梯度可以使两股旋流得到充分混合，随着流向距离的增加，旋流沿径向逐渐向外扩张.图9 是计算得到的X＝0 截面的轴向速度云图，展示的是图7 中黑色方框区域.在轴向速度云图中可以清晰地看到回流区的边界(vmean＝0)，其中内回流区是由双旋流器的存在而产生的，因为它位于“V 形”喷射流之间，预期会影响燃料蒸气与预热空气的混合，因而具有重要的意义.外回流区靠近燃烧室壁边，由于存在突扩，在拐角处与主流形成一个漩涡区.混合气体的燃烧主要受外回流区和中心回流区的控制，在二者之间的剪切区进行燃烧，因此回流区的位置可以预测火焰产生的位置.

图9 X＝0截面轴向时均速度云图和等值线Fig.9 Time-averaged velocity map in axial direction and isolines at section X＝0

2.2 局部高温区分析

在30 万步的基础上继续模拟，以200 步为一间隔输出模拟结果，总共输出30 000 步，用于分析瞬时流场的局部高温区.计算时间步长为1.18×10-7s，统计非定常流动时间为0.003 54 s，即统计非定常流动时间开始于0.034 s 终止于0.037 54 s.

图10 是在X＝0 截面不同时间的瞬时温度云图，展示的也是图7 中粉色方框区域，图中的轮廓线均为燃料蒸气(C12H23)质量分数在10%的等值线.结合图9、10 中可以看出瞬时温度场一直在变化，燃烧主要发生在速度剪切层.

图10 X＝0 截面瞬时温度云图和燃料(C12H23)质量分数等值线Fig.10 Instantaneous temperature map and isolines of mass fraction of fuel (C12H23) at section X＝0

本文在分析时，将温度高于1 700 K 的区域定义为高温区，在轴向Y＝0 至Y＝0.13 m(二分之一轴向长度)的流场中温度高于1 700 K 区域占全部区域的15%左右.在图11 中X＝0 的中心截面高温区主要存在燃烧室的前段 0～0.075 m(燃烧室轴向长度为0.264 m)即Y/D＝0～5 部分，其中D 为中心旋流器出口直径，D＝0.015 m，Y 是轴向方向.在Y/D＝0～5区域，分析垂直轴向截面上的高温区变化规律.

图11 X＝0截面瞬时高温分布和流线Fig.11 Instantaneous high temperature map and streamlines at section X＝0

2.2.1 进动涡核对高温区的影响

从图11 中可以看出筛选出的瞬时高温区区域与瞬时流线中的涡旋结构位置接近，说明高温区的产生与剪切层中的涡团有关，瞬时流线图揭示了剪切层中大尺度旋涡的存在.

在图12 中发现剪切层左右两侧中流线所突出的漩涡中心并不位于同一高度，用黑色折线将旋涡中心连接，漩涡中心排列成锯齿状，这与旋转对称有关，此现象是旋流中PVC 形成的强烈标志.进动涡核是旋流流场中在高旋流数下典型的不稳定性现象，位于IRZ 并绕中心轴进动，由剪切层的不稳定性产生.PVC 是不稳定的大尺度旋转流动结构，可以促进流场的混合，周围大部分流体绕其作旋转运动，使之出现进动特征[14].

图12 X＝0截面瞬时流线(轴向速度染色)Fig.12 Instantaneous streamlines at section X＝0(colored by axial velocity)

图13 中展示的是X＝0 截面瞬时放热率云图，文中瞬时放热率的定义是单位时间内单位体积所释放的热量.公式为：

式中：i 代表化学反应中涉及的组分；Hi是对应组分的焓值；ki是对应组分的化学反应速率.

在图13 中发现瞬时放热率比较高的区域与筛选出的高温区前半段区域形状类似，此处也是漩涡比较密集的区域，结合图14 的轴向涡量图，表明在涡量比较高的区域，PVC 对周围流体的作用显著，加强了已燃的高温气体和未燃的高温气体掺混，促进燃烧和放热.而在图11 筛选出的后半段高温区域放热率很小，说明高温区域的形成主要原因不是燃烧放热，而是PVC 旋转流动产生的流体掺混.有研究表明放热增强会导致PVC 频率的增加和PVC 振荡的幅值减小[14].从图14 中也可以看出随着轴向距离的增加，涡量值减小，PVC 对周围流体的作用减弱.

图13 X＝0截面瞬时放热率云图Fig.13 Instantaneous heat release rate map at section X＝0

图14 X＝0截面瞬时轴向涡量图Fig.14 Instantaneous axial vorticity map at section X＝0

2.2.2 高温区的演化

图15 是相同时间不同方向的PVC 瞬时视图，采用旋流器出口附近的瞬时压力p＝0.426 MPa 等值面显示并以瞬时轴向速度染色.中央喷嘴的旋流器所产生的旋流运动驱动着旋涡芯，像施加的涡流一样，整个结构围绕自身涡轴顺时针旋转.

图15 进动涡核展示(p＝0.426 MPa 等值面，轴向速度染色)Fig.15 Precession vortex core visualized(iso-surface p ＝0.426 MPa，colored by axial velocity)

为了展示PVC 在时间上的演变，图16 进一步展示了相同的瞬时压力等值面(0.426 MPa)XZ 方向上在4 个不同时间的PVC，螺旋所包围的流动在轴向速度的作用下进行再循环，图16 中漩涡变化了一个周期，PVC 的频率在1 000～2 500 Hz.

图16 不同时间下进动涡核展示(p＝0.426 MPa 等值面，轴向速度染色)Fig.16 Precession vortex core visualized at different time(iso-surface p＝0.426 MPa，colored by axial velocity)

在分析垂直轴向截面上的局部高温区变化时，首先筛选温度高于1 700 K 的区域，并将局部高温区域旋转一圈的时间定义为一个周期，在垂直轴向截面Y/D＝0～5 上，高温区的变化频率均在 1 000～2 500 Hz，频率近似.图17 中展示了轴向Y/D＝1.4截面上t＝0.034～0.034 72 s 的局部高温区变化图，图中高温区从1 到6 旋转变化了一个周期，频率为1 388 Hz 左右.图18 中展示了轴向Y/D＝3.3 截面上t＝0.034 26～0.035 086 s 的局部高温区变化图，黑色实线圈出了截面上尺寸比较大的高温区，可以发现高温区旋转变化了一圈，频率为1 210 Hz 左右.在本文的统计中，从t＝0.34 s 开始高温区变化的第1 个周期如下：Y/D＝1 截面高温区变化频率在1 410 Hz 附近；Y/D＝2 截面高温区变化频率在1 280 Hz 附近；Y/D＝4 截面高温区变化频率在1 140 Hz 附近.截面上高温区的变化频率在统计的时间内会有所差异，但总体范围与PVC 的振荡频率相近，均在1 000～2 500 Hz，说明流场中的涡旋结构对燃烧有影响，且PVC 会影响附近流体，使之出现进动特征，导致燃烧场出现不稳定性，使垂直轴向截面的局部高温区域呈现周期性旋转.

图17 垂直轴向Y/D＝1.4截面上高温区旋转变化Fig.17 Rotation change of high temperature zone at section Y/D＝1.4 in vertical axis

图18 垂直轴向Y/D＝3.3截面上高温区旋转变化Fig.18 Rotation change of high temperature zone at section Y/D＝3.3 in vertical axis

高温区的旋转变化现象随着轴向距离的增加而越来越不明显，这是由于随着轴向距离的增加，涡量减小，附近流场的进动特征沿流向逐渐衰减，但变化频率近似，说明这些位置的局部高温区在空间上处在同一个大尺度流动结构中.气动力学引起的涡旋形成及破碎过程造成了燃烧室中燃烧的不稳定，并且PVC 会通过多种方式影响旋流火焰的稳定性，如PVC 会加强燃料和空气的混合，强化未燃和已燃气体的混合、卷曲、拉伸以及局部熄灭的化学反应区，也有可能与火焰的热声耦合振荡相互作用[16-17].

取剪切层内点X＝0，Y＝0.021 m，Z＝0.021 m，即X＝0，Y/D＝1.4，Z/D＝1.4.图19 展示了在统计的时间内瞬时压力脉动和放热率脉动变化情况，在时域图中可以看到其压力脉动和放热率脉动的波动具有一定的周期性，在频域图中，压力脉动和放热脉动对幅值贡献比较大的频率均在1 000～2 500 Hz，且在1 000～2 500 Hz 内频谱近似.燃烧场中的高温区呈现径向旋转变化，存在变化频率，并且压力脉动和放热率脉动的波动具有一定的周期性，说明在该情况下，燃烧出现不稳定性现象.