基于Lyapunov稳定性原理的双Rijke管耦合控制策略

Coupling Control Strategy for Rijke Tube Systems Based on Lyapunov Stability Principles

摘要

为控制表征为场量高幅值脉动的热声不稳定现象，基于Lyapunov稳定性原理开展Rijke管双管耦合控制研究. 通过线性化和矩阵化耦合系统控制方程获得关键系数矩阵，依据特征值实部预测耦合系统实现振幅死亡的控制参数范围. 研究表明，最大瞬态能量增长率与特征值实部最小值密切相关，基于稳定性原理的双管耦合控制策略能有效预测控制结果. 特征值实部最小值越大，最大瞬态能量衰减速率越快，但使用热源强度比大于1.00的耦合管则会造成被耦合管系统瞬态能量的增长.

Abstract

To effectively control thermoacoustic instability represented as high-amplitude oscillations， the coupling control method for Rijke tube based on Lyapunov stability principles is explored. The key coefficient matrices are obtained through linearization and matrix control formulation of the coupled system. Predictions of the control parameter range strategy for achieving amplitude death of the coupling control system are based on the real part of the eigenvalues. The research demonstrates that the maximum transient energy growth rate is closely associated with the minimum real part of the eigenvalues， validating the efficacy of the coupling control strategy based on the Lyapunov stability principles. A larger minimum real part value results in a faster decay rate of transient energy， but employing a coupling system with the ratio of heat intensity larger than one leads to a transient energy growth.

Graphical abstract

关键词

热声 / 稳定性 / Rijke管 / 双管耦合控制 / 振幅死亡

Key words

thermoacoustics / stability / Rijke tube / coupling control strategy / amplitude death

热声不稳定现象广泛存在于热动力机械系统中，对装置造成一系列负面影响^［1］. 如造成装置过度振动，导致部件的结构损毁和电控系统的失效；其伴随的放热波动，会破坏装置的正常燃烧循环而引起燃烧过程的熄火或回火，造成较大的循环变动，燃烧效率降低，大量噪声以及过度传热^［2］；即使是低振幅的振动也会导致发动机部件由于长时间磨损而疲劳失效；增加质量、动量和能量输运，改变发动机的推力水平，使其偏离设计值^［3］. 因此，在产品设计早期阶段，确定热声不稳定性的可能性并设计有效的热声不稳定性控制方法是一项重要且艰巨的任务^［4-5］.

研究热声不稳定过程及其影响机制，旨在开发有效的控制手段，来减轻高幅值压力、速度的振荡对热动力机械带来的负面影响. 因此，需要以某种方式打破非定常放热和流动扰动之间的动态耦合，以抑制、防止燃烧的不稳定. 通常，依据是否向系统中输入额外的能量，分为主动控制和被动控制策略. 主动控制方面，Dharmaputra等^［6］基于贝叶斯优化提出的自适应控制策略具有良好的工况适应性，但由于要安装反馈控制装置并针对系统设计确切的执行器，过程复杂而面临诸多困难，其相关研究主要集中在实验室阶段^［7］. 被动控制方面，Henderson等^［8］研究的电声控制方法，具有结构简单、可靠性高等优势，但同时也存在改装复杂、设计成本高、设计周期长、泛化性能差的弊端.

近期提出的双管耦合控制方法，具有结构简单、控制简便的优势，通过调节针阀和乙烯基管的长度可实现两个热声不稳定系统的同时控制. 因此，本文基于Lyapunov稳定性原理和瞬态能量分析方法，预测双管耦合控制策略.

1 Rijke管双管耦合系统理论模型

1.1 Rijke管系统

由于热动力机械系统涉及热力学、声学、流体力学和化学动力学等多物理场的耦合作用，基于实际系统的热声不稳定机理研究与控制策略研究具有较大难度，因此结构简单、热声不稳定现象显著的热声管成为该研究领域重点关注的研究对象.本文以Rijke管为例，研究热声不稳定现象的耦合控制策略.

Rijke管^［9］是一种广泛用于热声不稳定性研究的基础模型，通常为直径均匀且较长的圆柱形管道，如图1所示. 图中x和L分别表示管内轴向位置和热声管长度，声学紧凑的热源被放置在轴向位置为x_f处并将热声管分隔为热源前、后两个区域，分别用下标1和2表示.变量u、M、T、ρ和c分别表示平均流速、平均马赫数、平均温度、平均密度和平均声速，其中平均特性为沿管的常数.

