目前对燃烧稳定性研究表明:燃烧不稳定是火焰、流场和声场间非线性耦合的结果,而液体火箭发动机的高频燃烧不稳定大多是火焰与声场之间的热声耦合导致的^[1]。所以,燃烧不稳定低阶动力学模型通常主要考虑声场与火焰热释放之间的耦合,而将湍流等流场影响简化为随机白噪声项,则热声耦合振荡模型可用如下随机振子方程描述。

式中:p'(x,t)为时间t位置x处的脉动压力; η_i和ψ_i分别为第i阶声模态的幅值和振型; Dε(t)为燃烧室压力振荡中湍流引起的白噪声扰动,D为噪声强度,ε(t)为单位高斯白噪声; f为火焰热释放对压力振荡幅值的反馈函数; α_i为声模态衰减系数; ω为角频率。

当处于稳定燃烧阶段时,假设火焰热释放反馈函数为一阶导数线性项f=2βη^·,记v=α-β为线性衰减系数,则系统模型方程可分解为单个特征模态的幅值方程,即

式(2)中,当v<0时系统处于线性稳定状态。由此可知,根据辨识得出的衰减系数v与0之间的距离即可评定燃烧装置的燃烧稳定性裕度。

衰减系数辨识方法中最重要的就是特征模态的获取,而前期研究表明:带通滤波参数选取对辨识精度有较大影响,特别是当处于稳定燃烧阶段,无法精准获取燃烧室内热声耦合频率,进而影响滤波宽度的选择^[13]。现有辨识方法采用的是将燃烧室固有声学频率作为特征模态中心频率,进而测算合适的滤波宽度,以此来获取特征模态信息。

然而燃烧噪声阶段,燃烧室压力监测信号为宽频振荡,较难选择中心频率。另外带通滤波宽度也对衰减系数辨识有严重影响。当带通滤波宽度较小时,带通滤波提取的模态信息会丢失部分热声耦合模态从而导致辨识精度下降; 而带通滤波宽度过大时,带通滤波提取的模态信息会包含多个特征模态信息从而导致辨识精度下降。所以得到准确可靠的热声耦合模态信息不仅要求包含确定的中心频率,还要求提取模态信息在频谱上存在一定宽度,以包含所有与该模态耦合的湍流噪声信息。

随着信号处理方法的发展,经验模态分解^[15](empirical mode decomposition,EMD)、变分模态分解^[16](variational mode decomposition,VMD)、局部均值分解^[17](local mean decomposition,LMD)等方法逐渐应用于时间序列分解中,然而EMD、LMD方法都存在模态混叠、边界效应和噪声鲁棒性弱等不足,而VMD方法虽然具有良好的噪声鲁棒性,但由于其需要预先设定好分解层数,致使分解结果可能存在过分解或欠分解情况。本文采用文献[18]提出的基于峭度的VMD分解方法计算分解模态中的最大相关分量的峭度值,当峭度值达到最大时停止分解,最终得到最优模态分解结果。

变分模态分解方法将原始信号分解成n个限制带宽的本征模态函数(intrinsic mode function,IMF),表示为

u_n(t)=A_n(t)cos[ω_n(t)] (3)

式中:u_n(t)为分解得到的n个IMF; A_n(t)为u_n(t)的瞬时幅值; ω_n(t)为u_n(t)的瞬时频率。每个分量在中心频率处集中,可用高斯平滑调制信号来估算带宽,由于VMD分解的稀疏性,分解时可将其转换为

式中:(σ(t)+(j)/(πt))u_n(t)为原始信号基于希尔伯特变换的解析; e^-jω_nt为模态函数的频谱调制修正。

为解决上述的约束最优化问题,将约束变分问题转换为非约束变分问题,添加二次惩罚项和拉格朗日因子,得到如下增广拉格朗日函数以求上述问题的最优解。

式中:α为乘法因子; λ为拉格朗日因子。

利用变分模态分解方法将原始信号分解为K个IMF,过程如下。

1)初始化{u_k}、{ω_k}、λ₁为0。

2)对所有ω≥0,更新泛函u^{^}_k,即

3)对所有ω≥0,更新泛函ω_k,即

4)对所有ω≥0,更新拉格朗日因子λ,即

5)重复循环步骤2)～4),直到满足如下约束条件

式中:τ为噪声容限; K为预设分解层数。

