聶雪媛，劉中玉，楊國偉

（中國科學(xué)院 力學(xué)研究所 流固耦合系統(tǒng)力學(xué)重點實驗室，北京 100190）

對飛行器氣動彈性不穩(wěn)定性預(yù)測在設(shè)計中非常重要。隨計算流體力學(xué)（Computational Fluid Dynamics，CFD）、計算結(jié)構(gòu)力學(xué)（Computational Structure Dynamics，CSD）及計算機(jī)性能的發(fā)展，直接耦合CFD／CSD的數(shù)值模擬方法成為非線性氣動彈性分析中可信度最高方法［1］。然而在氣動彈性研究中，非定常氣動力計算占據(jù)大部分計算時間，使基于CFD的氣動彈性模擬難以在實際工程應(yīng)用。如何既能保證計算精度又能提高計算效率成為氣動彈性的研究熱點。非定常流場模型降階技術(shù) （ReducedOrder Models，ROMs）的出現(xiàn)為CFD／CSD數(shù)值模擬方法用于工程設(shè)計成為可能。降階模型由全流場全階CFD模型近似投影獲得低階模型，該模型能捕捉非線性流場動力學(xué)特征，計算效率遠(yuǎn)高于CFD。非定常氣動力ROMs代替CFD與結(jié)構(gòu)模型耦合，能進(jìn)行氣動彈性分析或氣動伺服彈性分析及設(shè)計［2－3］。氣動力降階模型基于動態(tài)線性假設(shè)［4］，在給定非線性定常流場基礎(chǔ)上利用小擾動理論推導(dǎo)獲得到線性模型。在基于系統(tǒng)辨識的CFD降階模型方法中，Silva等［5］提出基于脈沖響應(yīng)的Volterra級數(shù)模型降階方法，并與CFL3D結(jié)合進(jìn)行氣動彈性顫振分析。但此方法需獲得系統(tǒng)全部Markov參數(shù)，而所用辨識數(shù)據(jù)因通過“一次激勵結(jié)構(gòu)一個模態(tài)”獲得，會造成結(jié)構(gòu)模態(tài)個數(shù)較多時所需脈沖響應(yīng)個數(shù)呈冪函數(shù)增長。為此，文獻(xiàn)［6－7］開發(fā)出同時激勵結(jié)構(gòu)多個模態(tài)方法提高辨識效率。認(rèn)為輸入多個激勵信號間須存在相位差；但未給出如何確定相位差，需多次調(diào)試。Gupta等［8－10］提出采用飛行試驗中易實現(xiàn)的速度為3211信號作為輸入進(jìn)行系統(tǒng)辨識。利用自回歸滑動平均模型（Autoregressive Moving Average，ARMA）作為替代CFD的降階模型。該方法亦需對各模態(tài)同時輸入的3211信號進(jìn)行反復(fù)試驗、調(diào)整頻率，以期輸出響應(yīng)能含希望激發(fā)的頻帶。為進(jìn)一步提高效率，Raveh［11］用多個白噪聲信號模擬結(jié)構(gòu)多模態(tài)同步輸入，辨識三種多輸入、輸出（Multiple Input Multiple Output，MIMO）模型，即頻域模型、ARMA模型及狀態(tài)空間模型，分析AGARD445.6機(jī)翼在跨聲速區(qū)域的顫振邊界。該方法雖可避免對輸入信號頻譜范圍要求，但仍未涉及多個輸入同步信號間相位差的確定。徐敏等［12－13］以一次激勵一個模態(tài)方式的階躍響應(yīng)信號為輸入，基于Volterra級數(shù)方法建立非定常氣動力降階模型，并用于二維機(jī)翼氣動彈性主動控制及氣動彈性分析。張偉偉等［14］用3211信號建立基于CFD的 ARMA模型用于 Isogai機(jī)翼及AGARD445.6機(jī)翼氣動彈性分析，并認(rèn)為輸入信號需反復(fù)調(diào)試。劉曉燕等［15］采用基于觀測器的Volterra核辨識方法，通過每個模態(tài)單獨(dú)激勵氣動力系統(tǒng)方式，在每個模態(tài)輸入含結(jié)構(gòu)感興趣頻率正弦信號建立氣動力降階模型用于氣動彈性非線性氣動力預(yù)測及顫振抑制。該方法本質(zhì)仍屬基于Volterra級數(shù)的ROMs。

