王洪波 赵长见 廖选平 王亮

(1.国防科学技术大学航天科学与工程学院, 长沙 410073) (2.中国运载火箭技术研究院, 北京 100076)

基于飞行工作模态分析的飞行器动载荷识别研究*

王洪波1,2赵长见1,2廖选平1,2王亮2†

(1.国防科学技术大学航天科学与工程学院, 长沙 410073) (2.中国运载火箭技术研究院, 北京 100076)

基于飞行遥测振动数据,本文提出了基于工作模态分析的飞行器动载荷识别方法.首先,详细介绍了ERA环境激励模态辨识方法的理论.其次,给出了飞行器结构动力学建模方法.再次,提出了基于工作模态辨识的飞行器动载荷辨识计算工作流程,详细分析了其中的注意点.最后,通过算例验证了方法的可行性,其中基于飞行器飞行振动遥测数据,采用环境激励模态辨识方法辨识其各时刻的模态,包括模态频率和模态振型,再利用振动响应的模态叠加原理和模态正交理论,获取各时刻飞行器低阶模态的响应,再结合模态剪力和模态弯矩进行动载荷识别.

模态辨识, ERA, 工作模态, 动载荷, 振动, 载荷识别

引言

高超声速飞行器在空中飞行过程中,不可避免地会受到大气扰动的影响,由于该扰动使得高超声速飞行器出现低频整体弹性响应,由此产生的飞行器上的分布载荷称为动载荷,该载荷与扰动激起的模态数目、阶次以及该阶模态响应有关,一般是由低阶模态影响,且从低阶到高阶影响逐渐减弱,因此对于模态频率较低的飞行器,其动载荷在整体载荷中所占成分较大,必须考虑.基于此,基于各类实测数据的飞行器载荷辨识成为了重要研究内容.

动载荷辨识作为结构动力学的一个逆问题,是根据结构动力学特性和实测响应估计结构所受的动态载荷,学者们在该方面进行了大量的研究.周盼和张权[1]对动载荷识别时域方法的研究现状与发展趋势进行了综述,分析了各自的优缺点.Kazemi[2]提出了利用应变的时间积分识别动载荷,研究结果表明利用应变第三时间积分求解的载荷精度很高,且无需应用正则化,其对测试噪声和传感器位置均不敏感,识别结果稳定.程良彦等人[3]采用小波反卷积方法对拱结构的冲击激励辨识进行实验研究.识别载荷与实测载荷整体吻合很好,但由于小波的紧支性构造的数据较离散,造成冲击激励的上升沿、下降沿识别结果差.常晓通和闫云聚[4]为避免传统载荷识别过程中产生的矩阵求逆病态和对初值敏感及累计误差等问题,将遗传算法应用到载荷识别过程中,将此动力学的反问题转化为结构动力学正计算,并且利用遗传算法进行寻优得到最优参数,从而得到待识别载荷时间历程的估计.通过仿真计算,所识别载荷计算振动响应与测量响应的误差为7%.杨帆和张方[5]针对工程结构中的薄板系统,基于正交小波级数分解理论,构造正交小波基函数,将待识别的分布动态载荷在正交小波基空间中进行小波级数分解,从模态理论出发给出了小波级数系数与结构测点响应信息的线性关系,计算出小波系数,通过小波重构,还原薄板结构表面承受的分布动载荷.本文的载荷识别理论适用于各种不同载荷类型.

另外,对于高超声速飞行器的低频工作模态辨识,需要采用环境激励模态辨识技术.该方法不需要测量激励型号,而仅依靠各通道的时域响应数据进行系统的模态辨识.环境激励模态辨识技术的程序是：首先,进行数据采样；然后对采样数据进行自相关和互相关计算,在进行多个测点的模态参数识别处理中,需要选取某个测点做参考点.一般情况下,选取响应较小的测点做参考点,计算其它测点与该参考点的互相关函数.然后,将计算出来的互相关函数作为输入数据,利用如ITD法、STD法、复指数法、ARMA模型时序法以及ERA法等时模态参数辨识方法进行参数识别.目前,该方法已应用于桥梁、高层建筑、汽轮机、飞机和汽车等的模态参数识别[6-11].