从分解过程可以看出分解层数与噪声容限对分解结果影响最大。基于峭度的分解层数K确定方法是以最大相关IMF的峭度作为对照依据,以峭度最大作为循环截止条件来确定K的最优解。具体步骤如下。

1)预设分解层数K^{^}_n,对原始信号进行变分模态分解,得到。

2)计算u_i(t)与原始信号x(t)的相关系数,得到最大相关分量u_j(t)。

3)计算u_i(t)的峭度值α_n=∑u⁴_j(t)。

4)重复步骤1)～3),直到找到峭度最大时的分解层数K。

在确定分解层数后,利用残差指数REI确定噪声容限,在τ∈[0,1]范围下选择残差最小的噪声容限τ。残差指数计算式为

对式(2)进行傅里叶变换并取模的平方后,可得到幅值方程的功率谱函数为

S_η(ω)=D/((ω²_i-ω²)²+4v²ω²) (11)

式中:ω_i为特征模态的角频率。

系统参数D、v可由压力时间序列辨识得到,即频域衰减系数辨识方法。

根据维纳—辛钦定理可得时域计算的自相关函数与频域计算的功率谱S_η(ω)的关系为

则可得到幅值方程的自相关函数为

由于线性衰减系数v远小于角频率ω_i,则自相关函数第二项可忽略,由此可得

由式(14)可知自相关函数幅值在频率ω_i附近振荡,且幅值成指数衰减,衰减率为线性衰减系数v,由此可得线性衰减系数,即时域的衰减系数辨识方法。

理论验证信号采用非线性耦合振子模型,模型方程为式(2),设定模型参数如表1所示,需辨识的信号由3种不同衰减率模拟信号叠加而成,数值模拟时间为10 s,采样频率为10 kHz。

表1 数值模拟信号系统参数
Tab.1 System parameters of numerical analog signal

图1为3种指定参数模拟信号叠加后的幅值与频谱图。

图1 数值模拟信号及其频谱图
Fig.1 Numerical analog signal and spectrograms

图2为带通滤波中心频率选择影响,当带通滤波宽度选取为中心频率的20%时,改变带通滤波中心频率,可见随着中心频率选取从1.0～1.2 kHz,频域方法辨识结果越来越向模态3靠近; 时域方法辨识结果先减小,然后趋于模态3的预设衰减系数值(-0.05),误差越来越大,逐渐偏离相对应的模拟信号预设衰减系数。其中,时域方法先降低是因为2个相近模态(0.9 kHz、1.0 kHz)之间的混叠效应。

图2 中心频率对系统辨识结果影响
Fig.2 Effect of domain frequency on attenuation coefficient identification

图3为带通滤波宽度对模拟信号的辨识结果影响,可见时域方法辨识结果先增大后减小,最后仍存在一定误差,频域方法辨识结果一直比较稳定。这是因为时域方法受滤波宽度影响较大,起初滤波宽度较小,模态2(1.0 kHz)只有部分信息进行辨识; 进而因模态1(0.9 kHz)与模态2(1.0 kHz)两者中心频率相近,滤波宽度增大会将模态1的影响引入辨识结果中; 最后随着滤波宽度继续加大,模态3(1.2 kHz)也影响着辨识结果。频域方法辨识结果较好是因为中心频率不变,滤波宽度影响程度较小。

在传统模态提取方法中,频域方法对中心频率的选取比较敏感,时域方法则需要合适的中心频率与滤波宽度。然而在燃烧噪声阶段,燃烧室压力表现为以固有频率为中心的宽频振荡,燃烧噪声自身的随机性会对辨识结果有较大影响。因此,准确可靠的模态提取方法不仅需要准确的中心频率,也需要一定的频带宽度以包含所有模态信息。

图3 滤波宽度对系统辨识结果影响
Fig.3 Effect of filter width on attenuation coefficient identification

图4为本文VMD算法处理的结果。处理VMD分解得到的结果如表2所示,可见VMD算法可以高效分离得到中心频率与模态信号,所得的模态经过时/频域方法辨识后与带通滤波结果误差较小,可以使用基于峭度最大的变分模态分解算法处理压力信号,与带通滤波方法相比,变分模态分解可以自动计算得到中心频率及信号,无需先验的热声谐振频率。

表2 VMD与带通滤波方法辨识结果
Tab.2 Identification results of VMD and bandpass filtering method