為方便輸入信號選擇、避免同步輸入多信號時反復(fù)測試方可確定相位差、提高辨識計算效率，本文以隨機(jī)白噪聲信號為系統(tǒng)輸入，用一次激勵一個模態(tài)方式同時運(yùn)行多個CFD數(shù)值模擬程序計算系統(tǒng)響應(yīng)。利用所得多組訓(xùn)練數(shù)據(jù)辨識出多個基于ARMA的單輸入多輸出（Single Input Multiple Output，SIMO）部分氣動力降階模型。因非定常氣動力建模建立于動態(tài)線性假設(shè)理論基礎(chǔ)上，故利用線性系統(tǒng)疊加原理獲得基于CFD的多輸入、輸出非定常氣動力降階模型，并將其與結(jié)構(gòu)模型耦合用于AGARD445.6機(jī)翼顫振邊界預(yù)測。通過將預(yù)測結(jié)果與風(fēng)洞試驗結(jié)果比較表明，本文所提建模方法準(zhǔn)確、有效。

1 基于ARMA的氣動力降階模型

顫振邊界預(yù)測為氣動彈性分析中的重要計算及研究彈性體在小擾動激勵下系統(tǒng)響應(yīng)，符合動態(tài)線性模型分類［16］，因此可用線性模型描述非定常氣動力，進(jìn)行氣動彈性顫振分析。

本文所用氣動力降階模型進(jìn)行氣動彈性顫振邊界分析過程可描述為將給定的結(jié)構(gòu)模態(tài)位移一次一個模態(tài)輸入到非定常流場求解器，利用N－S方程求出各階模態(tài)對應(yīng)的廣義氣動力，設(shè)關(guān)心結(jié)構(gòu)的前n個模態(tài)，則會獲得n組“單輸入－n輸出”（記為1－n）的模型訓(xùn)練數(shù)據(jù)，用系統(tǒng)辨識技術(shù)可得n個“1－n”的線性模型，利用線性系統(tǒng)疊加原理，最終獲得“n個輸入－n個輸出”（記為n－n）的氣動力降階模型，見圖1下方虛線框。用所得降階模型替代CFD與結(jié)構(gòu)模型耦合，在結(jié)構(gòu)振幅小擾動下通過改變動壓，可快速獲得結(jié)構(gòu)位移響應(yīng)，從而迅速找到某馬赫數(shù)下的顫振臨界參數(shù)，見圖1上方虛線框。此過程所耗時間與CFD／CSD耦合計算相比可忽略。

圖1 基于氣動力辨識模型的氣動彈性分析Fig.1 Aeroelastic analysis based on aerodynamic force identification model

本文所用氣動力降階模型為自回歸滑動平均ARMA模型，即時間域離散差分方程，形式為

式中：ya為n個廣義氣動力輸出組成的列向量；u為單模態(tài)位移輸入信號標(biāo)量；na，nb為輸出、輸入延遲階數(shù)；k為離散時間變量；A，B為欲辨識系統(tǒng)模型參數(shù)。

設(shè)結(jié)構(gòu)模態(tài)數(shù)為n（以下n均為結(jié)構(gòu)前n階模態(tài)），則式（1）所描述的為圖1下方虛線框中每個“單模態(tài)位移輸入－n個模態(tài)廣義氣動力輸出”的1－n ARMA模型。據(jù)文獻(xiàn)［14］取4～8階延遲模型即可精確描述非定常效應(yīng)，因此利用MATLAB的arx函數(shù)可較快辨識最好的表征訓(xùn)練數(shù)據(jù)（式（1））最優(yōu)的1－n ARMA模型。獲得1－n ARMA模型后用文獻(xiàn)［17］方法在Matlab中編程實現(xiàn)將離散差分方程轉(zhuǎn)換為離散狀態(tài)空間方程形式。第i階模態(tài)激勵下狀態(tài)空間方程為

上述系數(shù)矩陣均為分塊矩陣。式（2）為第i階模態(tài)激勵下所得“單輸入－n輸出”的SIMO系統(tǒng)狀態(tài)空間方程。為便于氣動彈性分析，需獲得“n個模態(tài)輸入激勵下n個模態(tài)廣義氣動力輸出”的多輸入多輸出MIMO系統(tǒng)狀態(tài)空間模型。因此對n個 “單輸入－n輸出”狀態(tài)空間方程線性疊加，獲得所需MIMO狀態(tài)空間方程為

各參數(shù)中上標(biāo)為模態(tài)階次。

2 激勵信號

