燃烧不稳定是指燃烧器中燃烧的周期性振荡现象, 其伴随着压力和速度的大幅度自激振荡, 这一振荡带来的循环负荷会造成结构过度振动和燃烧室壁面的热量过度传递等问题, 对结构组件造成致命的破坏。长久以来, 燃烧不稳定一直都是火箭发动机^[1]、航空发动机^[2]、燃气轮机^[3-4]等许多高性能燃烧装置面临的重大挑战, 引发了各方面的深入研究^[5-7]。热声不稳定是燃烧不稳定最主要的一种形式, 其本质是热源的放热与声场环境存在正反馈机理^[8]。在拥有特定声学特性的燃烧装置中, 系统的热释放取决于压力和速度波动, 热释放波动又会作用于燃烧室声场, 从而形成反馈回路, 使得系统产生不稳定性。Rijke管是典型的热声不稳定系统, 也是研究燃烧与声场相互作用机理及控制的基础, 因此许多研究都是针对Rijke管来开展的^[9-16]。目前为止, Rijke管热声不稳定现象的研究理论主要分为线性理论、非线性理论和非线性动力学理论。现有线性理论研究表明, 在Rijke管1/4长度处放置热源, 系统是线性不稳定的。Rott等人^[17]利用各阶模态的经典线性稳定性分析, 将系统视为网状模型, 每一个元素用线性传递函数建模, 通过检测特征值来判断系统的稳定性。利用线性化模型方法的研究, 对热声系统发生不稳定的条件有了一定的理解。利用非线性理论对Rijke管热声系统中极限环振荡的非线性效应也已有所研究。Heckl^[18]通过对Rijke管热声系统非线性效应进行研究, 建立了一个可以预测极限环振幅的经验模型。Hantschk等人^[19]研究表明极限环振荡由放热元件释放到流场的热通量的非线性决定。Matveev等^[20-21]采用能量平衡方法证明了系统的平衡状态处于热声能量的输入和声场损失平衡的时刻, 这一结论再次证实不稳定放热的非线性是极限环振幅的主导因素。

非线性动力学理论可对热声系统的非线性系统特征进行更深入的刻画、分析和诊断, 利用非线性动力学来研究热声不稳定是目前研究的一大热点。Jahnke和Culick^[22]最早将现代动力系统理论引入燃烧不稳定研究中以分析非线性不稳定燃烧。他们提出了一种可对各稳态下极限环的解进行体系化高效计算的延拓算法, 并由此得到系统在各稳态和极限环下的稳定性, 通过分岔分析判定得到不稳定发生点。Balasubramanian等^[23]针对水平Rijke管建模, 研究了非正则性和非线性效应在简单热声系统中的作用和影响。研究发现, 控制方程的非正则性会使极限环振荡瞬态增长, 触发系统的非线性行为。Subramanian等^[24-25]使用DDE-BIFTOOL分析了Rijke管热声系统中亚临界分岔性质的问题, 描述了带精确滞后的热声系统模型, 采用多尺度方法对Rijke管模型进行弱非线性分析, 重点讨论了Rijke管热声不稳定的亚临界Hopf分岔特性。

非线性动力学理论可深入描述和解释热声不稳定现象, 但目前利用非线性动力学对Rijke管热声系统非线性行为的研究相对匮乏。尽管对Rijke管在各稳态的极限环振荡基本清楚, 但对其在多参数变化下的稳定性边界和双稳态区域尚不明晰。此外, 系统在双稳态区域内存在的触发现象, 尚未有深入研究。为了进一步探究Rijke管热声系统在不同系统参数作用下的稳定性边界, 获得两参数变化下的双稳态区域变化特点和系统在双稳态区域内触发现象的特性, 本文首先建立水平Rijke管热声系统简化模型, 利用Galerkin方法对其一维声场控制方程进行求解, 结合非线性动力学理论, 探究Rijke管热声系统的稳定性边界和稳定性区域, 随后对双稳态区域内的触发现象进行分析, 获得系统的临界触发值, 并研究触发值对极限环振荡的影响。本文研究工作可加深对Rijke管热声不稳定机理的认识, 为火箭发动机等高性能燃烧装置的热声不稳定研究提供重要理论基础。