依据文献［10］，热声管内流体满足质量和黏性不可压缩流体动量守恒方程，且在不考虑黏性耗散和传热效应时密度变化率满足：

\frac{\mathrm{d}\rho }{\mathrm{d}t}=\frac{\mathrm{1}}{c^{\mathrm{2}}}\frac{\mathrm{d}p}{\mathrm{d}t}-\frac{\left(\gamma -\mathrm{1}\right)}{c^{\mathrm{2}}}q（1）

式中：p、γ和t分别表示流动气体的压力、比热比和时间. 代入线性声学假设可获得系统波动方程：

\frac{\mathrm{1}}{{\overline{c}}^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}p\mathrm{\text{'}}}{\partial t^{\mathrm{2}}}-\nabla \cdot \left(\nabla \cdot p^{\text{'}}\right)+\frac{\zeta }{{\overline{c}}^{\mathrm{2}}}\frac{\partial p\mathrm{\text{'}}}{\partial t}=\frac{\gamma -\mathrm{1}}{{\overline{c}}^{\mathrm{2}}}\frac{\partial q\mathrm{\text{'}}}{\partial t}（2）

式中：\zeta表示系统通用阻尼系数；p\mathrm{\text{'}}和q\mathrm{\text{'}}分别表示压力和热量的扰动量.

热声系统中包含复杂的非线性多物理场耦合作用，如气体动力学非线性等，使热声不稳定具有高度非线性. 因此，通常非线性热声不稳定系统通过两种方式引入非线性，其一是将声学扰动参数描述为非线性函数，其二是引入非线性源项. 本文基于伽辽金展开对物理场量进行线性展开，采用具有非线性特性的改进的King’s Law热源模型作为系统的源项^［11］. 虽然热声模态是通过非线性非稳定放热反馈作用与声场扰动的耦合结果，但仅使用一个主特征频率也可近似系统的稳定性动力学^［12］，包括极限环振荡特性等. 因此，获得Rijke管单模态控制方程：

\frac{{\mathrm{d}}^{\mathrm{2}}\eta \left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}+\mathrm{2}\zeta \omega \dot{\eta }\left(t\right)+{\omega }^{\mathrm{2}}\eta \left(t\right)=