综上所述,对于高超声速飞行器动载荷设计,缺少飞行过程中的动载荷辨识,因此无法直接获得其对飞行器飞行载荷的贡献比例,也无法验证设计时的计算结果；另外,对于高超声速飞行器,其飞行过程中的动载荷测量较为复杂,往往加速度响应较容易获取,因此如何利用常规的加速度测量结果进行动载荷辨识需要进行研究.基于此,本文提出了基于工作模态分析的飞行器动载荷识别方法,其中利用飞行器飞行振动遥测数据,采用环境激励模态辨识方法辨识其各时刻的模态,包括模态频率和模态振型,再利用振动响应的模态叠加原理和模态正交理论,获取各时刻飞行器低阶模态的响应,再结合模态剪力和模态弯矩进行动载荷识别反算.

1 模态辨识技术

本文采用ERA法进行飞行器工作模态辨识.ERA法属于一种多输入多输出的时域整体模态参数辨识方法.原理是利用实测的脉冲响应数据或自由响应数据来构造Hankel矩阵,采用奇异值分解的方法,求得系统的特征值与特征向量,从而求得模态参数.该方法由于使用了现代控制理论中的最小实现原理,使得计算量大大减小,有很好的精度,是目前最完善、最先进的模态参数辨识方法之一.

对于n维线性系统,当振动系统响应采用加速度传感器测量时,振动方程用向量可表示为:

(1)

设离散时间点为k=0,1,2,…,采样时间间隔为Δt,则t=t0+kΔt,对公式(1)Z变换并整理后有:

(2)

构造Hankel矩阵:

(3)

整理得:

(4)

令k=1,对H(0)做奇异值分解H(0)=U∑VT,可推导出:

∑1/2VTEL

(5)

又A1=∑-1/2UTH(1)∑-1/2；B1=∑1/2VTEL；

设系统矩阵A的特征值矩阵为Λ,特征矢量矩阵为ψ′,由指数矩阵的性质,从而知A1的特征矢量与A的相同,A1的特征值矩阵为:

Z=eΛΔt=diag(z1,z2,…,z2n)

(6)

式中Z的对角矩阵元素为zi=eλiΔt,i=1,2,…,2n；A的特征值矩阵Λ=diag(λ1,λ2,…,λ2n),且

(7)

由此可确定各个模态振动的固有频率、阻尼比和模态矩阵:

模态矩阵 Φ=Gψ

2 飞行器结构动力学模型

飞行器采用梁质量块动力学模型,如图1所示,各分站质量以带质量的0维单元连接在各节点上,模型各舱段的刚度由该舱段的材料以及等效厚度的形式给出,梁单元采用变截面铁木辛克梁有限元模型,由于飞行器为连续气动外形,因此各节点处截面半径不相等,因此采用如公式(8)的型函数导数与材料参数矩阵乘积在单元长度上积分的方法,得到各单元的刚度矩阵.

图1 飞行器结构动力学模型Fig. 1 Dynamic model of aircraft

(8)

梁单元的质量矩阵采用集中质量单元,将各分站质量以及转动惯量填入对应节点位置,其中对应节点的转动惯量计算公式如(9)所示.

(9)

其中,me为该站上的集中质量；R为该站处的弹体半径；h为处于端点时取它与相邻点间距之半；处于中间点时取它相邻的前后两点间距之半.

基于以上给出的单元刚度矩阵和质量矩阵,进行总刚度矩阵和质量矩阵的组装,从而得到飞行器的结构动力学模型.

3 飞行器动载荷辨识计算流程

图2给出了基于某一段遥测数据的飞行器动载荷辨识的计算流程.首先是基于振动遥测数据的飞行工作模态辨识,其中对多个通道的飞行振动遥测数据进行预处理,使用带通滤波筛选出低频模态所在频带的信号,再对信号进行重采样,减少信号长度,通过环境激励下模态辨识方法进行模态辨识；其次是由于飞行过程中的飞行器分布质量和模态质量是无法辨识的,因此根据理论计算模型得出对于时刻的质量分布以及模态质量计算值,使用辨识后的模态频率对低阶模态剪力和模态弯矩进行计算；再次,根据辨识的模态利用模态正交性分离出低阶模态的响应；最后结合模态弯矩和模态剪力计算结果得到飞行器动载荷的反算结果.

图2 计算流程Fig.2 Flow chart of dynamic load identification

计算过程中有几个处理注意点：

a)模态剪力和模态弯矩的计算方法使用前提

由于计算过程中是混合了基于遥测数据辨识和理论计算模型的,尤其是模态剪力和模态弯矩的计算方法,因为天地差异的存在,模态频率辨识值一般与计算值存在一定的差异,本文计算模态剪力和模态弯矩时模态振型、分布质量和模态质量为计算值,模态频率为辨识值,这样可以提高计算准确度.综上所述,辨识模态与计算模态间的差异需要进行分析,尤其是模态振型的差异,因此在使用该计算方法前,需要对两模态振型的正交性和一致性进行检查和分析.

