本发明属于航空航天推进系统领域，具体涉及一种燃烧稳定性裕度评估方法。
背景技术：
在火箭发动机研制阶段，通常引入外部脉冲扰动来评估发动机的燃烧稳定性。如果扰动触发了燃烧不稳定，说明系统处在双稳定区，系统稳定性裕度不足。如果扰动产生的振荡幅值随时间逐渐衰减，通过计算衰减系数，可判断系统距离分叉点的距离(理论上，系统在分叉点的衰减系数为零)，即系统稳定性裕度。但是，引入外部脉冲扰动评估系统稳定性需要额外的激励装置，同时存在触发不稳定进而破坏发动机的风险，因此需要发展无外部激励装置时量化系统稳定性裕度的方法。
现有文献可查到的方法，如lieuwen,t,onlinecombustorstabilitymarginassessmentusingdynamicpressuredata,asmej.eng.gasturbinespower,2005,127(3):478-482.lieuwen基于非线性耦合振子模型，利用稳定工作阶段压力时间序列的自相关函数计算系统衰减系数；如stadlmair,n.v.,hummel,t.,andsattelmayer,t.thermoacousticdampingratedeterminationfromcombustionnoiseusingbayesianstatistics,asmeturboexpo2017:turbomachinerytechnicalconferenceandexposition(50848),v04at04a025.stadlmair基于lieuwen方法，同时采用贝叶斯统计识别方法，可以同时识别多个声模态的衰减系数；如yi,t.,gutmark,e.j.,onlinepredictionoftheonsetofcombustioninstabilitiesbasedonthecomputationofdampingratios.j.soundvib.,2008,310(1–5):442–447，yi将lieuwen的方法转换到频域，通过拟合压力时间序列功率谱获得衰减系数。
上述方法中，lieuwen方法需要在中心频率处进行带通滤波，以消除不同声模态间的相互作用。由于噪声阶段压力时间序列为燃烧室声模态附近的宽频振荡，判断中心频率存在较大的任意性，会对计算结果产生较大影响。lieuwen通过系统发生振荡燃烧后的共振频率作为中心频率进行带通滤波，这种方法存在两个问题，首先系统参数变化时(如混合比)，共振频率会发生变化；其次在实际应用中，希望在稳定燃烧(燃烧噪声)阶段判断系统的稳定性裕度，而此时无法“精确”知道系统发生振荡燃烧后的共振频率。
理论上，系统在衰减系数为零时，说明系统会发生燃烧不稳定，即衰减系数的阈值为零。实际应用中，燃烧室在衰减系数未达到零便已发生燃烧不稳定，因此对于衰减系数的阈值需要根据已有的实验结果进行修正。
技术实现要素：
本发明的目的是解决现有方法无法精确获取系统稳定工作时固有模态频率和无稳定性判断准则的问题，提出一种燃烧稳定性裕度评估方法。
为实现以上发明目的，本发明技术方案是：
一种燃烧稳定性裕度评估方法，包括以下步骤：
步骤一、获取燃烧室压力脉动时间序列；
通过动态压力传感器获取燃烧室或燃烧室喷前的压力脉动时间序列
步骤二、计算压力时间序列功率谱；
采用周期图法计算压力时间序列功率谱，即对压力脉动时间序列
步骤三、通过参数识别获得压力时间序列功率谱的峰值频率；
利用下式拟合压力时间序列功率谱
其中，spp(ω；ωi,γ,νi)为理论压力时间序列功率谱；ω为角频率，γ为噪声的强度，ωi为峰值频率，νi为衰减系数；
步骤四、以峰值频率ωi为中心进行带通滤波；
确定滤波带宽ωh，以峰值频率ωi为中心进行带通滤波，滤波后压力时间序列记为pt；
步骤五、计算滤波后压力时间序列的自相关函数ρτ；
自相关函数ρτ按下式计算，
其中，e为数学期望计算，pt τ为pt延迟为τ的时间序列，μ为pt的平均值或数学期望，
步骤六、通过hilbert变换获取自相关函数的包络曲线ht；
通过hilbert变换获取自相关函数的包络曲线ht，计算如下：
其中，t为压力脉动时间序列采样时间，τ为自相关函数的延迟时间；
步骤七、利用指数函数去拟合包络曲线获得衰减系数；
利用指数函数
步骤八、比较衰减系数与阈值的距离，对燃烧室内燃烧过程稳定性进行评估。
进一步地，步骤五中，方差
其中，n为压力时间序列长度。
进一步地，步骤八中，阈值为0.016。
本发明技术方案与现有技术相比，具有如下优点：
本发明提出一种燃烧室燃烧稳定性裕度评估方法，利用系统识别方法，从压力时间序列频域功率谱中识别峰值频率作为带通滤波的中心频率，减少了带通滤波计算中选择中心频率时的任意性，能够更准确的计算系统固有模态频率的衰减系数；同时，本发明建立衰减系数与经验公式的换算方法，获得了不稳定发生时的衰减系数阈值，确定了稳定性判断准则，提高了燃烧稳定性裕度评估方法的工程适用性和可靠性。
附图说明
图1为本发明方法实施例中自相关函数的振荡曲线示意图。
具体实施方式
以下结合附图和具体实施例对本发明的内容作进一步详细描述。
由于现有技术中没有确定带通滤波中心频率的定量方法，基于此，本发明提出一种定量确定中心频率的方法，即提供一种燃烧稳定性裕度评估方法，采用本发明方法能够在燃烧噪声阶段评估燃烧室内燃烧稳定性的裕度，对于评定发动机的工作可靠性和评估改进方案的有效性有指导作用。该方法包括以下步骤：包括采集燃烧室压力脉动时间序列；计算压力时间序列功率谱；通过参数识别获得功率谱的峰值频率；以峰值频率为中心进行带通滤波；计算滤波后压力时间序列的自相关函数；通过hilbert变换获得自相关函数的包络；利用指数函数去拟合包络曲线获得衰减系数；基于所述衰减系数对燃烧室内燃烧过程稳定性进行评估。
本发明方法原理如下：随机宽带信号(白色噪声)激励的振荡系统(燃烧室)可以用下述方程描述，
其中，pi为第i阶模态压力振荡幅值，ωi为第i阶模态峰值频率，νi为第i阶模态衰减系数，ξ(t)为噪声激励项，求解上式压力振荡幅值的功率谱密度spp(ω)，得到下式，
其中，ω为角频率，γ为噪声ξ(t)的强度，试验中得到的压力时间序列满足式(1)，其频域功率谱形式满足式(2)，因此可以通过拟合上式获得峰值频率。
本发明提供的燃烧稳定性裕度评估方法具体包括以下步骤：
步骤一、获取燃烧室压力脉动时间序列；
通过动态压力传感器监测燃烧室或燃烧室喷前的压力脉动时间序列
步骤二、计算压力时间序列功率谱；
采用经典周期图法计算压力时间序列功率谱，即对压力脉动时间序列
步骤三、通过参数识别获得压力时间序列功率谱的峰值频率ωi；
利用下式拟合压力时间序列功率谱
其中，spp(ω；ωi,γ,νi)为理论压力时间序列功率谱；ω为角频率，γ为噪声的强度，ωi为峰值频率，νi为衰减系数；
步骤四、以峰值频率为中心进行带通滤波；
确定滤波带宽ωh，以峰值频率ωi为中心进行带通滤波，滤波后压力时间序列记为pt；
带通滤波算法成熟，以butterworth滤波器为例，一般流程为：首先根据参数要求设计出相应的模拟滤波器，然后通过脉冲响应不变法或双线性变换法将模拟滤波器变换成数字滤波器，滤波后压力时间序列记为pt；
步骤五、计算滤波后压力时间序列的自相关函数；
自相关函数ρτ按下式计算，
其中，e为数学期望计算，pt τ为pt延迟为τ的时间序列，μ为pt的平均值或数学期望，
其中，n为压力时间序列长度；
步骤六、通过hilbert变换获取自相关函数的包络曲线ht；
通过hilbert变换获取自相关函数的包络曲线ht，计算如下：
其中，t压力脉动时间序列采样时间，τ为自相关函数的延迟时间；
步骤七、利用指数函数去拟合包络曲线获得衰减系数；
利用指数函数
该步骤中，等价于用线性表达式y＝-ωiνit b去拟合曲线lnht，采用线性最小二乘方法，具体为寻找一组待定参数(νi,b)，使下述残差的平方和最小，线性最小二乘方法为现有成熟算法；
步骤8)比较衰减系数与阈值的距离，对燃烧室内燃烧过程稳定性进行评估。
衰减系数离阈值越远，说明燃烧室内燃烧状态越稳定，相应地，衰减系数离阈值越近，说明燃烧室内燃烧状态越不稳定。
文献m.l.dranovsky,v.(ed.)yang,f.(ed.)culick,andd.(ed.)talley.combustioninstabilitiesinliquidrocketengines:testinganddevelopmentpracticesinrussia.aiaaprogressinastronauticsandaeronautics,volume221,2007.定义了脉动压力衰减率δt，通过大量试验得到δt≈0.1…0.3范围内时，燃烧室可以稳定工作。通过计算，衰减系数νi与衰减率δt关系为νi＝δt/2π，因此可以确定衰减系数阈值为0.016。衰减系数离阈值越远，说明燃烧室内燃烧状态越稳定。相应地，衰减系数离阈值越近，说明燃烧室内燃烧状态越不稳定。
下面以对本发明方法做举例性的进一步详细描述。
通过动态压力传感器监测燃烧室或燃烧室喷前压力脉动时间序列
对压力脉动时间序列
通过拟合获得功率谱
确定滤波带宽ωh，以峰值频率ωi为中心进行带通滤波，滤波后压力时间序列记为pt，滤波带宽取峰值频率的10％～20％，本实例中滤波带宽ωh为峰值频率ωi的20％，即580hz；
计算滤波后压力时间序列的自相关函数ρτ，计算公式如步骤五所示，为了下一步有效识别包络曲线，自相关曲线计算长度为3～10个振荡周期(如图1所示，步骤五自相关函数的计算结果是一条振荡曲线，但其长度受延迟时间τ的控制，延迟时间τ提前给定，为了下文有效识别，通过给定合适的延迟时间τ的值，控制自相关函数的长度，使其包含3～10个振荡周期)，本实例中计算了自相关函数ρτ约8个振荡周期；
通过hilbert变换计算自相关函数的包络，计算公式如步骤六所示，利用指数函数
1.一种燃烧稳定性裕度评估方法，其特征在于，包括以下步骤：
步骤一、获取燃烧室压力脉动时间序列；
通过动态压力传感器获取燃烧室或燃烧室喷前的压力脉动时间序列
步骤二、计算压力时间序列功率谱；
采用周期图法计算压力时间序列功率谱，即对压力脉动时间序列
步骤三、通过参数识别获得压力时间序列功率谱的峰值频率；
利用下式拟合压力时间序列功率谱
其中，spp(ω；ωi,γ,νi)为理论压力时间序列功率谱；ω为角频率，γ为噪声的强度，ωi为峰值频率，νi为衰减系数；
步骤四、以峰值频率ωi为中心进行带通滤波；
确定滤波带宽ωh，以峰值频率ωi为中心进行带通滤波，滤波后压力时间序列记为pt；
步骤五、计算滤波后压力时间序列的自相关函数ρτ；
自相关函数ρτ按下式计算，
其中，e为数学期望计算，pt τ为pt延迟为τ的时间序列，μ为pt的平均值或数学期望，
步骤六、通过hilbert变换获取自相关函数的包络曲线ht；
通过hilbert变换获取自相关函数的包络曲线ht，计算如下：
其中，t为压力脉动时间序列采样时间，τ为自相关函数的延迟时间；
步骤七、利用指数函数去拟合包络曲线获得衰减系数；
利用指数函数
步骤八、比较衰减系数与阈值的距离，对燃烧室内燃烧过程稳定性进行评估。
2.根据权利要求1所述的燃烧稳定性裕度评估方法，其特征在于：步骤五中，方差
其中，n为压力时间序列长度。
3.根据权利要求1或2所述的燃烧稳定性裕度评估方法，其特征在于：步骤八中，阈值为0.016。
技术总结