1 Rijke管模型的建立与求解

本文所研究的水平Rijke管模型如图 1所示。

在该模型中, 管长L是固定不变的, 圆管半径为r, 热源位置x_h的值可以通过改变加热器在管道内的位置来改变, 为了方便以后的建模计算, 我们定义热源相对位置为x_f=x_h/L。在实验条件下, 电热丝加热装置可以通过改变通过的电流改变加热强度。同时为了达到建模的目的, 热释放区域被认为是紧凑型的, 管内流动为绝热环境下的无黏流动。由于Rijke管的管长远远大于声波的波长, 故可将Rijke管热声振荡问题简化为一维问题。忽略平均流量与温度梯度的影响, 控制管内声场的动量和能量方程组如下:

(1)

(2)

为了简化和方便分析, 对(1)式和(2)式进行无量纲化处理, 得到无量纲化的系统控制方程组:

(3)

(4)

式中, 无量纲化尺度如下:

(5)

(1)~(5)式中, x是沿轴向方向的距离, t是时间, 为管内平均压强, u为管内平均流动速度, u′为声振速度, p′为声波压强, γ是介质的比热比, 声速为c₀, Ma为管内平均流马赫数, ρ为介质密度。此外, ξ为阻尼系数, 为单位面积上电加热丝的热释放速率, x_f为热源相对位置。符号‘^~’和‘^-’分别表示有量纲量和平均量。

本文采用的热源模型为Heckl经验模型^[18]。热量从热源到声速波动响应由相关性分析来量化, 以衡量从加热器到其周围流动的准稳态传热, 并引入一个时间延迟τ以考虑热传导的惯性, 热释放速率表达式如下

(6)

式中, L_w和d_w代表电热丝的长度和直径, T_w为电热丝的温度, T代表周围空气的平均温度, S为管道截面积, λ为空气的热传导率, 为标准狄拉克分布, c_v表示单位质量空气的恒容比热容, ρ为气体平均密度。

热惯性是调制热释放与声振荡间关系的重要参量。在本文所选热源模型中, 由于热惯性的存在, 传热和流动之间存在时间延迟τ。Lighthill指出在雷诺数大于10并且频率小于时, 时间延迟的大小与电热丝直径和管内平均流速有关, 无量纲时间延迟τ可用(7)式来表示, 其中d_w为电热丝直径, (c₀/L)为附加因子以使时间延迟τ无量纲化。根据不同的热源设置和流速, 无量纲时间延迟τ可在0~1左右变化。

(7)

使用Galerkin投影和模态展开的方法可以求解Rijke管热声振荡问题。Galerkin方法的原理是任何域中的函数都可以表达成该域内基函数的叠加。基函数的选取原则是它们须满足边界条件, 且基函数的选择不是唯一的。通过将声学变量扩展到基础函数系, 可以把Rijke管热声模型一维无量纲控制方程(3)~(4)简化为常微分方程组^[27]。本文中选择的基函数是线性化后的系统自共轭部分的特征函数。因此, 速度场和压力场可以表示成管道固有模态的形式:

(8)

(9)

式中, k_j=jπ为j阶管内声模态的无量纲波长。在极限N→∞时, 这些基函数形成一个完整的基础函数系。本文中水平Rijke管热声模型两端的声边界条件为开放声边界, 显然, 本文的速度和压力信号展开满足声边界条件的要求。

将(8)式和(9)式代入到(3)式和(4)式中, 并沿基础函数系投影, 得到:

(10)

(11)

式中, ξ=2ξ_jω_j, w_j=jπ为第j阶管内声模态的无量纲频率。

本文建模中所考虑的频率相关阻尼是由管端不完全反射造成声损失以及边界层损失造成的。这种阻尼先前已被Subramanian^[25]和Juniper等^[26]在相同的设置下对热声不稳定建模。阻尼ξ_j可以用(12)式表示, c₁和c₂为阻尼系数, 在某一实验条件下为保持不变的常数。在实验条件下, 阻尼大小可以通过改变管末端条件来实现。在本文计算中, 阻尼对系统响应的影响通过改变阻尼模型中阻尼系数c₁的值来实现。无量纲加热功率K由(13)式给出：

