高速飞行器热结构工作时变模态参数辨识

周思达， 刘莉， 李昱霖， 周小陈

北京理工大学 宇航学院， 北京 100081

高速飞行器热结构工作时变模态参数辨识

周思达*， 刘莉， 李昱霖， 周小陈

北京理工大学 宇航学院， 北京 100081

高速飞行器由于其很高的飞行速度而无可避免地受到气动加热作用的影响，进而引起结构特性的时变。采用理论或有限元方法(FEM)进行数值分析，难以获取反映结构在飞行(工作)状态下的真实模态参数。通过辨识获取高速飞行器热环境下的时变结构模态参数是一项十分具有挑战性的任务。针对此问题，引入参数化时频域的最大似然方法，对气动加热作用下的高速飞行器升力面结构的时变模态参数进行了辨识。通过模拟真实飞行状态的数值算例研究，说明参数化时频域的最大似然方法能够很好地辨识出低信噪比(SNR)情况下的模态频率和模态振型，验证了参数化时频域最大似然方法适用于具有显著时变特征的高速飞行器热结构的时变结构模态参数辨识，可为将来相关的工程研究和应用提供良好的理论支持。

时变； 模态参数； 参数辨识； 高速飞行器； 热结构； 参数化； 最大似然

高速飞行器在速度上的显著优势，因此在应用领域受到了广泛关注。关于高速飞行器的研究已经成为相关领域的热点。当然，高速飞行器也因其高飞行速度给设计和分析带来了众多挑战。如何准确地获取高速飞行器在气动加热环境下的热结构时变模态参数就是挑战之一[1]。首先，高速飞行器在飞行中才会产生相应的气动加热效应，因此，相关的模态分析是工作状态下的；第二，由于热环境随时间改变，高速飞行器结构在热作用下的动力学特性，即模态参数，是随时间变化的。另外，虽然有限元方法(FEM)等数值方法能够对高速飞行器热环境下的结构动力学特性进行预测，但无论是理论还是有限元等数值方法，无法最为真实地反映复杂结构在环境作用下的工作模态参数。这就需要基于实验的工作模态分析手段，来获取结构真实的模态参数[2]。最近，也有学者进行了高速飞行器热结构动力学特性的实验研究。吴振强等[3]利用实验得到了不同温度下的飞行器壁板结构动特性；白云鹤[4]采用非参数化时频域的方法对热结构模态频率进行了辨识。

在基于响应的时变结构模态参数辨识领域，近年来产生了许多新方法，其中，Poulimenos和Fassois[5]提出的泛函序列时间相关自回归滑动平均(FS-TARMA)方法、Zhou等[6]提出的参数化时频域两步最小二乘方法、Xu等[7]提出的基于小波变换的非参数化时频域方法等具有代表性。

本文将利用所提出的参数化时频域最大似然方法[8]对高速飞行器热结构在工作状态下的结构模态参数进行辨识，建立模拟真实情况的数值算例，并通过算例讨论方法在不同测量噪声情况下的性能，说明方法对高速飞行器热结构模态参数辨识的适用性。为将来高速飞行器热结构的实验研究，特别是模态参数辨识提供一定的理论参考和技术支持。

1 时变结构的时频域参数化动力学模型

线性时不变结构的传递函数可利用参数化的公分母模型来表达[9]。假设作用在结构上的激励为白噪声，结构响应的功率谱函数也可以表示为参数化的公分母模型。将时不变结构功率谱函数的公分母模型拓展为时间相关的公分母模型，以表示时变结构的时间相关功率谱函数，如：

(1)

式中：tτ和ωf分别为采样时刻和角频率，下标τ=1,2,…,Nτ表示采样时间，f=1,2,…,Nf表示采样频率；k=1,2,…,NoNr(即k=(i-1)No+o)，No和Nr分别为输出点和参考点数量；分子多项式和分母多项式皆定义为

(2)

(3)

(4)

进一步写为

(5)

为方便以下推导，将式(5)的参数向量写为

(6)

式中：θ为参数化的公分母模型的待定参数向量。

