旋流预混燃烧室火焰动力学奇异谱分析研究

杨 旸,余志健

(1.中国科学院工程热物理研究所 南京未来能源系统研究院，南京 210000；2.中国科学院工程热物理研究所 先进燃气轮机实验室，北京 100190；3.中国科学院大学，北京 100190)

重型燃气轮机主要采用贫预混燃烧技术降低氮氧化物排放。贫预混燃烧火焰传播速度降低，同时燃烧室上取消了气膜孔，减弱系统声学阻尼。当热释放率波动与压力波动相位差小于90°，并且注入系统的能量大于系统的耗散时，波动随时间不断增大，最终达到极限环状态或过载，出现燃烧不稳定问题[1]。

在燃烧室设计阶段需对火焰动力学分析以解决上述问题。目前主要有两种方法：一是建立燃烧室闭环模型，利用快速傅里叶变化分析自激热声数据特性，但该方法需要燃烧系统详细的阻抗边界条件；二是将燃烧系统作为线性或弱线性的黑箱模型，利用脉冲响应，在开环状态下计算火焰传递函数[2]。后者考虑纯火焰的特性，在火焰和燃烧室设计上有更大的自由度，更适合工业应用[3]。

获得和分析热释放率数据是应用上述方法的起点。由于预混燃烧中主要为大尺度湍流使火焰前缘发生褶皱，因而大涡模拟适用于预混燃烧模拟[4]。然而，由于大涡模拟时间步长小、网格尺度细，求解时间长，结果难以为设计工作提供足够快的反馈。此外，模拟获得的热释放时间序列有时间短、非平稳和趋势被噪声所掩盖等特性。采用传统的火焰传递函数法，计算的响应特性结果差，难以分析火焰未来趋势，例如是否熄灭。本文将一种数据驱动算法——奇异谱分析(Singular Spectrum Analysis,SSA)应用于上述热释放率时间序列分析中，以进一步提取火焰动力学特性。

奇异谱分析最初由Colebrook[5]提出，用于描述海洋浮游动物浓度变化。而后Broomhead和King从一个测量的时间序列中重构了非线性动态系统的吸引子[6-7]。奇异谱分析方法对短时间、有噪声、混沌时间序列具有良好的计算效果[8]。利用该方法的降噪能力，可有效重构自激管式燃烧室由热声不稳定引起的稳态和非稳态压力时间序列[9]。与传统方法相比，奇异谱能够提取高度非线性的特征。采用动态模态分解分析和奇异谱分析相结合的方法，研究了高压预混燃烧室燃烧不稳定性与火焰涡相互作用的关系。结果表明，纵向模态与轴对称涡旋脱落有关[10]。

目前采用奇异谱方法对火焰动力学研究较为缺乏。本文对一个旋流预混燃烧室进行了大涡模拟，通过在大涡模拟中添加进口速度激励以获得热释放率响应时间序列。而后采用奇异谱分析方法重构并降解了该序列，定义了传递路径响应参数(激励与降解后热释放率之间的传递函数)，以分析火焰响应特性。此外，奇异谱模态结果还与动力学模态降解(DMD)结果比较，以形象揭示响应特性起源。本文致力于探索奇异谱方法在燃烧室设计阶段对火焰动力学的分析，弥补现有方法的不足，以指导燃烧稳定的低污染燃烧室设计。

1 数值方法

1.1 大涡模拟

1.1.1 几何及边界条件

对旋流预混燃烧室进行大涡模拟以获得热释放率时间序列。大涡模拟计算域如图1所示。气流方向为从左到右。图1(b)为采用的8叶片轴向旋流器。叶片吸力面设有两个燃料孔，压力面设有一个燃料孔。几何尺寸如表1所示。采用中心体半径R1作为归一化尺寸，掺混区长度为l，中心体末端半径为R1_TE，扩张角为α。旋流气流随扩张角向燃烧室扩散时，当气流离心力与径向压力梯度平衡时，根据涡破碎理论形成回流区[11]。回流区的形成通常由标准旋流数控制，旋流器速度基旋流数为0.4。边界条件如表2所示。当前几何条件下，当量比为0.514。

(a)模型燃烧室

表1 经R1归一化几何参数

表2 边界条件

1.1.2 网格划分和网格分辨率

网格采用star ccm生成的多面体网格，边界层采用三层棱柱层网格。计算域为45°扇形，采用周期性边界条件降低计算成本。燃烧室不同区域采用不同空间分辨率的网格尺寸。旋流器、掺混区、燃料孔及燃烧室上游采用更小的网格尺寸以捕捉流动和火焰细节。整体网格划分和局部位置网格如图2所示，总网格数为460.5万。

图2 全局及局部网格划分