模型辨識主要為激勵信號選取，要求所用激勵信號能激起被辨識系統(tǒng)所有頻段信號。隨機(jī)信號為含頻帶范圍最寬的信號，理論上用其作為激勵信號能獲得系統(tǒng)所有頻段范圍內(nèi)響應(yīng)。為實際氣動彈性工程應(yīng)用，僅研究某特定頻率附近模態(tài)運(yùn)動所致氣動力，對隨機(jī)信號進(jìn)行帶通濾波處理。本文主要研究AGARD445.6機(jī)翼前四階模態(tài) 9.60 Hz、38.16 Hz、48.35 Hz、91.54 Hz引起的廣義氣動力。對輸入信號進(jìn)行濾波，保證激勵信號包含該頻段即可。

輸入信號類型確定后，數(shù)據(jù)采樣頻率及長度將決定辨識精度。由于該激勵信號用于CFD計算，故其采樣時間與CFD計算時間步相同，本文為無量綱時間0.05，對應(yīng)有量綱頻率 fs為 11 523.6 Hz（馬赫數(shù)為0.96時），其 Nyquist頻率為5 761 Hz，遠(yuǎn)高于所需最高頻率91.54 Hz。對白噪聲激勵信號采用帶通方式濾波，保留頻段范圍為0.000 1fs～0.1fs。采樣長度據(jù)頻率分辨率定義需在CFD上運(yùn)行約2 500個時間步，即可辨識系統(tǒng)最低頻率9.60 Hz。因此可獲得精度較高模型。雖運(yùn)行2 500個時間步需耗費(fèi)一定計算時間，但與頻繁修改參數(shù)進(jìn)行CFD計算及氣動彈性分析的傳統(tǒng)方法相比，計算效率改善仍相當(dāng)可觀。

圖2 機(jī)翼近壁面網(wǎng)格圖Fig.2 Mesh grid of near the wing

3 算例與分析

為驗證該辨識方法的有效性及準(zhǔn)確性，本文對AGARD445.6進(jìn)行模型辨識，利用辨識所得模型進(jìn)行顫振邊界、顫振速度預(yù)測，并與風(fēng)洞試驗結(jié)果比較。計算網(wǎng)格單元數(shù)約40萬。近壁面網(wǎng)格見圖2。CFD網(wǎng)格計算中首層網(wǎng)格厚度約為根弦長的8×10－5，結(jié)構(gòu)振型歸一化見文獻(xiàn)［18］。采用非定常 N－S方程計算機(jī)翼彈性運(yùn)動時非定常流場，空間采用基于結(jié)構(gòu)網(wǎng)格技術(shù)的有限體積法，選HLLE格式進(jìn)行空間離散，用代數(shù)無限插值法實現(xiàn)網(wǎng)格運(yùn)動。用雙時間步進(jìn)行時間離散，擬時間步為隱式格式LUSGS方法。施加模態(tài)位移激勵時需保證產(chǎn)生的結(jié)構(gòu)實際位移、速度滿足線性擾動假設(shè)。

本文采用高斯白噪聲信號，在輸入CFD前先進(jìn)行濾波處理。所選帶通范圍 0.000 1fs～0.1fs，其中 fs為CFD計算用采樣頻率。對應(yīng)有量綱帶通約1.152～1 152 Hz。理論上為辨識出最低頻率，需至少2 500個時間步的CFD計算結(jié)果作為模型訓(xùn)練數(shù)據(jù)，但實際應(yīng)用中的ARMA模型，用200個時間步足夠，多余數(shù)據(jù)可用于校驗辨識的模型。由于針對結(jié)構(gòu)前四階模態(tài)，因此只要辨識四個“單輸入－4輸出”的ARMA模型。模型辨識結(jié)果階次na＝4，nb＝5。

圖3為馬赫數(shù)0.96、輸入第一階模態(tài)位移為過濾高斯白噪聲信號時辨識模型與CFD計算所得氣動力見圖3，其中虛線為模型計算結(jié)果，實線為CFD計算結(jié)果。兩種方法計算結(jié)果吻合。

圖4為馬赫數(shù)0.96、給定第二階模態(tài)位移50 Hz正弦信號時辨識所得ARMA模型計算的氣動力與CFD直接計算的氣動力結(jié)果比較見圖4，其中虛線為模型計算結(jié)果，實線為CFD計算結(jié)果，可見二者吻合較好。

圖3 第一階模態(tài)激勵下辨識模型與CFD計算結(jié)果對比Fig.3 Comparison of ARMA model simulated responses and CFD response to excitation of the 1st mode