另外，如式(2)所示，分子多项式和分母多项式都含有时频基函数，其数值特性直接影响后续辨识的数值条件，某些情况下甚至影响辨识的成败。本文采用一种正交多项式与z域多项式的混合形式(Orthogonal Polynomial andz-domain Mapping hybrid，OPZM)的时频基函数[8]，即时间多项式采用正交多项式(如Chebyshev多项式、Jacobi多项式等)，频率多项式采用z域多项式：

(8)

式中：ω为角频率；Ts为采样长度。

2 参数化时频域的最大似然辨识

根据频域的最大似然费用函数和时频域时变结构动力学参数化的时间相关公分母模型(式(1))，时频域的最大似然辨识方法的费用函数为

(9)

为得到较大的动态范围，采用以下对数形式的总方程误差：

(10)

矩阵形式的总方程误差可写为

(11)

与线性最小二乘辨识不同，由于如式(10)所示的总方程误差是关于待估参数θ的非线性函数，因此最大似然辨识方法的费用函数在最小化时需求解非线性代数方程，可采用Gauss-Newton法。

根据Gauss-Newton迭代，在第m个迭代步，待估参数的增量Δθm可由如下线性方程求得：

(12)

式中：J为m个迭代步时的Jacobi矩阵。

在每一个迭代步后对待估参数θ进行更新，如：

θm+1=θm+Δθm

(13)

待估参数θ的初始值由参数化时频域的最小二乘辨识给出[10]。利用更新的待估参数θ，可以得到更新的Jacobi矩阵和更新的总方程误差(式(10))。对应Jacobi矩阵(如式(12)定义)的分块矩阵分别为

(14)

(15)

其中，

(16)

(17)

根据式(14)和式(15)重写式(12)，得

(18)

式(18)中的前NoNr个方程可写为

(19)

消去Δβk，由式(18)可得

(20)

求解方程式(20)可得当前迭代步中待估参数α的增量Δα，并利用式(19)回代得到待估参数β的增量Δβ。按照式(6)将Δα和Δβ写成待估参数θ的增量Δθ。每一个迭代步都是一个求解缩减的正则方程的过程，迭代的终止条件为前后两个迭代步的所有响应的功率谱函数估计值与测量值之间相对误差的最大值之差小于ζ，定义为

(21)

一旦得到待估参数向量θ，便可计算出各个时刻的模态参数，具体方法参考文献[11]。

3 气动加热下的升力面结构模态参数辨识

本文提出的方法能够针对工作状态下的结构，辨识其时变的结构模态参数。下面通过一个数值仿真算例，说明其针对高速飞行器结构在时变热环境下，如何辨识结构随时间变化的模态参数。

为实现与“冻结”的结构进行比较，从而验证所提方法，算例将不考虑结构变形对气动力和热传导的影响，结构受到热的影响是单向的。另外，不考虑热应力导致的结构微分刚度对结构整体刚度的影响，仅考虑材料受温度的影响。这些假设在实验中都是难以满足的，但是对于方法的验证至关重要。当然，以上忽略的效应对结构的影响，在实际中可以包含在结构特性随时间变化的整体效应中，并且本文采用的参数化时频域最大似然辨识方法基于结构响应，反映结构与外部因素作用的宏观模态特性，而不依赖复杂耦合问题内在机理的建模，因此，本文数值算例的假设仍不失一般性。

本文算例参考文献[12]，建立如图1所示的升力面。

图1 升力面外形尺寸Fig.1 Profile of lifting surface

3.1 气动加热计算

对于高速飞行器，气动加热是其显著特点。为模拟真实状态下升力面的气动加热效应，采用高精度CFD计算方法，获取结构表面的温度场边界，并用于后续计算。假设飞行器的飞行状态曲线如图2所示(速度和攻角曲线)[13-14]。

图2 飞行状态曲线[13-14]Fig.2 Flight state curves[13-14]

在一系列时刻点，利用Fastran气动加热分析软件进行CFD计算仿真，获得升力面上下表面的温度场分布，参考文献[14]。以10 s时刻为例，上表面(迎风面)的温度分布如图3所示。