\frac{\beta }{E}\phi \left(x_{\mathrm{f}}\right)\left[{\left|\frac{\mathrm{1}}{\mathrm{3}}-{\phi }^{\text{'}}\left(x_{\mathrm{f}}\right)\eta \left(t-\tau \right)\right|}^{\mathrm{1}/\mathrm{2}}-{\left(\frac{\mathrm{1}}{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{2}}\right]（3）

式中：η（t）为时变函数；φ（x_f）为模态函数；ω为特征频率；β为表征热源特性的热源强度；E为模态函数的正交算子；τ为热源延迟时间.

1.2 Rijke管耦合系统

实验证明^［13］，利用针阀和乙烯基管耦合Rijke管系统，可以在适当的控制参数下完全抑制热声不稳定的发生. 耦合系统实验装置示意图如图2所示，其中T_c和T_h分别表示热源冷端与热端的温度，p₁和p₂分别表示振子A和振子B的压力扰动量. 利用针阀和连接管使被控振子系统A与被控振子系统B相连接. 耦合系统中主要包含时间延迟耦合作用和扩散耦合作用实现耦合管振幅的控制，其控制方程表达为：

\left\{\begin{array}{l}\frac{\mathrm{d}{\dot{\eta }}_{\mathrm{a}}}{\mathrm{d}t}+{\omega }_{\mathrm{a}}^{\mathrm{2}}{\eta }_{\mathrm{a}}+{\zeta }_{\mathrm{a}}{\dot{\eta }}_{\mathrm{a}}=\frac{{\beta }_{\mathrm{a}}}{E_{\mathrm{a}}}{\phi }_{\mathrm{a}}\left(x_{\mathrm{f}}\right)\times \\ \left[{\left|\frac{\mathrm{1}}{\mathrm{3}}-{\dot{\phi }}_{\mathrm{a}}\left(x_{\mathrm{f}}\right){\eta }_{\mathrm{a}}\left(t-\tau \right)\right|}^{\mathrm{1}/\mathrm{2}}-{\left(\frac{\mathrm{1}}{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{2}}\right]+\\ K_{\mathrm{d}}\left({\dot{\eta }}_{\mathrm{b}}-{\dot{\eta }}_{\mathrm{a}}\right)+K_{\mathrm{\tau }}\left[{\dot{\eta }}_b\left(t-{\tau }_{\mathrm{2}}\right)-{\dot{\eta }}_a\left(t\right)\right]\\ \frac{\mathrm{d}{\dot{\eta }}_{\mathrm{b}}}{\mathrm{d}t}+{\omega }_{\mathrm{b}}^{\mathrm{2}}{\eta }_{\mathrm{b}}+{\zeta }_{\mathrm{b}}{\dot{\eta }}_{\mathrm{b}}=\frac{{\beta }_{\mathrm{b}}}{E_{\mathrm{b}}}{\phi }_{\mathrm{b}}\left(x_{\mathrm{f}}\right)\times \\ \left[{\left|\frac{\mathrm{1}}{\mathrm{3}}-{\dot{\phi }}_{\mathrm{b}}\left(x_{\mathrm{f}}\right){\eta }_{\mathrm{b}}\left(t-\tau \right)\right|}^{\mathrm{1}/\mathrm{2}}-{\left(\frac{\mathrm{1}}{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{2}}\right]+\\ K_{\mathrm{d}}\left({\dot{\eta }}_{\mathrm{a}}-{\dot{\eta }}_{\mathrm{b}}\right)+K_{\mathrm{\tau }}\left[{\dot{\eta }}_a\left(t-{\tau }_{\mathrm{2}}\right)-{\dot{\eta }}_b\left(t\right)\right]\end{array}\right.（4）

式中：下标a与b分别表示系统A与系统B；τ、K_τ和K_d分别表示耦合延迟时间、时间延迟耦合强度和扩散耦合强度.

2 双管耦合控制策略预测

2.1 双管耦合系统矩阵表达

为获得式（4）所示的双管耦合系统的矩阵形式以求解耦合控制策略，首先对其进行线性化和参数代换，获得：

\left\{\begin{array}{c}{\eta }_{\mathrm{a}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+m_{\mathrm{0}}{\eta }_{\mathrm{b}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+m_{\mathrm{a}\mathrm{1}}{\eta }_{\mathrm{a}}+m_{\mathrm{a}\mathrm{2}}{\eta }_{\mathrm{a}}^{\mathrm{\text{'}}}+m_{\mathrm{3}}{\eta }_{\mathrm{b}}^{\mathrm{\text{'}}}=\mathrm{0}\\ {\eta }_{\mathrm{b}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+m_{\mathrm{0}}{\eta }_{\mathrm{a}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+m_{\mathrm{b}\mathrm{1}}{\eta }_{\mathrm{b}}+m_{\mathrm{b}\mathrm{2}}{\eta }_{\mathrm{b}}^{\mathrm{\text{'}}}+m_{\mathrm{3}}{\eta }_{\mathrm{a}}^{\mathrm{\text{'}}}=\mathrm{0}\end{array}\right.（5）

\left\{\begin{array}{l}{m}_{\mathrm{0}}={\tau }_{\mathrm{2}}{K}_{\mathrm{\tau }}\\ m_{\mathrm{a}\mathrm{1}}={\omega }_{\mathrm{a}}^{\mathrm{2}}+\frac{\sqrt[]{\mathrm{3}}{\beta }_{\mathrm{a}}}{\mathrm{2}{E}_{\mathrm{a}}}\\ m_{\mathrm{a}\mathrm{2}}={\zeta }_{\mathrm{a}}-\frac{\sqrt[]{\mathrm{3}}{\beta }_{\mathrm{a}}}{\mathrm{2}{E}_{\mathrm{a}}}{\phi }_{\mathrm{a}}^{\mathrm{\text{'}}}\left(x_{\mathrm{f}}\right){\phi }_{\mathrm{a}}\left(x_{\mathrm{f}}\right)\tau +K_{\mathrm{d}}+K_{\mathrm{\tau }}\\ m_{\mathrm{3}}=-\left(K_{\mathrm{d}}+K_{\mathrm{\tau }}\right)\end{array}\right.（6）

参数m_b1、m_b2与m_a1、m_a2表达式具有相同形式，因此，可将式（4）化为线性矩阵形式：

\frac{\mathrm{d}\mathit{Y}}{\mathrm{d}t}+\mathit{H}\mathit{Y}=\mathrm{0}（7）

\mathit{H}=\left[\begin{array}{cc}\begin{array}{cc}\mathrm{0}& -\mathrm{1}\\ \frac{-m_{\mathrm{a}\mathrm{1}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}& \frac{m_{\mathrm{0}}{m}_{\mathrm{3}}-m_{\mathrm{a}\mathrm{2}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}\end{array}& \begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \frac{m_{\mathrm{0}}{m}_{\mathrm{b}\mathrm{1}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}& \frac{m_{\mathrm{0}}{m}_{\mathrm{b}\mathrm{2}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}\end{array}\\ \begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \frac{-m_{\mathrm{b}\mathrm{1}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}& \frac{m_{\mathrm{0}}{m}_{\mathrm{3}}-m_{\mathrm{b}\mathrm{2}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}\end{array}& \begin{array}{cc}\mathrm{0}& -\mathrm{1}\\ \frac{m_{\mathrm{0}}{m}_{\mathrm{a}\mathrm{1}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}& \frac{m_{\mathrm{0}}{m}_{\mathrm{a}\mathrm{2}}}{m_{\mathrm{0}}^{\mathrm{2}}-\mathrm{1}}\end{array}\end{array}\right]（8）

\mathit{Y}={\left[{\eta }_{\mathrm{a}},{\eta }_{\mathrm{a}}^{\mathrm{\text{'}}},{\eta }_{\mathrm{b}},{\eta }_{\mathrm{b}}^{\mathrm{\text{'}}}\right]}^{\mathrm{T}}（9）

求解式（7）获得变量矩阵 Y 关于系数矩阵 H 的解：

\mathit{Y}(t)=\mathrm{e}\mathrm{x}\mathrm{p}(-\mathit{H}t)\mathit{Y}(\mathrm{0})（10）

2.2 基于Lyapunov稳定性原理的耦合参数预测

在渐近稳定性理论中，有多种用于分析动态系统微分或差分方程解的方法，其中研究平衡点附近解稳定性的Lyapunov稳定性分析方法是最重要的一类. 因此，基于该稳定性原理可预测获得耦合双管线性稳定的控制参数集，即当矩阵 H 的所有特征值具有非负实部时，系统在平衡点邻域内稳定.

耦合系统Lyapunov稳定分析图（τ=0.02）如图3所示，系统A~D分别为具有不同特征值实部最小值的系统.利用 H 矩阵特征值实部最小值（以下称最小值）来衡量系统稳定性. R_β 为热源强度比，定义为耦合系统与被耦合管系统的热源强度之比. 由图3可知，当延迟及扩散耦合强度较小时，最小值低于零，随耦合强度增大，最小值平面先升高后降低. 当R_β >1.00时［图3（a）］，随热源强度比增大，最小值平面逐渐升高，使耦合系统稳定的耦合强度参数范围逐渐增大. 当R_β <1.00时［图3（b）］，随热源强度比增大，最小值平面逐渐降低，耦合强度参数范围逐渐减小. 当R_β =1.00时，最小值平面均低于零平面，该现象说明，相同振子耦合时，扩散耦合和时间耦合无法使系统实现线性稳定，该现象归因于具有相同热源强度的振子系统声能无法扩散传递，即产生文献［14］中耦合分插图振幅死亡区域关于调谐比R_β =1.00对称分布的现象. 同时，当耦合振子热源强度足够小时，耦合振子无法激发，此时耦合振子相当于一个阻尼器，可有效衰减被耦合振子的振荡幅度，实现振幅死亡，故当R_β =0.25时，最小值平面位于零平面的上方. 由此，获得了基于Lyapunov稳定性原理的双管耦合控制策略参数.

2.3 耦合系统瞬态能量增长率

基于稳定性理论获得的参数范围准确性需进一步验证，耦合系统整体渐近稳定性能否代表单个被耦合振子稳定性的问题要进一步探究，当控制需求中考虑控制速率和控制时间时如何定量定义其评判标准，需引入系统最大能量增长率分析方法，以直观、定量地分析耦合系统或被耦合管系统的不稳定控制效果，预测耦合控制参数，为系统控制提供准确的理论依据.

单位截面积声能密度的有量纲形式为：

\mathrm{\wp }(x_{\mathrm{f}},t)={\int }_{\mathrm{0}}^L\left(\frac{{{\tilde{p}}^{\mathrm{\text{'}}}}^{\mathrm{2}}}{\mathrm{2}{\rho }_{\mathrm{0}}{c_{\mathrm{0}}}^{\mathrm{2}}}+\frac{{\rho }_{\mathrm{0}}{{\tilde{u}}^{\mathrm{\text{'}}}}^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{d}\tilde{x}（11）

将场量的伽辽金展开和能量矩阵形式进一步代入，依次求得耦合系统和被耦合管系统能量矩阵表达式为：

\left\{\begin{array}{c}{\mathrm{\wp }}_{\mathrm{总}}\left(x_{\mathrm{f}},t\right)=\frac{\mathrm{1}}{\mathrm{2}\gamma {\overline{p}}_{\mathrm{1}}}\mathit{M}\mathit{Y}{\left(t\right)}^{\mathrm{2}}\\ {\mathrm{\wp }}_{\mathrm{被耦}\mathrm{合管}}\left(x_{\mathrm{f}},t\right)=\frac{\mathrm{1}}{\mathrm{2}\gamma {\overline{p}}_{\mathrm{1}}}\mathit{N}\mathit{Y}{\left(t\right)}^{\mathrm{2}}\end{array}\right.（12）

\begin{array}{l}\mathit{M}=\left[\begin{array}{cccc}{\omega }_{\mathrm{a}}\sqrt[]{E_{\mathrm{a}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sqrt[]{E_{\mathrm{a}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\omega }_{\mathrm{b}}\sqrt[]{E_{\mathrm{b}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \sqrt[]{E_{\mathrm{b}}}\end{array}\right]\\ \mathit{N}=\left[\begin{array}{cccc}{\omega }_{\mathrm{a}}\sqrt[]{E_{\mathrm{a}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sqrt[]{E_{\mathrm{a}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\end{array}（13）

耦合和被耦合管系统最大瞬态能量增长率为：

\left\{\begin{array}{l}{g}_{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{总}}\left(x_{\mathrm{f}},t_{\mathrm{0}}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[\frac{{\mathrm{\wp }}_{\mathrm{总}}\left(x_{\mathrm{f}},t_{\mathrm{0}}\right)}{{\mathrm{\wp }}_{\mathrm{总}}\left(x_{\mathrm{f}},\mathrm{0}\right)}\right]=\\ \mathrm{e}\mathrm{x}\mathrm{p}{\left(-\mathit{M}\mathit{H}{\mathit{M}}^{-\mathrm{1}}{t}_{\mathrm{0}}\right)}^{\mathrm{2}}\\ g_{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{被耦}\mathrm{合管}}\left(x_{\mathrm{f}},t_{\mathrm{0}}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[\frac{{\mathrm{\wp }}_{\mathrm{被耦}\mathrm{合管}}\left(x_{\mathrm{f}},t_{\mathrm{0}}\right)}{{\mathrm{\wp }}_{\mathrm{被耦}\mathrm{合管}}\left(x_{\mathrm{f}},\mathrm{0}\right)}\right]=\\ \mathrm{e}\mathrm{x}\mathrm{p}{\left(-\mathit{N}\mathit{H}{\mathit{N}}^{-\mathrm{1}}{t}_{\mathrm{0}}\right)}^{\mathrm{2}}\end{array}\right.（14）

当R_β =3.00和R_β =0.25时，系统实部最小值等高线图和耦合系统最大瞬态能量增长率如图4所示. 由图4（a）可知，系统C位于零平面以下，耦合系统Lyapunov不稳定，因此在图4（c）中最大瞬态能量增长率变化曲线始终大于1，且随时间周期性上涨，系统瞬态能量始终大于初始能量，系统线性不稳定，任意初始扰值都将导致系统高幅值振荡. 而位于零平面以上的系统A、B、D、E和F，其最大瞬态能量增长率逐渐减小，且减小速率随特征值实部最小值的逐渐增大而逐渐增大.如位于最小值平面“脊”［图3（a）］处的系统A，其最大瞬态能量增长率衰减速度最快，在无量纲时间30内下降3个数量级以下；系统B衰减速度则十分缓慢.

仅通过总系统的增长率来表达单个系统的特性并不准确，尤其当具有较大热源强度的振子耦合时，被耦合管系统瞬态能量变化情况则尤为重要.因此，绘制了耦合管系统A、B、C、D、E、F中单一振子的最大瞬态能量增长率变化曲线，如图4（d）所示.由图4（d）可见，单一系统的最终能量变化趋势与耦合总系统的趋势相一致，但在耦合初期，最大瞬态能量增长率会因热源强度比大于或小于1.00而发生改变. 系统A与系统E最大瞬态能量增长率对比如图5所示.

依据最大瞬态能量增长率的定义，当其小于1时，表明瞬态能量小于系统初始能量. 由图5可见，系统A和系统E中，被耦合管系统的最大瞬态能量增长率表现出了与总系统不一致的特性. 对于系统A，当使用小于被耦合管系统热源强度的振子进行耦合时，被耦合管系统瞬态能量周期性减小；但当使用热源强度更大的系统进行耦合时，系统出现暂时性瞬态能量周期性增大，而在经过较长时间Δt_A>Δt_E后，系统能量周期性减小.

为验证双管耦合控制策略预测结果，利用式（4）进行数值求解，得到被耦合管系统扰动曲线，如图6所示. 由图6可知，被耦合管系统A、B和C在耦合开始时刻接收来自更高热源强度系统的能量而瞬间增大扰动幅值. 具有不同最小值的3个系统，有3种耦合扰动结果. 相较于系统B，具有更大最小值的系统A在瞬间增大扰动幅值后快速衰减，在无量纲时间330前获得振幅死亡，而系统B幅值下降较为缓慢，在无量纲时间500内无法获得振幅死亡. H 矩阵特征值实部最小值小于零的系统C为线性不稳定系统，扰动曲线在瞬态扰动幅值增长后逐渐形成新的极限环扰动，且由于耦合系统扰动能量的扩散，耦合后稳定极限环幅值高于原极限环幅值. 观察系统E，其采用热源强度小于被耦合管系统的耦合系统进行耦合，在图5中未观测到单管瞬态能量的增长. 图6中，系统E的扰动幅值在耦合开始时刻逐渐衰减，并快速获得振幅死亡.

综上所述，当基于耦合控制理论设计主动控制器选择合适的耦合控制参数时，耦合系统总特性与单一系统特性相一致，为更快获得稳定的系统，考虑选择特征值实部最小值更大的系统，如系统A和系统E. 系统B等临近零平面的控制策略虽能最终趋于稳定，但其实现极限环抑制耗费的时间过于漫长. 对于需要严格控制的系统，则需要选择热源强度较小的系统进行耦合控制以避免瞬态能量的周期性增长对系统造成短期的负面影响.

3 结 论

本文针对Rijke管理论模型开展双管耦合控制策略预测研究. 通过对耦合系统控制方程的线性化和矩阵化，获得了耦合系统的系数矩阵. 基于Lyapunov稳定性原理求解不同耦合控制参数下系数矩阵特征值实部最小值分布图，获得了耦合系统线性稳定区. 主要结论如下：

1）当延迟耦合强度及扩散耦合强度较小时，耦合控制策略无法实现振幅死亡；随延迟和扩散耦合强度增大，系数矩阵特征值最小值先升高后降低；线性稳定范围随热源强度比增大而增大. 当R_β =0.25时，最小值平面均位于零平面的上方，此时耦合振子相当于阻尼器，可有效实现振幅的完全控制.

2）最大瞬态能量增长率变化趋势与系数矩阵特征值最小值相关，最小值越大，最大瞬态能量衰减速率越快. 耦合总系统与被耦合管系统瞬态能量增长率变化趋势一致，但使用R_β >1.00的耦合管进行控制会造成被耦合管系统瞬态能量的上升.

3）数值求解的各系统扰动曲线表明，基于Lya- punov稳定性和最大瞬态能量增长率预测的双管耦合控制策略是有效的，可依据耦合控制策略的不同实现系统扰动的控制，包括快速获得振幅死亡、新幅值极限环振荡和瞬态幅值上升等控制结果.

参考文献

[1] 吉龙娟， 李德立， 胡光亚， 等. 预混火焰中热声振荡非线性特征的实验研究［J］. 燃烧科学与技术， 2024， 30（2）：119-127．

[2] JI L J， LI D L， HU G Y，et al. Nonlinear characteristics of thermoacoustic oscillation in premixed flame［J］. Journal of Combustion Science and Technology，2024，30（2）：119-127．（in Chinese）

[3] POINSOT T ．Prediction and control of combustion instabilities in real engines［J］．Proceedings of the Combustion Institute，2017，36（1）： 1-28．

[4] LIU X Q， ZHOU H X， LAI Y F，et al ．Equivalence ratio independence and dependence ranges of system responses for a nonlinear thermoacoustic oscillation in a Rijke tube［J］．Journal of Sound and Vibration， 2023， 547： 117545．

[5] 柯恩雷， 季晨振， 潘登， 等 ．结合系统辨识方法的旋流部分预混燃烧不稳定预测及控制研究［J］．工程热物理学报， 2023，44（8）： 2319-2329．

[6] KE E L， JI C Z， PAN D，et al ．The prediction and control of partially premixed swirl combustion instability by system identification method［J］．Journal of Engineering Thermophysics，2023， 44（8）： 2319-2329．（in Chinese）

[7] 周昊， 王恒栋， 黄燕， 等. 横向射流对Rijke型燃烧器热声不稳定控制效果的影响［J］. 中国电机工程学报， 2015， 35（12）：3075-3080．

[8] ZHOU H， WANG H D， HUANG Y，et al ．Investigation of the control effects of the transverse jet on the thermoacoustic instability in the Rijke combustor［J］．Proceedings of the CSEE，2015，35（12）：3075-3080．（in Chinese）

[9] DHARMAPUTRA B， RECKINGER P， SCHUERMANS B，et al ．BOATS：Bayesian optimization for active control of thermoacoustics［J］. Journal of Sound and Vibration，2024，582：118415．

[10] ZHAO D， LU Z L， ZHAO H，et al. A review of active control approaches in stabilizing combustion systems in aerospace industry［J］．Progress in Aerospace Sciences， 2018， 97： 35-60．

[11] HENDERSON B R， XU K G. Electric field damping of Rijke tube combustion instabilities［J］. Journal of Propulsion and Power，2018， 34（1）： 85-96．

[12] LI X Y， ZHAO D， YANG X L，et al. Unity maximum transient energy growth of heat-driven acoustic oscillations［J］. Energy Conversion and Management，2016，116：1-10．

[13] JUNIPER M P. Triggering in the horizontal Rijke tube：non-normality， transient growth and bypass transition［J］. Journal of Fluid Mechanics， 2011， 667： 272-308．

[14] ZHAO X Y， ZHAO D， CHENG L，et al. Predicting thermoacoustic stability characteristics of longitudinal combustors using different endpoint conditions with a low mach number flow［J］. Physics of Fluids， 2023， 35（9）： 094122．

[15] BONCIOLINI G， NOIRAY N. Synchronization of thermoacoustic modes in sequential combustors［J］. Journal of Engineering for Gas Turbines and Power， 2019， 141（3）： 031010．

[16] BIWA T， TOZUKA S， YAZAKI T. Amplitude death in coupled thermoacoustic oscillators［J］. Physical Review Applied， 2015， 3（3）： 034006．

[17] THOMAS N， MONDAL S， PAWAR S A，et al ．Effect of time-delay and dissipative coupling on amplitude death in coupled thermoacoustic oscillators［J］．Chaos， 2018， 28（3）： 033119．

基金资助

国家重点研发计划资助项目(2021YFB2500505)

National Key Research and Development Program of China(2021YFB2500505)