b)低阶模态位移计算方法

低阶模态坐标下的响应计算是基于多个通道的遥测数据,使用模态正交性计算出模态坐标下的响应,但由于动载荷是由模态剪力和模态弯矩与模态位移的乘积计算出来的,而飞行振动遥测数据一般是使用加速度传感器测得,因此采用模态正交性计算出的模态响应是模态加速度的响应,因此需要对其进行处理从而得到模态位移值,本文采用了对短时间的模态加速度的响应进行频域分析,对分析结果转化为模态位移的频域分析结果,以此得到模态位移的均方根值,再假设其为正态分布,对均方根值放大3倍得到99%概率的模态位移响应值.

4 算例

4.1 飞行器工作模态辨识

图3给出了针对某一段遥测数据的模态参数辨识的计算流程.首先对数据进行预处理,再使用带通滤波筛选出待辨识频带的信号,再对信号进行重采样,减少信号长度,最后通过模态辨识方法进行模态辨识.

图3 工作模态辨识计算流程Fig. 3 Flow chart of operational mode analysis

图4 处理前后参考点信号的功率谱密度曲线对比Fig.4 Comparison of power spectral density of the signal before and after filiter

图5 时频分析结果Fig.5 Result of time-frequency analysis

选择一段1.5s长度的信号,对其进行预处理,由于原信号采样频率为5120Hz,而飞行器动载荷主要由低阶模态构成,因此低通滤波器频带选择为200Hz进行滤波.

对信号处理后的功率谱密度分析结果和时频分析结果分别如图4和图5所示,从图上可以发现低频的谐振峰位置在45Hz和88Hz左右.

原始信号共4个通道,分别对其使用以上预处理方法进行处理后,采用ERA模态辨识方法对飞行器的工作模态进行辨识,模态频率和振型的辨识结果如表1所示.对比可以发现,辨识的低阶模态频率与计算值相差最大2.25%,从而互相验证了辨识结果和计算模型的准确性.

表1 模态频率辨识结果Table 1 Result of operational modal frequency analysis

4.2 飞行器低阶模态坐标下的响应计算

根据以上的飞行器工作模态辨识的结果,利用模态振型的正交性,求解前两阶模态的模态加速度响应.在此之前,首先对比辨识模态和理论计算模态,典型模态对比结果如图6所示.从图上可以发现,辨识结果与理论计算结果基本一致,吻合较好.据此计算出前两阶模态的模态加速度响应曲线如图7所示,根据模态加速度计算结果,对其进行频域分析后,得到一阶模态位移均方根值,再乘以3得到一阶模态位移99%概率估计值如表2所示.由于二阶模态位移响应过小,因此后续反算飞行器动载荷时仅计算一阶模态动载荷.

图6 辨识模态和理论计算模态对比Fig. 6 Comparison of the results from operational modal shape analysis and theoretical calculation

图7 前两阶模态的模态加速度响应曲线Fig. 7 Modal acceleration response of first two order modes

OrderRMSofmodaldisplacementresponse13.26×10-623.26×10-8

4.3 飞行器动载荷反算

根据理论模型计算得到的一阶模态剪力和弯矩,根据模态位移响应99%概率的估计值,计算得到飞行器分布动载荷分布结果,包括动剪力和动弯矩反算结果,分别如图8和图9所示.

图8 动剪力反算结果Fig. 8 Identification result of shear

图9 动弯矩反算结果Fig. 9 Identification result of moment

5 结论

对于高超声速飞行器动载荷设计,缺少飞行过程中的动载荷辨识,因此无法直接获得其对飞行器飞行载荷的贡献比例,也无法验证设计时的计算结果；另外,对于高超声速飞行器,其飞行过程中的动载荷测量较为复杂,往往加速度响应较容易获取,因此如何利用常规的加速度测量结果进行动载荷辨识需要进行研究.针对该问题,本文提出了一种基于工作模态辨识的飞行器动载荷辨识的方法,其无需采用动应变测量和标定等特殊方法,利用最易获取也是最普遍测量的飞行振动遥测数据,采用环境激励模态辨识方法辨识其各时刻的模态,再利用振动响应的模态叠加原理和模态正交理论,获取各时刻飞行器低阶模态的响应,再结合模态剪力和模态弯矩进行动载荷反算.

通过研究,可以得出以下结论：

a)低阶模态振型辨识结果与理论计算结果较为一致；

b)通过本方法可以获得飞行器动载荷分布反算结果,为获得动载荷占飞行载荷比例以及验证动载荷设计方法提供支撑.

