基于时序模型的颤振预测技术研究
2018-09-10刘立坤张春蔚王东森
刘立坤 张春蔚 王东森
摘要：本文主要研究了颤振边界预测应用技术，针对湍流激励下的结构响应数据，采用遗忘因子最小二乘法建立时序模型，结合Jury稳定性判据构造不同飞行速度下的稳定性参数，通过外推得到颤振临界速度预测结果。以仿真试验为例，讨论了临近边界速度大小和响应数据时间长度等因素对预测精度的影响。最后结合颤振风洞试验实测数据验证了该方法的合理性和工程实用性。
关键词：颤振边界预测;时序模型;遗忘因子最小二乘法
中图分类号：V217 文献标识码：A
颤振飞行试验是新型或结构有重大改变的飞行器都必须进行的试飞项目之一，它使用真实飞行器在实际飞行条件下进行试验，具有典型的高风险、高耗费和长周期的特点。颤振边界预测技术是颤振试飞中确定颤振边界的主要手段。由于颤振飞行试验条件的复杂性，试验数据的品质较低，颤振数学模型存在诸多不确定因素，颤振包线一般很难直接从试验获得，因此如何基于亚临界响应测量与分析进而完成对包线的预估十分重要[1]。
国内在风洞颤振试验中，常用的亚临界颤振边界预测方法主要包括阻尼外推法和稳定性参数法，但多数都需要通过模态参数辨识，获取不同条件下结构的固有频率和阻尼比。阻尼外推法通过直接外推临界参数确定颤振边界;稳定性参数法则由模态频率和阻尼换算得到特征多项式的复数特征根，通过根和系数的关系构建劳斯判据或Jury判据稳定性参数，然后进行临界参数外推确定颤振边界。稳定性参数法在国内外已有广泛的应用，常见的著作和论文主要有管德院士的《气动弹性试验》[2]、宇航出版社出版的《防空导弹结构与强度》[3]、中国空气动力研究与发展中心郭洪涛等编写的《颤振试验数据处理技术研究》[4]等。然而，由于存在激励力不足、测量点有限、试验条件复杂及测试条件恶劣等问题，颤振试验特别是飞行试验测量的试验数据品质较低，数据往往具有模态密集、有效样本短和数据信噪比低等特征，这些问题的存在大大增加了模态参数识别的技术难度，进而影响到颤振边界预测结果的准确性。
目前，国内基于时序模型的颤振边界预测方法较为少见，本文通过直接识别时序模型多项式系数获得稳定性参数，不必经过模态参数辨识，提供了一种基于时序模型的颤振边界预测方法，以仿真试验为例，讨论临近颤振边界不同速度大小和响应数据时间长度等因素对预测精度的影响。最后结合颤振风洞试验实测数据验证了该方法的合理性和工程实用性。
1 方法介绍
1.1 遗忘因子最小二乘算法
遗忘因子最小二乘算法最初由Lennart Ljung于1999年提出[5]，广泛应用于控制系统辨识领域，其思想是利用前一时刻的响应对下一时刻进行预测，根据当前的预测误差设置适当的增益值，来影响预测的趋势，经过迭代运算，使得最终的预测误差达到最小，并以误差最小时的参数作为系统辨识的结果。具体算法如下，假设时间序列{y（t）;t=1，2，…，N}是测量到的一组响应值，则有：式中：ε（t）是t时刻预测误差Ψ班（t）是t时刻斜率，θ（t）是t时刻的参数估计值，y（t）是t时刻的响应值（t）是由t-1时刻响应得出的t时刻响应预测值，增益K（t）为当前预测误差对参数估计的影响程度，λ是遗忘因子取值为1，通过不同时刻的循环迭代使得J（θ）最小，并以θ值作为系统辨识的结果[6～8]。
1.2 ARMA气动弹性系统辨识[9～12]
假设{y（t）;t=1，2…，N}为N自由度气动弹性系统经大气湍流激励的响应，将其表示成自回归滑动平均（ARMA）模型的差分方程形式：式中：e（t）是零均值、方差桅σ2的白噪声;变量z-1为后移算子，z-1y（t）=y（t-1），多项式系数A（z-1）和B（z-1）分别为：式中：n和m分别为自动回归（AR）和滑动平均（MA）多项式阶次。
由于白噪声项e（t）无法测量，由过去时刻的数据y（t）得到的预测旬（t）为：
定义预测误差为：
联立式（8）和式（11）得：
将式（13）带入遗传因子最小二乘算法，经过迭代运算，即可得到ARMA多项式系数。
1.3 稳定裕度估计
假设经过遗忘因子算法辨识出ARMA多项式系数，将其表示成特征多项式形式：
对于离散系统使用Jury稳定性判据判断系统的稳定性，对于k=1，3，…，n-1，定义k阶方阵Xk、Yk为：Jury判据要求：
Jury判据要求：
计算不同速度、不同模态结构的稳定性因子F-（n-1），n为ARMA模型阶次，绘制稳定性因子随速度的变化曲线，并使用二次曲线拟合的方法进行外推，就能得到稳定性参数为零时对应的速度值，即为颤振临界速度值。
2 仿真数据验证分析
为验证该方法在湍流激励情况下的有效性，以得克萨斯A&M大学航空航天工程系的俯仰一滚转气动弹性风洞模型为算例[3]，主要研究了不同速度区间下的数据点数和响应数据时间长度对预测结果的影响。
在MATLAB中进行数值仿真试验，其动力学模型见式（18），具体取值见参考文献[6]。該模型为二自由度系统，模型阶次真值为4，颤振速度Vflutter为12.7m/s，对应的速度频率和速度阻尼曲线如图1和图2所示。
采用零均值高斯白噪声模拟真实大气湍流激励形式进行激励，首先分析了响应信号时长为60s、无噪声情况下不同速度区间数据的预测结果，见表1。对于新型或重大改型飞机，速度在0.8Vd～0.98Vd之间须经过颤振飞行试验验证，（Vd为最大设计俯冲速度，根据CCAR25部相关条款Vd≤Vflutter/1.15），表1给出了8种典型速度区间的预测结果。其中，fz-fit为时序模型预测结果，g-fit为阻尼外推法预测结果。图3和图4给出了两种情况的颤振速度预测结果图[7-9]。
分析結果表明，低速试验点时，阻尼较为分散，阻尼外推法预测误差较大，时序模型预测结果比阻尼外推法精度高且较为保守。随着速度逐渐接近临界颤振速度，两种方法预测结果趋于一致并且都较为准确。
在颤振飞行试验中，为了得到有足够置信度的信息，大气湍流激励方法常要求飞行员持续保持飞机稳定平飞状态一定时间，采集足够长的结构响应信号，这大大增加了试验的成本和风险。选取速度区间0.8Vd～0.98Vd，针对不同时间长度的响应信号数据进行了分析，表2给出了具体分析结果。
分析结果表明，随着响应信号数据样本的减少，基于时序模型的预测方法预测精度逐渐降低。但即使是5s的短时数据样本情况，该方法仍然具有一定的精度。
3 颤振飞行试验数据处理中的应用
某型飞机风洞试验，采用大气湍流激励，利用布置在平尾、垂尾等位置的加速度传感器测量结构响应。在速度201～381m/s的部分试验点进行试验，图5～图7给出了部分速度下左平尾结构响应自功率谱曲线，在速度381m/s时，频域两模态频率归一、响应幅值明显放大，速度381m/s已非常接近颤振速度。考虑到试验的安全性，在该速度点终止试验。分别使用基于时序模型的预测方法和阻尼外推法进行颤振速度预测[10～12]。
图8给出了部分试验点的颤振速度预测结果，fz为稳定性因子值，g为结构阻尼比，fz-fit为稳定性因子外推拟合曲线，g-fit为阻尼外推拟合曲线。图中，阻尼外推法颤振速度预测值约为411m/s，基于时序模型的颤振预测方法预测值约为403m/s，两种方法预测结果一致性较好。
4 结论
采用基于时序模型的颤振边界预测方法，对湍流激励下的气动弹性模型进行颤振边界预测，并与传统阻尼外推法进行了对比分析。仿真试验结果表明，基于时序模型的颤振预测方法在接近颤振临界速度时与阻尼外推法预测结果一致，并且由于稳定性因子随速度成近似线性变化的关系，在低速试验点时其预测精度较阻尼外推法有一定的优势。
将该方法应用于某型飞机风洞试验数据，对平尾进行颤振速度预测，结果表明，该方法能够得出工程适用的颤振速度预测值。此外，与阻尼外推法相比，该方法对于湍流激励得到的短时响应信号也有一定的适应能力，在实际飞行试验中，能够有效保障科研试飞的安全，降低飞行试验的风险，具有一定的工程使用价值。
参考文献
[1]张伟伟，钟华寿，肖华，等.颤振飞行试验的边界预测方法回顾与展望[J].航空学报，2015，36（5）：1367-1384.
[2]管德.气动弹性试验[M].北京：北京航空航天大学出版社，1986.
[3]许椿荫.防空导弹结构与强度[M].北京：中国宇航出版社，2006.
[4]郭洪涛.颤振试验数据处理技术研究[D].绵阳：中国空气动力研究与发展中心，2005.
[5]Jung L.System identification：Theory for the user[M].Beijing：Tsinghua University Press，2002.
[6]Baldelli D H，Zeng J，Lind R，et al.Flutter-prediction toolfor flight-test-based aeroelastic parameter-varying models[J].Journal of Guidance Control&Dynamics，2009，32（1）：158-171.
[7]Matsuzaki Y Flutter boundary prediction based on Jury，s stabil-ity criterion：a review[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and Materials Conference，2011.
[8]Matsuzaki Y.Flutter boundary prediction of digitalized smartmultimode systems using steady-state responses[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and MaterialsConference，2006.
[9]Torii H，Matsuzaki Y Flutter margin evaluation for three-modediscrete-time systems[J].Journal ofAircraft，2013，38（1）：42-47.
[10]AlAA.Flutter prediction of multimode systems from their tur-bulent responses by Jurys criterion[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and Materials Conference，2013.
[11]Matsuzakir Y Flutter boundary prediction of an adaptive smartwing during process of adaptation[R].Proc Spie，2005：5764.
[12]Matsuzaki Y，Torii H.Flutter boundary prediction of a smartwing in process of structural adaptation[C]//Aiaa/asme/asce/ahs/asc Structures，Structural Dynamics and Materials Conference，2007.
