钟华寿，张伟伟，肖华，叶正寅

(1.西北工业大学 航空学院，陕西 西安710072;2.中国飞行试验研究院 中航工业飞行仿真航空科技重点实验室，陕西 西安710089)

0 引言

颤振是由于弹性结构在空气动力、惯性力和弹性力的耦合作用下，产生的一种具有破坏性的气动弹性不稳定状态，整个飞行包线内都不能出现颤振。当前颤振研究主要有理论计算、风洞试验和颤振试飞三种手段。作为在飞行包线内对这些结果的最终验证手段，许多新型号的颤振飞行试验必不可少。

颤振试飞需要花费巨大的人力和物力，同时具有一定的风险，提出一种高效的、低成本、低风险的颤振预测方法对减少颤振试飞的周期和费用，降低颤振试飞风险具有重大的意义。Dimitriadis等［1］将目前主要的颤振预测方法分为两类:第一类主要是辨识结构运动方程，如NG法［2］，其采用系统辨识方法，通过两个不同速度下的响应，获得运动方程的系数，然后求解不同速度下已被辨识的系统获得颤振临界速度;第二类主要是基于稳定性判据的曲线外推，包括阻尼外推法［3］、颤振余度法［4］、包线函数法［5］以及 ARMA 方法［6］等，其中，阻尼外推法和颤振余度法在国内型号试飞中使用较多［7］。以上两类方法都无法通过一个亚临界速度点预测颤振临界特性，第二类方法甚至需要多个点，由此给颤振试飞工作带来了额外的工作量。当亚临界响应动压逐渐远离临界动压时，这两类方法的颤振预测精度都有所降低，如阻尼外推法要获得较高的颤振临界点预测精度，要求使试飞点尽量接近颤振临界点，由此给颤振试飞工作带来了一定的风险性。

本文研究之前，文献［8-9］中曾采用相似的研究方法对具有沉浮和俯仰两自由度机翼的颤振问题进行了分析和研究。本文在文献［8］的基础上进行改进和推广，进一步对具有多模态的三维弹性机翼展开分析和研究，即利用三维弹性机翼的一个亚临界响应，采用系统辨识方法获得该马赫数下的非定常气动力模型，而该模型与高度(动压)无关，接着耦合结构状态方程和气动力模型，建立闭环系统的气动弹性稳定性分析方程，该方程中的动压以参数形式出现，通过分析系统稳定性随动压的变化规律，求解系统的颤振临界特性。

1 研究方法

1.1 获得结构模态参数

通过地面模态试验或有限元计算获得模态振型、模态质量和频率等相关参数。

1.2 测试模态位移响应

设Δt为采样时间间隔，利用加速度与位移的二阶中心差分关系，解算出模态位移响应为:

1.3 模态气动力响应的求解

求解模态坐标下的结构运动方程，可获得模态气动力系数响应Fa。

式中:M为广义质量矩阵;G为广义阻尼矩阵，忽略;K为广义刚度矩阵;q为来流动压;Fb为模态空间的外激力，本文Fb=0。

1.4 非定常气动力建模

若结构振动位移较小，流场动力学特性主要表现为线性动态特征。选用离散型输入输出差分模型进行亚临界响应的气动力建模:

式中:矩阵Ai和Bi为待辨识的矩阵，以模态位移作为输入，模态气动力系数作为输出，采用最小二乘法进行估计。

将式(4)差分模型转化为状态空间模型，转化过程和各参数含义可详见文献［10-11］:

1.5 获得临界颤振特性

其中:

耦合结构状态方程式(5)和气动状态方程式(6)得到如下气动弹性系统的稳定性分析状态方程:

将气动弹性稳定性分析转化为求解状态矩阵的特征值后，通过求解不同动压下状态矩阵的特征值，绘出系统阻尼和频率随动压变化的V-g图和V-ω图，从而判读出系统的颤振临界特性。

2 算例与分析

本文采用数值方法［12］获得有关数据，对所发展的颤振预测方法的有效性进行研究和分析。

以AGARD 445.6机翼为例，取前三阶模态进行分析。通过改变同一马赫数下的来流动压(本文通过改变来流密度)，分别对马赫数0.499，0.678，0.901，0.960，1.072进行临界动压的寻找，求解时间步长为4×10－3s，颤振速度和颤振频率以无量纲化形式给出，机翼详细参数及无量纲化过程详见文献［13］。图1为颤振速度边界的试验值和计算值的对比。由图可以看出，所使用的程序和网格计算的颤振速度与试验结果基本吻合。

图1 计算值和试验值的对比Fig.1 Comparison between calculated and experimental values

在Ma=0.499，q∞=0.85q*时(q*表示对应马赫数下的颤振临界动压，下同)，对每阶模态给定一个初始扰动，模态位移响应如图2所示，求解出的模态气动力系数响应如图3所示。

图2 模态位移响应Fig.2 Modal displacement responses

图3 模态气动力系数响应Fig.3 Modal aerodynamic coefficient responses

经辨识后，耦合结构和气动状态方程，建立闭环系统的稳定性分析方程，通过特征分析获得V-g图和V-ω图，如图4所示。根据频率靠近原则，可判别出第一阶模态和第二阶模态发生耦合，第一阶模态首先失稳，图中的无量纲颤振速度为0.439 1，频率比为0.574 9。对前三阶位移响应信号施加白噪声，信噪比为20 dB，采用带通滤波，滤波前后的响应对比如图5所示。经滤波后分析得到V-g图和V-ω图，如图6所示。其无因次颤振速度为0.440 7，频率比为0.582 5，与之前的预测结果极为接近。