Fluent在大涡模拟求解中采用隐式求解器，无网格无关性检验，因为无限减小网格尺寸将导致求解趋向DNS。本文采用基于能量指数法的后检验来判断网格在LES计算中是否足够好。即当求解的湍流动能占总湍流动能85%时，认为网格求解精度达到[12]。由于应用边界中心差分格式进行空间离散，因此总湍流动能中忽略由数值格式引起的能量耗散。中心平面后估计指数M分布云图如图3所示，可以看出大部分区域能量占比都大于85%,由于本研究集中在火焰区，出口小部分区域解析度不够也可满足本文研究目的。

图3 中心平面后估计指数M分布云图

1.1.3 数值方法

生成的网格导入ANSYS fluent 18.0进行数值计算。大涡模拟采用WALES亚网格尺寸模型[13]，燃料和空气进口无湍流波动设置。瞬态求解采用二阶差分格式。CFL数在全域内保证小于1。采用部分预混燃烧模型，该模型通过求解进程变量的平均输运方程模拟预混燃烧反应。控制方程如下所示：

(1)

(2)

ut=Au′Da1/4

(3)

甲烷燃烧采用GRI 3.0反应机理[15]，比热容与温度的关系采用分段多项式模型。模拟中采用绝热边界条件。

1.2 火焰传递函数和传递路径响应

(4)

传递函数计算步骤如下：1)计算输入信号的自功率谱；2)计算输入和输出的互功率谱；3)在拉普拉斯域上获得输入与输出的幅值和相位值。

1.3 奇异谱分析

奇异谱分析是处理平稳或非平稳时间序列的一种无模型方法。时间序列在无先验假设前提下，可分解为趋势、周期及噪声等可解释分量[21]。首先将燃烧室大涡模拟获得的热释放率时间序列

(6)

(7)

窗长度为L(L需小于n/2)，奇异谱分析算法步骤如下[21-22]：

(1)嵌入

首先获得轨迹矩阵Y(Hankel矩阵)。

(8)

窗长度L越大，原时间序列信息损失越大；L越小，计算量越大，存储所需内存越大。

(2)奇异值分解

对矩阵Y进行奇异值分解，Y=UΣVT。其中UL×L,V(n-L+1)×(n-L+1)为左右奇异向量，ΣL×(n-L+1)为对角线矩阵，其中特征值λi按照降序排列。根据谱分析理论，矩阵Y可用前d阶模态重构，如公式(9)所示。

(9)

(3)重构

(4)筛选

由于奇异值分解的特征值已经以递减的方式排序，因此可选择前k个模态进行分析。引入回归模型中协方差R2值，R2越接近1，表明自变量对因变量的解释越有效。

2 结果与分析

2.1 火焰结构和火焰传递函数

中心平面时均和均方根(RMS)进程变量(Progress Variance，PV)分布云图如图4所示。从时均云图可看出，中心回流区PV值为1，角回流区PV值小于0.8，在火焰内层下游，PV也小于0.8，火焰面有效形成。进程变量RMS云图表明火焰内外层之间有强烈的混合，而中心回流区与火焰内层混合作用弱。

图4 中心平面时均(上)和RMS(下)进程变量(PV)云图

激励响应大涡模拟过程，热释放率采样频率为10 000 Hz。图5为从大涡模拟中采集的激励后热释放率时间序列，可以看出热释放率呈下降趋势。Kwiatkowski-Phillips-Schmidt-Shin (KPSS)测试表明该序列趋势上是非平稳的，但序列的一阶导数是平稳的，因而可进行奇异谱分析。图6为计算的火焰传递函数幅值和相位。传递函数幅值在零频率附近有较高的幅值，较高的幅值掩盖了其它频率带的响应信息，因而火焰准确的响应频率信息难以从火焰传递函数中获得。需发展新方法如奇异谱分析进一步揭示热释放率内在响应信息，以指导燃烧室设计。

图5 激励后热释放率时间序列

图6 火焰传递函数

2.2 热释放率时间序列奇异谱重构

窗长度L选为1 200，图7为各阶重构矩阵对轨迹矩阵的单独和累计能量占比。可以看出采用前20阶模态的重构矩阵累计能量占原矩阵能量的80%以上。第7阶之后的单独模态能量占比都保持较低水平。采用前20阶模态重构后，累积协方差达到0.9以上，表明前20阶模态重构原始热释放率时间序列有效。

图7 各阶重构矩阵Y占轨迹矩阵的单独和累计能量

图8为前k(k=10,20)阶奇异谱模态重构后的热释放率时间序列与原始时间序列对比。可以看出采用前20阶奇异谱模态重构的时间序列可有效的跟踪原始时间序列，并可准确捕捉热释放率波动特性。后续主要以前20阶奇异谱模态进行分析。

图8 采用前k阶奇异谱模态重构后的热释放率时间序列和原始序列对比

2.3 奇异谱模态和传递路径响应