(12)

(13)

式中, L_w为电热丝长度, T_w为电热丝温度, T为周围空气平均温度。

在极限N→∞条件下, (7)式和(8)式这些基函数形成一个完整的基础函数系, 但是考虑到计算的可行性, 模态数目只能设置为有限个。精确捕捉到系统的线性和非线性行为所需Galerkin模态数目被称为模态收敛。通过Matlab编程实现对(10)式和(11)式的数值求解, 并比较了不同数目Galerkin模态下的系统演化。对于t阶声学模态, 增加一阶Galerkin模态所形成的误差不超过1.4%^[28], 其误差在可接受范围内。因此, 在本文以后的所有计算中, 采用十阶Galerkin模态的模型来进行线性和非线性稳定性分析以确保收敛性。

在本文计算中, Rijke管热声模型的几何参数和物理特性与之前典型实验条件设置是一致的^{[20, 29]}。特别地, L=1.0 m, r=0.025 m, γ=1.4, λ=0.032 8 W/(m·k), ρ=1.025 kg/m³, C_v=720 m²/(s²·K)。相应地, 根据典型Rijke管实验条件下, 无量纲加热功率K在0~3的范围内变化, 无量纲时间延迟τ的变化范围为0~1左右, 阻尼系数c₁的取值范围为0~0.5, 无量纲处理后的热源相对位置x_f可在0~0.5内变化。

2 结果与讨论 2.1 稳定边界和双稳态区域

线性稳定性分析可研究无穷小扰动对系统稳定性的影响和系统在稳定状态附近的局部稳定性变化。非线性稳定性分析则是研究有限振幅扰动对系统稳定性的影响, 用于表征系统的渐近状态。通过对系统进行线性稳定性分析和非线性分析, 可得到系统参数间的线性稳定边界和非线性稳定边界。线性稳定边界和非线性稳定边界将参数平面划分成3个区域:全局不稳定区域, 双稳态区域和全局稳定区域。参数在线性稳定边界内的区域, 对于任意小的初始扰动, 热声系统最终都发展为振荡状态, 该区域称为全局不稳定区域。参数在非线性稳定边界以外区域, 对于任意幅值的扰动, 系统都最终趋于稳定状态, 称为全局稳定区域。双稳态区域位于线性稳定边界和非线性稳定边界之间。

Rijke管热声系统的主要控制参数为无量纲加热功率、热源位置、阻尼系数和时间延迟。通过对Rijke管热声模型进行线性和非线性分析, 分别获得了无量纲加热功率K、热源相对位置x_f、阻尼系数c₁和无量纲时间延迟τ间的线性边界、非线性边界和各稳定区域。

Rijke管热声系统中的热源为电热丝加热, 在实验条件下可以通过改变电流大小来改变加热器功率。在本文的数值计算中, 通过改变无量纲加热功率K值以及无量纲时间延迟τ值来探究二者对系统稳定性产生的影响。在非线性动力学理论中, 无量纲加热功率的增大代表着系统驱动力的增加, 从而使得系统更加不稳定。无量纲加热功率K与无量纲时间延迟τ间的稳定边界如图 2所示。

可见, 全局稳定区域位于非线性稳定边界外, 全局不稳定区域位于非线性稳定边界内。当无量纲加热功率K值足够小时, 系统的稳定状态不受时间延迟的影响, 恒处于全局稳定状态。全局不稳定区域内, 对于任意小的扰动, 系统最终都趋近于不稳定振荡状态。在非线性稳定边界和线性稳定边界之间的双稳态区域内, 系统的最终状态取决于外界扰动的大小。当扰动大于某一阈值时, 系统才会产生热声振荡; 相反, 系统将呈现稳定状态。在τ < 0.5时, 双稳态区域随着系统无量纲时间延迟τ的增加而逐渐增大。可以发现, 加热产生的非线性影响在时间延迟较大时更加严重。在无量纲时间延迟取某一定值(如τ=0.4)时, 系统在无量纲加热功率较小时处于全局稳定状态, 随无量纲功率不断增大, 系统状态变化为全局不稳定。通过无量纲功率和无量纲时间延迟的稳定区域分布也可以证明增大无量纲加热功率K使得Rijke管热声振荡更加不稳定, 这与杨亚晶等^[29]的实验结果相一致。