图4 通过特征分析获得的V-g图和V-ω图Fig.4 V-g and V-ω diagrams obtained by eigen-analysis

图5 含有20 dB噪声的前三阶模态位移滤波前后响应对比Fig.5 Comparison of modal displacement responses of the 1st three order mode with 20 dB noise before and after filtering

图6 经加噪滤波后分析得到的V-g图和V-ω图Fig.6 V-g and V-ω diagrams by using the signals with noise

重复上述方法和过程，计算Ma=0.678，0.901，0.960，1.072 对应的q∞(依次为0.875q*，0.857q*，0.85q*，0.82q*)的响应，加上 Ma=0.499，q∞=0.85q*的结果，5个马赫数点的预测结果如图7所示。由图可以看出，跨声速区出现明显的凹坑，预测效果较为理想。

图7 预测值和计算值的对比Fig.7 Comparison between predicted and calculated values

在Ma=0.678时，依次降低来流动压，分别计算0.875q*，0.750q*，0.625q*，0.500q*下的亚临界响应，以研究亚临界点逐渐远离颤振临界点时对预测精度的影响。以0.500q*为例，分析得到的V-g图和V-ω图，如图8所示。各动压比下计算出的无因次颤振临界速度和频率比如图9所示。从图9可以看出，随着亚临界点逐渐远离颤振临界点，颤振临界速度和频率比的预测精度都有所下降，但下降的值较小。从图9中还可以看出，对于本算例，在Ma=0.678时，在0.625q*～q*范围内，所预测的颤振临界速度和频率比随动压比近似呈线性变化，并且其斜率相当小，即在距颤振临界点一定范围内，利用亚临界响应点对颤振临界点预测的结果精度较高。当来流动压小于0.625q*时，预测的误差增大得较快;但当来流动压为0.500q*时，其预测的无因次颤振速度与计算值之间的误差也仅为7.62%。

图8 通过特征分析获得的V-g图和V-ω图Fig.8 V-g and V-ω diagrams obtained by eigen-analysis

图9 不同动压比下响应预测的和Fig.9 Predicted values of and at different dynamicpressure ratios

3 结论

(1)该方法只需利用三维弹性机翼一个速度较低的亚临界响应即可获得对应马赫数下的颤振临界特性，适用于两个及以上复杂模态的颤振分析，可极大降低颤振试飞或风洞试验的成本和风险;

(2)当响应测试动压与颤振临界动压的动压比降低时，颤振临界特性的预测精度有所下降，但即使动压比在0.5时，对于本算例，所预测无因次颤振临界速度的误差仍可控制在8%以内;

(3)可利用两个或多个不同动压下的亚临界响应预测和验证颤振临界特性;

(4)初步研究只是在仿真环境中进行，还需要在试验环境中，噪声强度较大的情况下施加适当的激励后验证其实用性。

［1］ Dimitriadis G，Cooper J E.Flutter prediction from flight flutter test data［J］.Journal of Aircraft，2001，38(2):355-367.

［2］ Nissim E，Gilyard G B.Method for experimental determination of flutter speed by parameter identification［R］.AIAA Paper 89-1324，1989.

［3］ Rhlin C L，Watson J J，Ricketts R H，et al.Evaluation of four subcritical response methods for on-line prediction of flutter onset in wind tunnel tests［J］.Journal of Aircraft，1983，20(10):835-840.

［4］ Zimmerman N H，Weissenburger J T.Prediction of flutter onset speed based on flight testing at subcritical speeds［J］.Journal of Aircraft，1964，1(4):190-202.

［5］ Cooper J E，Emmet P R，Wright J，et al.Envelope function:a tool for analyzing flutter data［J］.Journal of Aircraft，1993，30(5):785-790.

［6］ Matsuzaki Y，Ando Y.Estimation of flutter boundary from random responses due to turbulence at subcritical speeds［J］.Journal of Aircraft，1981，18(10):826-868.

［7］ 俱利锋，梁坤鹏，梁海洲，等.颤振边界预测技术在颤振试飞中的应用研究［J］.飞行力学，2010，28(5):79-83.

［8］ 张伟伟，于俊杰，全景阁，等.一种基于亚临界响应的颤振边界预测新方法［J］.航空工程进展，2012，3(4):390-396.

［9］ Jieun S，Taehyoun K，Seung J.Experimental determination of unsteady aerodynamic coefficients and flutter behavior of a rigid wing［J］.Journal of Fluid and Structures，2012(29):50-61.

［10］ Zhang W W，Ye Z Y.Reduced-order-model-based flutter analysis at high angle of attack［J］.Journal of Aircraft，2007，44(6):2086-2089.

［11］ Zhang W W，Ye Z Y.Effect of control surface on airfoil flutter in transonic flow ［J］.Aeta Astronautica，2010，66(7-8):999-1007.

［12］张伟伟.基于CFD技术的高效气动弹性分析方法研究［D］.西安:西北工业大学，2006.

［13］卢学成，叶正寅，张陈安.基于ANSYS/CFX耦合的机翼颤振分析［J］.计算机仿真，2010，27(9):88-91.