图3 t=10 s时刻升力面上表面的温度分布Fig.3 Temperature distribution on upper face of lifting surface at t=10 s

3.2 分布温度边界下的瞬态热传导计算

为得到升力面结构在上下表面分布温度边界情况下的结构各处温度，采用有限元方法进行瞬态热传导分析。用于瞬态热传导的有限元模型如图4所示，单元为四节点的热传导单元(DS4)。

图4 升力面有限元模型Fig.4 Finite element model of lifting surface

由于在表面处，CFD和有限元的网格并不一一对应，因此，需利用CFD的计算结果(以图3为例)，在结构上下表面有限元节点上分别进行二维插值，从而得到各个时刻(前20 s)上下表面的分布温度边界，结构材料属性如表1所示。

表1 结构材料属性

由于材料的热传导系数和热容系数与温度无关，因此，在材料特性和分布温度边界条件已知的情况下，利用线性有限元的瞬态热传导分析，计算得到结构各节点处的温度分布。翼稍处某节点温度随时间的变化曲线如图5所示。

图5 翼稍处某节点温度随时间的变化曲线Fig.5 Temperature vs time at a certain point of wing tip

3.3 状态“冻结”下的参考热模态计算

由于材料的弹性模量与温度相关，导致结构刚度与温度相关，因此，为验证本文采用的参数化时频域的最大似然辨识方法，利用一系列状态“冻结”的结构获取其模态参数与前者比较。状态“冻结”指某时刻下给定由热传导分析得到的分布温度场，边界条件为固定翼根部(固定节点的平动和转动自由度)。另外，为得到显著的时变特性，本文将结构的弹性模量设置为一个较大的范围，如表1所示。

结构有限元模型网格与热传导有限元模型一致，如图4所示。单元类型为四节点减缩积分的Mindlin壳单元，厚度方向积分点个数为5。利用有限元模态分析(Lanczos特征值提取方法)，获取各时刻下的结构模态。模态频率如图6所示；10 s时刻第一、二阶模态振型如图7所示。

图6 状态“冻结”下的各时刻模态频率Fig.6 “Frozen”-configuration modal frequency at different time instants

图7 10 s时刻的第一、二阶模态振型Fig.7 Mode shapes for the first and second modes at 10 s

3.4 温度变化下的结构动力学响应计算

为辨识提供数据，采用基于有限元方法的结构瞬态响应分析，获取结构的响应信号。结构边界条件与有限元模态分析一致。激励点位于如图4所示的三角指示处，激励方向为沿坐标z轴方向，激励形式为白噪声(辨识时假设未知)。利用Newmark-β法计算响应(Newmark-β积分参数γ=0.5，β=0.25)，积分步长为0.000 5s(约为第二阶模态频率最大值的37倍)。

为模态参数辨识提供响应数据，分别在结构表面的6个位置(如图4中的圆点所示)处z方向“测量”速度响应，采样频率为200 Hz。为模拟与真实情况较为贴近的测量噪声污染的情况，对“测量”的响应加入信噪比(SNR)分别为10 dB (10%测量噪声)和3 dB (50%测量噪声)的噪声。信噪比分别为10 dB和3 dB时，响应点1处的响应如图8所示；信噪比分别为10 dB和3 dB时，响应点1处的时间相关自功率谱函数(采用平滑伪Wigner-Ville分布(SPWVD)计算[15])如图9所示。

图8 响应点1处的速度响应历程Fig.8 Velocity history at point 1

图9 响应点1处的时间相关自功率谱函数Fig.9 Time-dependent auto power spectra at point 1

3.5 工作状态热结构模态参数辨识

3.4节通过有限元的瞬态响应分析，得到了“测量”的结构响应。这些响应将为结构模态参数辨识提供最基本的数据。

利用第2节的参数化时频域的最大似然方法进行辨识。首先，利用时频分析(此处用SPWVD)估计各响应的时间相关功率谱函数，并且为了估计振型，使用响应的互功率谱，即以响应点1为参考点，响应点1～6分别与其计算互谱，如图9所示(自谱为自身的互谱的特例)。

