高速飛行器熱結構工作時變模態參數辨識

周思達， 劉莉， 李昱霖， 周小陳

北京理工大學 宇航學院， 北京 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