圖4 第二階模態(tài)激勵下辨識模型與CFD計算結(jié)果比較Fig.4 Comparison of ARMA model simulated response and CFD response to sinusoidal excitation of the 2nd mode

圖5 輸入3211信號Fig.5 3211 Excitation of the first four modes

三、四階模態(tài)激勵的模型準(zhǔn)確性檢驗與一、二階類似。在檢驗辨識出的四個“單輸入－4輸出”ARMA模型準(zhǔn)確性基礎(chǔ)上對其分別進(jìn)行狀態(tài)空間轉(zhuǎn)換獲得（式（2））方程再進(jìn)行線性疊加，最終獲得“4輸入－4輸出”氣動力狀態(tài)空間模型。對該空間模型精確度采用3211信號進(jìn)行檢驗。輸入具有一定相位差的4個3211信號見圖5。在圖5信號輸入下“4輸入－4輸出”狀態(tài)空間模型計算的廣義氣動力與CFD直接計算結(jié)果比較見圖6，其中實線為CFD計算結(jié)果，方塊為“4輸入－4輸出”狀態(tài)空間模型計算結(jié)果。由圖6看出，MIMO系統(tǒng)狀態(tài)空間模型計算結(jié)果與CFD直接計算結(jié)果吻合較好為評價辨識模型計算精度，本文用常用評價模型準(zhǔn)確度參數(shù)即模型擬合度為

圖6 模型輸出與CFD計算結(jié)果比較Fig.6 Comparison of identification model simulated results and CFD results

式中：為CFD計算結(jié)果；為辨識模型計算結(jié)果；m為數(shù)據(jù)長度；η＝1為ARMA模型與N－S方程計算結(jié)果完全一致。

白噪聲信號、50 Hz正弦信號及3211信號作為結(jié)構(gòu)模態(tài)位移輸入時，ARMA的模型擬合度見表1。由表1看出，所辨識模型精度較高。

表1 不同輸入信號激勵下模型擬合度Tab.1 Comparison of model fitness to different excitation signals

采用辨識模型預(yù)測的AGARD445.6機(jī)翼顫振頻率邊界及速度邊界見圖7、圖8。由二圖看出，耦合結(jié)構(gòu)的N－S方程及ROM計算結(jié)果與風(fēng)洞試驗結(jié)果均較吻合，且機(jī)翼顫振速度在跨聲速區(qū)域出現(xiàn)特有的“凹坑”現(xiàn)象。表明本文所用非定常氣動力模型辨識方法能用于跨聲速氣動彈性穩(wěn)定性分析。

圖7 AGARD 445.6機(jī)翼無因次顫振頻率邊界Fig.7 Nondimensional flutter frequency boundary of AGARD 445.6 wing

圖8 AGARD 445.6機(jī)翼無因次顫振速度邊界Fig.8 Nondimensional flutter speed boundary of AGARD 445.6 wing

基于系統(tǒng)辨識模型進(jìn)行氣動彈性分析效率較高。如尋找馬赫數(shù)為0.96時的顫振臨界參數(shù)。采用求解非定常N－S方程方法搜尋一次參數(shù)至少需CFD運(yùn)行2000個周期，若每個周期含20次內(nèi)迭代（1次內(nèi)迭代約3 s），運(yùn)行時間約需2000 min（Intel Xeon 2.4 GHz雙CPU），而改變參數(shù)重新計算仍費(fèi)時2 000 min。用非定常氣動力建模技術(shù)時以產(chǎn)生2 500個辨識用數(shù)據(jù)而言，需CFD運(yùn)行2 500 min，而參數(shù)辨識時間在Matlab下約需10 s。故模型確定后搜尋一次參數(shù)，Matlab僅需1～2 s。

4 結(jié) 論

（1）本文以過濾隨機(jī)白噪聲信號為激勵信號，利用動態(tài)線性假設(shè)及線性系統(tǒng)疊加原理將辨識所得MISO系統(tǒng)ARMA模型線性疊加獲得MIMO系統(tǒng)模型，實現(xiàn)非定常氣動力建模。將辨識所得模型用白噪聲信號、正弦信號、3211信號等進(jìn)行檢驗，并與直接CFD計算結(jié)果比較，從而驗證辨識模型的準(zhǔn)確性。該方法與現(xiàn)有辨識方法相比，具有對激勵信號要求少、易于施加、辨識效率高等優(yōu)點。