图10 时变模态频率的辨识结果Fig.10 Estimated modal frequency of time-varying structure

图11 辨识得到的10 s时刻的模态振型Fig.11 Estimated shapes for the first and second modes at 10 s

如图10所示，采用参数化时频域的最大似然方法，能够很好地辨识出在气动加热环境下，结构时变的模态频率，甚至在信噪比较低(50%噪声，信噪比为3 dB)的情况下，除很少的间断外，辨识值和有限元方法得到的理论参考值十分吻合。另外，比较图7和图11，发现辨识得到的模态振型也十分接近于理论振型，取得较好的结果，同样可以得到其余任意时刻的模态振型，在此限于篇幅，仅给出10 s时刻的模态振型。

4 结 论

本文利用参数化时频域的最大似然方法对气动加热作用下的高速飞行器升力面结构的时变模态参数进行了辨识。不同于传统短时的方法，需要将响应数据分段进行一系列时不变模态参数辨识，本方法在时频域上整体地一次辨识出结构随时间变化的模态参数，数值算例说明其不仅能够很好地辨识出低信噪比情况下的模态频率，还能准确地辨识出结构的模态振型，并验证了参数化时频域最大似然辨识方法适用于具有显著时变特征的高速飞行器热结构的时变结构模态参数辨识。另外，为实现对比，本文采用的算例仅考虑了气动热作用对结构的单向作用，但是辨识方法本身是基于结构响应的，辨识结果可反映结构与外部因素作用的宏观模态特性，能够包括结构与气动和热耦合等。因此，本方法仍可为实际的的工程研究和应用提供良好的理论支持。