为了研究阻尼对系统稳定区域的影响, 我们通过改变阻尼模型中的一个阻尼系数c₁来实现。在实验条件下, 系统的阻尼是通过改变管道的末端条件来实现^[17]。根据非线性动力学理论, 增加阻尼可以增强系统的稳定性。图 3为阻尼系数和无量纲时间延迟件的稳定边界和稳定区域分布图, 可以发现, 阻尼系数c₁与无量纲时间延迟τ间的双稳态区域为带状。随着无量纲时间延迟τ的增大, 双稳态区域先增加, 在τ=0.5处达到最大值, 然后逐渐变小。系统全局不稳定区域处于阻尼系数较小段, 而当阻尼系数较大时, 系统为全局稳定状态。在c₁>0.42时, 系统的稳定状态不再受时间延迟的影响, 并一直处于全局稳定状态。如预期的一样, 增加阻尼对系统的不稳定振荡有抑制作用。

热源位置对系统的动力学有显著影响, 热源位置的变化通过把电热丝加热器放置在距管口不同位置处来实现。不同于无量纲加热功率、阻尼系数二者和无量纲时间延迟间的稳定边界分布, 热源位置的变化以一种特殊的方式影响着系统稳定区域的分布。在系统其他参数一定时, 热源相对位置x_f与无量纲时间延迟τ间的系统稳定区域如图 4所示。系统全局不稳定区域位于双稳态区域之间。系统双稳态区域为2个带状区域, 且关于x_f=0.25近似呈对称分布。另外, 系统的双稳态区域随无量纲时间延迟τ的增加呈现先减小后增大的变化规律, 与阻尼系数c₁和无量纲时间延迟τ间的双稳态区域分布正好相反。全局不稳定区域被分割为两部分, 在0.06 < x_f和x_f>0.44时, 系统恒处于全局稳定区域, 不受时间延迟的影响。由热源相对位置和无量纲时间延迟间的稳定区域分布, 可以发现, 热源相对位置对系统稳定区域的影响具有近似对称性。

2.2 触发

对于一个非线性系统, 存在一些初始条件使振荡逐渐衰减, 同时也存在初始条件, 使得振荡增长形成极限环, 这种现象称为“触发”。简单热声系统中, 触发现象发生的必要条件:系统要有不稳定周期解, 亚临界Hopf分岔点和折叠分岔点。本文所研究的Rijke管热声系统满足该必要条件^[30]。在其双稳态区域内, 系统可能发展为极限环振荡或最终达到稳定状态, 触发现象就发生在系统的双稳态区域内。我们探究了在热源无量纲加热功率K、热源相对位置x_f和阻尼系数c₁双稳态区域的触发现象, 并得到了系统触发所需的最小扰动, 即“临界触发值”。

无量纲功率K双稳态区域内, 系统在不同扰动下的非线性行为随时间的变化如图 5所示。当t < 100时, 系统保持稳定状态, 在t=100时, 系统受到U₁=0.2的扰动后振荡逐渐衰减, 并最终达到稳定状态, 而当t=600时, 系统受到U₁=0.45的扰动后振荡逐渐增长, 并最终达到极限环状态。因此扰动值的大小对系统最终的稳定状态产生不同的影响。表 1给出了对应于图 5参数设置下, 系统在无量纲加热功率K变化时的临界触发值。U₁和P₁分别为速度扰动值和压力扰动值。可以看出, 随着无量纲加热功率K的增大, 临界触发值U₁和P₁均呈减小趋势, 这与Matveev等^[20]的实验发现规律一致。