（2）將辨識所得時域信號變換為狀態(tài)空間模型用于AGARD445.6機(jī)翼顫振邊界預(yù)測。氣動力降階模型與結(jié)構(gòu)模型耦合計算結(jié)果與CFD／CSD耦合計算結(jié)果及風(fēng)洞試結(jié)果一致，表明本文方法不僅效率顯著提高，能用于氣動彈性穩(wěn)定性分析，并為氣動伺服彈性分析及設(shè)計提供可能。

［1］陳剛，李躍明.非定常流場降階模型及其應(yīng)用研究進(jìn)展與展望［J］.力學(xué)進(jìn)展，2011，41（6）：686－701.CHEN Gang，LI Yueming.Advances and prospects of the reduced order model for unsteady flow and its application［J］.Advances in Mechanics，2011，41（6）：686－701.

［2］Dowell E H.Eigenmode analysis in unsteady aerodynamics：Reducedorder models［J］.AIAA Journal，1996，34（8）：1578－1583.

［3］Lucia D J，Beran P S，Silva W A.Reducedorder modeling：new approaches for computational physics［J］.Progress in Aerospace Sciences，2004，40（1／2）：51－117.

［4］Dowell E H.A modern course in aeroelasticity，3rd Revised and Enlarged Edition［M］.New York：Klewer Academic Publishers，1995.

［5］Silva W A.Bartels R E.Development of reducedorder models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code［J］.Journal of Fluids and Structures，2004，19：729－745.

［6］Silva W A. Simultaneous excitation of multipleinput multipleoutput CFDbased unsteady aerodynamic systems［J］.Journal of Aircraft，2008，45（4）：1267－1274.

［7］Kim T，Hong M，Bhatia K G，et al.Aeroelastic model reduction for affordable computational fluid dynamicsbased flutter analysis［J］.AIAA Journal，2005，43：2487－2495.

［8］Gupta K K，Voelker L S，Bach C，et al.CFDbased aeroelastic analysis of the X43 hypersonic flight vehicle［C］.Proceedings of the 39th Aerospace Sciences Meeting and Exhibit，2001.

［9］Cowan T J，Gupta K K.Accelerating CFDbased aeroelastic predictions using systems identification［C］.36th AIAA Atmoshperic Flight Mechanics Conference and Exhibit，1998：85－93.

［10］Cowan T J，Gupta K K. Development of discretetime aerodynamic model for CFDbased aeroelastic analysis［C］.Proceedings of the 37th Aerospace Sciences Meeting and Exhibit，1999.

［11］Raveh D E. Identification of computationalfluiddynamic based unsteady aerodynamic models for aeroelastic analysis［J］.Journal of Aircraft，2004（41）：620－632.

［12］徐敏，陳剛，陳士櫓.基于非定常氣動力低階模型的氣動彈性主動控制律設(shè)計［J］.西北工業(yè)大學(xué)學(xué)報，2004，22（6）：748－752.XU Min，CHEN Gang，CHEN Shilu.Design method of active control law based on reducedorder model of unste ady aerodynamics for aeroelasticity［J］.Journal of Northwestern Polytechnical University，2004，22（6）：748－752.

［13］姚偉剛，徐敏.基于Volterra級數(shù)降階模型的氣動彈性分析［J］.宇航學(xué)報，2008，29（6）：1711－1716.YAO Weigang，XU Min.Aeroelastic numerical analysis via Volterra series approach［J］.Journal of Astronautics，2008，29（6）：1711－1716.

［14］張偉偉，葉正寅.基于非定常氣動力辨識技術(shù)的氣動彈性數(shù)值模擬［J］.航空學(xué)報，2006，27（4）：579－583.ZHANG Weiwei，YE Zhengyin.Numerical simulation of aeroelasiticity basing on identification technology of unsteady aerodynamic loads［J］.Journal of Astronautics，2006，27（4）：579－583.

［15］劉曉燕.基于非定常氣動力降階模型的氣動彈性研究［D］.北京：北京航空航天大學(xué)，2011：54－61.

［16］Dowell E H，Crawley E F.A modern course in aeroelasticity［M］.Kluwer Academic，1989.

［17］Cowan T J， Gupta K K. Development of discretetime aerodynamic model for CFDbased aeroelastic analysis［R］.AIAA，1999.

［18］Carson Y E.Agard standard aeroelastic configurations for dynamic responseiwing 445.6［C］.AGARD－R－756，1988.