[1] Yu Y H, Li S, Wang Y D. Study on dynamic characteristics of metallic thermal protection system with thermal environ[J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(1): 171-175 (in Chinese). 余艳辉, 李书, 王远达. 热环境下的金属热防护系统的动力学特性研究[J]. 振动测试与诊断, 2013, 33(1): 171-175.

[2] Zhang L, Brincker R, Andersen P. An overview of operational modal analysis: major development and issues[C]∥1st International Operational Modal Analysis Conference (IOMAC), 2005： 179-190.

[3] Wu Z Q, Cheng H, Zhang W, et al. Effects of thermal environment on dynamic properties of aerospace vehicle[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(2): 334-342 (in Chinese). 吴振强, 程昊, 张伟, 等. 热环境对飞行器壁板结构动特性的影响[J]. 航空学报, 2013, 34(2): 334-342.

[4] Bai Y H. Research of structural modal experiment and parameter identification for high speed vehicles under flying environments[D]. Harbin: Harbin Institute of Technology, 2013 (in Chinese). 白云鹤. 高速飞行器飞行环境下结构模态试验及参数识别方法研究[D]. 哈尔滨: 哈尔滨工业大学, 2013.

[5] Poulimenos A G, Fassois S D. Output-only stochastic identification of a time-varying structure via functional series TARMA models[J]. Mechanical Systems and Signal Processing, 2009, 23(4): 1180-1204.

[6] Zhou S D, Heylen W, Sas P, et al. Parametric modal identification of time-varying structures and the validation approach of modal parameters[J]. Mechanical Systems and Signal Processing, 2014, 47(1): 94-119.

[7] Xu X, Shi Z Y, You Q. Identification of linear time-varying systems using a wavelet-based state-space method[J]. Mechanical Systems and Signal Processing, 2012, 26(6): 91-103.

[8] Zhou S D, Heylen W, Sas P, et al. Maximum likelihood estimator of operational modal analysis for linear time-varying structures in time-frequency domain[J]. Journal of Sound and Vibration, 2014, 333(11): 2339-2358.

[9] Verboven P. Frequency-domain system identification for modal analysis[D]. Brussels: Vrije Universiteit Brussel, 2002.

[10] Zhou S D, Liu L, Yang W, et al. Matrix fraction polynomial model-based least square estimation of modal parameters for linear time-varying structures[J]. Journal of Vibration and Shock, 2014, 33(6): 118-123 (in Chinese). 周思达, 刘莉, 杨武, 等. 基于矩阵分式多项式时变结构模态参数最小二乘辨识[J]. 振动与冲击, 2014, 33(6): 118-123.

[11] Zhou S D. Study on theory and experiment of time-frequency-domain modal parameter estimation for linear time-varying mechanical structures[D]. Beijing: Beijing Institute of Technology, 2012 (in Chinese). 周思达. 线性时变结构时频域模态参数辨识理论及实验研究 D]. 北京: 北京理工大学, 2012.

[12] McNamara J J, Friedmann P P, Powell K G, et al. Aeroelastic and aerothermoelastic behavior in hypersonic flow [J]. AIAA Journal, 2008, 46(10): 2591-2610.

[13] Zhu H, Liu L, Zhou S, et al. Integrated aerodynamic thermal structure design optimization method of lifting surfaces [J]. Journal of Aircraft, 2012, 49(5): 1521-1526.

[14] Li Y L. Studies on aerodynamic-structural-thermal multidisciplinary design optimization and method of flight vehicles[D]. Beijing: Beijing Institute of Technology, 2014 (in Chinese). 李昱霖. 气动热结构多学科分析及高效优化策略研究 [D]. 北京: 北京理工大学, 2014.

[15] Roshan-Ghias A, Shamsollahi M B, Mobed M, et al. Estimation of modal parameters using bilinear joint time-frequency distributions[J]. Mechanical Systems and Signal Processing, 2007, 21(5): 2125-2136.

Tel: 010-68918752

E-mail: zhousida@bit.edu.cn

*Corresponding author. Tel.： 010-68918752 E-mail: zhousida@bit.edu.cn

Operational identification of time-varying modal parameters for thermal structures of high-speed aerial vehicles

ZHOU Sida*, LIU Li, LI Yulin, ZHOU Xiaochen

SchoolofAerospaceEngineering,BeijingInstituteofTechnology,Beijing100081,China

The high-speed aerial vehicles are unavoidably impacted by the aerodynamic heating effect due to their high speed, which may further result in the time-varying characteristics of their structures. The theoretical and finite element method (FEM)-based numerical approaches hardly acquire the real modal parameters of the in-flight (operational) aerial vehicles. Focusing on this problem, the parametric time-frequency-domain maximum likelihood method is introduced into the application of the high-speed aerial vehicles and the time-varying modal parameters are estimated by this parametric time-frequency-domain maximum likelihood method. Based on a high-fidelity simulation example, it is proved that the introduced method can identify the modal frequency and mode shapes even in the low signal noise ratio (SNR) cases and the method is suitable for the modal parameter estimation of thermal structures of high-speed aerial vehicles with notable time-varying natures, which can theoretically support the relative engineering studies and applications.

time-varying; modal parameter; parameter identification; high-speed aerial vehicles; thermal structures； parametric; maximum likelihood

2014-07-01； Revised： 2014-09-10； Accepted： 2014-09-23； Published online： 2014-09-26 10:41

s： National Natural Science Foundation of China (11372036, 11402022); Beijing Institute of Technology Foundation for Basic Research (20120142009)

2014-07-01; 退修日期： 2014-09-10; 录用日期： 2014-09-23; 网络出版时间： 2014-09-26 10：41

www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0267.html

国家自然科学基金(11372036， 11402022)； 北京理工大学基础研究基金(20120142009)

Zhou S D, Liu L, Li Y L, et al. Operational identification of time-varying modal parameters for thermal structures of high-speed aerial vehicles[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 373-380. 周思达, 刘莉, 李昱霖, 等．高速飞行器热结构工作时变模态参数辨识[J]．航空学报, 2015, 36(1): 373-380.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2014.0267

V214.1； TH113.1

A

1000-6893(2015)01-0373-08

周思达 男， 博士， 讲师。主要研究方向：飞行器结构动力学、时变结构动力学建模与辨识。

*通讯作者.Tel.： 010-68918752 E-mail: zhousida@bit.edu.cn

URL： www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0267.html