图4 VMD分离模态结果
Fig.4 VMD separation results

试验发动机采用气氧/乙醇作为推进剂,点火过程利用火花塞来启动发动机,喷嘴类型为针栓式,试验工况设置为:氧化剂(O₂,300 K)为112 g/s,燃料(质量分数为50%的乙醇水溶液,300 K)为218 g/s,氧燃比为0.51; 该工况下燃烧室内总温1 506 K,声速为910 m/s,模型燃烧室为圆筒形结构,燃烧室直径为70 mm,长度为350 mm,收缩段长度为20 mm,通过计算可得1阶纵向声学固有频率为1.264 kHz,2阶 为2.528 kHz,4阶为5.056 kHz。压力传感器选用一个采样频率为125 kHz的Kistler 6043A动压传感器,脉动压力测点距喷注面板平面15 mm。

图5为燃烧室脉动压力监测信号图,从6.1 s开始进行点火,在6.2～7.2 s之间为燃烧阶段,7.2 s时刻停止供应氧气,燃烧持续时间约1.0 s,期间燃烧室内压力约为1.0 MPa,由放大图可知,发动机处于稳定燃烧状态。为保证实验安全,发动机点火在富氧条件下进行,熄灭在富燃条件进行。

图6为燃烧室压力功率谱密度,可见在稳定燃烧阶段(6.2～7.2 s)存在2个凸频,在频率2.534 kHz、4.814 kHz附近显示宽频燃烧振荡,分别对应理论计算出的燃烧室2阶纵向、4阶纵向声模态。

图5 燃烧室动态压力监测图
Fig.5 Dynamic pressure of combustion chamber

图6 燃烧室压力信号功率谱密度
Fig.6 Power spectral density of combustion chamber pressure

图7为VMD算法对脉动压力数据的处理结果,经过计算:分离模态数K取值为9,噪声容限τ为0.92,共分离出9个本征模态函数,中心频率分别为0.195 kHz、2.510 kHz、3.710 kHz、4.882 kHz、7.421 kHz、10.156 kHz、20.210 kHz、22.613 kHz、32.738 kHz,其中模态2、模态4的中心频率(2.510 kHz、4.882 kHz)与燃烧室固有2阶纵向声频率(2.528 kHz)和4阶纵向声频率(5.056 kHz)相近,且与压力信号功率谱密度图中2个宽频凸峰相近(2.534 kHz、4.814 kHz),可以看出VMD算法可准确地分离得到谐振模态中心频率。

进而对IMF2和IMF4进行衰减系数辨识,辨识结果如表3所示。

图8为滤波宽度对线性衰减系数影响,随着滤波宽度增加,衰减系数随之减小。所以在带通滤波宽度为中心频率的30%时,衰减系数变化基本稳定,因此时域、频域方法中带通滤波宽度均选择为中心频率的30%。

图7 VMD分离模态结果
Fig.7 VMD separation results

表3 VMD辨识结果
Tab.3 Identification results of VMD method

图8 滤波宽度对系统辨识结果影响
Fig.8 Effect of filter width on attenuation coefficient identification

进而对滤波结果进行衰减系数辨识,所得结果与VMD方法所得结果对比如表4所示,对比结果可得无论2阶还是4阶模态的衰减系数,时域方法所得结果相差较小,频域方法所得结果差异相对较大,具体表现在基于VMD方法滤波后的频域结果更高,主要是因为VMD滤波所得的频带宽度更大。

表4 VMD与带通滤波方法辨识结果
Tab.4 Identification results of VMD and bandpass filtering method

图9为衰减系数随时间变化图,对发动机稳定燃烧阶段数据进行滑窗切片处理(滑窗宽度为0.20 s,前进步长为0.01 s),对切片数据进行衰减系数辨识,由此可见衰减系数在稳定燃烧阶段保持稳定,2阶纵向模态在开机表现出较强的不稳定趋势,4阶纵向模态在启停期间并未有明显不稳定趋势。出现这种现象的原因是在启动点火阶段,燃烧室燃料分布不均匀,还未形成稳定燃烧,所以未能形成热声耦合共振场^[19-20],得到的衰减系数本质上反映了燃烧室结构的频率选择性。

图9 衰减系数随时间变化
Fig.9 Variation of attenuation coefficients over time