图 6和图 7分别为阻尼系数c₁双稳态区域和热源相对位置x_f双稳态区域内, 系统在不同扰动下的行为变化。可以发现在系统较小扰动后, 虽然出现振荡, 但振荡随时间逐渐衰减至稳定状态, 而系统在受到超过触发临界值的扰动后出现振荡增长, 并最终形成极限环。表 2给出了对应于图 6参数设置下, 系统在阻尼系数取不同值时的临界触发值。可以发现, 随着阻尼系数的增大, 系统需要更大的扰动才能被触发形成极限环振荡。表 3为对应于图 7参数设置下, 当热源处于不同位置时, 系统的临界触发值。在0.11 < x_f < 0.14时, 热源距离管口越远, 系统的临界触发值越小; 在0.37 < x_f < 0.4时, 热源与管口距离增大, 系统触发需要更大临界值。在热源相对位置x_f的2个双稳态区域内, 系统的临界触发值差异较大, 呈现相反的增长趋势。

针对K, x_f和c₁参数的双稳态区域, 计算了系统对U₁和P₁扰动的触发情况(见表 1~3)。分析可知, 系统的临界触发值P₁具有与U₁一致的变化规律, 说明Rijke管热声系统对压力扰动和速度扰动具有相同的响应趋势。

2.3 极限环

在非线性热声系统中, 极限环是最典型和最具研究意义的周期振荡。在2.2节中, 讨论了在控制参数的双稳态区域内, 系统临界触发值的变化规律。当扰动小于临界触发值时, 振荡会逐渐衰减, 系统最终达到稳定状态; 而当扰动高于临界触发值时, 系统被触发后, 振荡会逐渐发展为稳定的极限环。

为了探究双稳态区域内扰动值对极限环的影响, 以无量纲加热功率双稳态区域的触发现象为例进行分析。在无量纲加热功率双稳态区域内选取不同的扰动值, 观察系统行为随时间的演化。由图 8的内插图可以发现, 当扰动值较小时, 系统振荡逐渐增长, 最后达到稳定极限环; 而当在较大的扰动值激励下, 系统热声振荡初期先逐渐衰减, 而后达到稳定的极限环状态。分别对系统的极限环振荡进行频谱分析, 获得其频谱图, 如图 8所示。在2种不同幅值扰动的情况下, 系统极限环的频率和振幅都保持不变。可见在双稳态区域内, 扰动幅值对系统最终达到的稳定极限环状态不会产生影响。值得注意的是, 扰动幅值的变化虽然不会影响最终的极限环, 但会改变系统初期的振荡状态。系统在U₁=0.4扰动下达到稳定振荡所需时间相比U₁=0.8延长了约3倍, 表明增大扰动会缩短热声振荡达到稳定的时间历程。

3 结论

本文对Rijke管热声系统进行建模, 利用Galerkin方法对控制方程进行数值求解, 分析了系统非线性动力学特性, 得出如下结论:

1) 获得无量纲加热功率K, 阻尼系数c₁、热源相对位置x_f和无量纲时间延迟τ之间的线性、非线性稳定边界以及双稳态区域。热源相对位置x_f的系统稳定性区域近似关于x_f=0.25对称, 阻尼系数c₁的双稳态区域在τ=0.5时达到最大值。

2)无量纲加热功率K对系统的非线性影响在延迟时间较大时更加严重。无量纲加热功率K的增大会促使系统产生热声振荡, 而阻尼系数c₁的增加对系统的热声振荡有抑制作用。

3)获得系统在无量纲加热功率K, 阻尼系数c₁和热源相对位置x_f取不同值时的临界触发值。无量纲加热功率K的增加和阻尼系数c₁的减小均会导致系统触发值的增大, 触发值在热源相对位置x_f的2个双稳态区域内具有相反的变化趋势。

4)在双稳态区域内, 系统对压力扰动和速度扰动具有相同的响应趋势, 扰动值的大小不影响系统最终的极限环状态, 但会改变系统初期的振荡状态。系统在U₁=0.4扰动下达到稳定振荡所需时间比U₁=0.8延长约3倍, 增大扰动幅值会缩短热声振荡达到稳态的时间历程。