1 周盼,张权,率志君,李玩幽. 动载荷识别时域方法的研究现状与发展趋势. 噪声与振动控制,2014,34(1):6～11 (Zhou P, Zhang Q, Shuai Z J, Li W Y. Review of Research and Development Status of Dynamic Load Identification in Time Domain.NoiseandVibrationControl,2014,34(1):6～11 (in Chinese))

2 Kazemi M, Hematiyan M R,Ghavami K. An efficient method for dynamic load identification based on structural response. In: EngOpt 2008-International Conference on Engineering Optimization. Rio de Janeiro,Brazil, 2008:1～5

3 程良彦,宋振华,王志华,马宏伟. 基于小波反卷积的冲击力识别. 暨南大学学报(自然科学版),2008,29(5):443～446 (Cheng L Y, Song Z H, Wang Z H, Ma H W. The identification of impact force using a wavelet deconvolution method.JournalofJinanUniversity(NaturalScienceEdition), 2008,29(5):443～446 (in chinese))

4 常晓通,闫云聚. 一种基于遗传算法的时域内稳态动载荷识别方法. 应用力学学报,2012,29(3):274～277 (Chang X T, Yan Y J.A Genetic Algorithm based menthod for Steady-state Dynamic Force Identification in Time Domain.ChineseJournalofAppliedMechanics,2012,29(3):274～277 (in Chinese))

5 杨帆,张方,姜金辉. 正交小波级数拟合法的薄板结构分布动载荷识别技术. 北京理工大学学报,2014,34(6):561～569 (Yang F, Zhang F, Jiang J H. Orthogonal Wavelet Series Fitting in Identification of Distributed Load for Thin Plate Model.TransactionsofBeijingInstituteofTechnology, 2014,34(6):561～569 (in Chinese))

6 Peeters B, De Roeck G, et al. Stochastic Subspace Techniques Applied to Paremeter Identification of Civil Engineering Structures. In: Proceeding of New Advances in Modal Synthesis of Large Structures:Nonlinear, Damped and Nondeterministic Cases,1995:151～162

7 Ibrahim S R. Efficient Random Decrement Computation for Identification of Ambient

8 Abdelgham M,Goursat M,et al. On-Line Modal Monitoring of Aircraft Structures Under Unknow Excitation.MechanicalSystemsandSignalProcessing,1999,13(6):839～853

9 De Roeck G, et al. Benchmark Study on System Identification through Ambient Vibration Measurements. In: 18th IMAC,2000:1106～1112

10Brincker R, et al. Modal Identification from Ambient Responses Using Frequency Domain Decomposition.In: 18th IMAC,2000:625～630

11Bonato B,Ceraavolo R,A De Stefano. Use of Cross-Time-Frequency Estimators for Structural Identification in Non-Stationary Conditions and Under Unknowm Excitation.JournalofSoundandVibration,2000,237(5):775～791

*The project supported by the Research Fund of State Key Laboratory of State Key Laboratory (MCMS-0115G01), National defense technology foundation for scientific research (JSZL2015203B002)

† Corresponding author E-mail: wangliang@nuaa.edu.cn

21 March 2016,revised 20 April 2016.

STUDY ON DYNAMIC LOAD IDENTIFICATION OF AIRCRAFT BASED ON OPERATIONAL MODE ANALYSIS*

Wang Hongbo1,2Zhao Changjian1,2Liao Xuanping1,2Wang Liang2†

(1.NationalUniversityofDefenseTechnology,CollegeofAerospaceScienceandEngineering,Changsha410073,China)(2.ChinaAcademyofLaunchVehicleTechnology,Beijing100076,China)

Based on the operational mode analysis of the aircraft and fighting vibration data, the dynamic load identification technology is investigated in this paper. Firstly, the theory of the ERA method is introduced. Secondly, the dynamic model of the aircraft is put forward. Thirdly, the strategy for implementation is put forward, where the notices are analyzed in detail. Eventually, an example case is studied, where the operational modes are identified based on the telemetry data including the mode frequency and shape. The mode response are then obtained by mode superposition method and mode orthogonal theory, before the dynamic load is derived.

mode identification, ERA, operational mode, dynamic load, vibration, load identification

*国家重点实验室2015年开放课题(MCMS-0115G01)；国防技术基础科研项目(JSZL2015203B002)

10.6052/1672-6553-2016-044

2016-03-21收到第1稿,2016-04-20收到修改稿.

† 通讯作者 E-mail: wangliang@nuaa.edu.cn