图9为1、4、9及12阶热释放率奇异谱模态时间序列。图10至图14为热释放率按响应特性分类的奇异谱模态传递路径函数幅值。从热释放率奇异谱模态幅值可以看出，模态可以划分为趋势模态和对特定频率的四种响应模态。趋势模态反映热释放率随时间变化的总体趋势，响应模态反映热释放率的动态响应特性。响应特性分为32 Hz、56 Hz和112 Hz、200 Hz以及128～136 Hz。图9为第一、三和四种响应类型中的首阶模态和第二种响应类型的第二阶模态(首阶和第二阶模态能量占比较大)的时间序列。

(a)1阶

1阶和2阶模态反映热释放率的趋势特性。图10为1阶和2阶模态传递路径响应幅值。从图9(a)可以看出，1阶模态的振幅随时间不断衰减，在0.15 s后达到低值并保持。传递路径响应幅值曲线可以看出在零频率处保持较高增益，表明热释放率时间序列趋势演变是低频的。

图10 1阶和2阶趋势模态火焰传递路径响应

3、4、5和6阶奇异谱模态对32 Hz频率有响应。图11为3、4、5和6阶模态传递路径响应幅值。从图9(b)，可以看出4阶模态的时间序列呈波动状，振幅随时间不断增大。从频谱特性看出，该部分各模态的特征频率在30 Hz左右，与响应特性保持一致。

图11 32 Hz响应模态(3、4、5和6阶)传递路径响应

9、10、11、18和19阶模态同时有56 Hz和112 Hz响应特性。图12为上述各阶模态传递路径响应幅值。从图9(c)所示，9阶模态时域上也为波动状，频谱上主频在110 Hz左右。传递路径响应幅值有56 Hz和112 Hz两个主频。

图12 56 Hz和112 Hz响应模态(9、10、11、18和19阶)传递路径响应

12和13阶热释放率模态对200 Hz有响应特性。图13为12和13阶模态传递路径响应幅值。从图9(d)所示，12阶模态时域上振幅不断增大，到最大值后衰减。传递路径响应曲线表明主频有200 Hz。

图13 200 Hz响应模态(12和13阶)传递路径响应

16阶模态对频率带128～136 Hz有响应。图14为16阶模态传递路径响应幅值，响应峰带宽较长。如图6所示，采用传统的火焰传递函数方法，难以解析出设计阶段大涡模拟获得的短时、非稳态和含噪声的非稳态热释放率数据响应特性。采用奇异谱分析处理这类数据，通过各阶模态传递路径响应特性，可有效分辨出热释放率的响应特性。

图14 128～136 Hz响应模态(16阶)传递路径响应

进一步对奇异谱各阶模态的时间序列进行分析，图15为各响应频率代表性模态时间序列递归分析图。可以看出对于1阶和2阶的趋势模态(无特定响应频率)，模态递归图为黑白图，无明显特征规律。对于具有32 Hz响应特性的3阶和4阶模态，模态递归图为黑白相间的棋盘状，且某一区域该特征极为明显。对于具有200 Hz响应特性的12和13阶模态，时间序列模态递归图也呈棋盘状，且棋盘尺寸更为密集，代表响应频率越高。

图15 各响应频率代表性模态时间序列递归图

表3为动力学模态降解(DMD)处理该热释放率时间序列获得的各阶模态特征频率和增长率值。与奇异谱分析方法获得的模态特征频率对比，可以发现奇异谱中的32 Hz、112 Hz和200 Hz左右的响应频率也可被DMD方法的11阶、13阶和18阶模态捕获。图16为动力学模态降解法获得的18阶模态(对应31 Hz响应)进程变量的分布云图，可以看出该模态下，火焰内外层有剪切运动，角回流区和中心回流区有混合作用。

表3 动力学模态降解法获得的各阶模态频率及增长率

图16 动力学模态降解法获得的18阶模态(对应31 Hz响应)

3 结论

(1)燃烧室设计阶段，对于某些工况下采用大涡模拟获得的短时、非稳态和含噪声的热释放率时间序列，采用传统火焰传递函数法获得的频率响应信息易被低频高幅值段掩盖。

(2)对热释放率时间序列采用奇异谱分析方法处理，重构后热释放率时间序列可有效吻合原始时间序列。

(3)奇异谱分析方法获得的各阶模态火焰传递路径响应参数表明，模态可分为趋势模态和响应模态；该火焰有4种典型的频率响应特性，分别为32 Hz、56 Hz和112 Hz、200 Hz以及128～136 Hz响应，揭示了火焰内在频率响应特性。

(4)动力学模态降解法捕获的31 Hz、98 Hz和187 Hz频率与奇异谱模态获得的四种响应特性基本吻合。31 Hz下DMD模态形状可以看出火焰内外层有剪切运动，回流区有混合作用。