基于ARMA/ROM的顫振邊界模擬和快速預測

馬艷峰，賀爾銘，李俊杰，曾憲昂，唐長紅

(1.西北工業大學 航空學院，陜西 西安710072;2.西安飛機設計研究所 強度設計研究所，陜西 西安710089;3.西安飛機設計研究所 總師辦，陜西 西安710089)

0 引言

使用CFD技術對氣動彈性進行分析能夠顯著提升計算精度和應用范圍。CFD技術在時空維上對氣體流動的刻畫愈細致，則非定常氣動力系統的維數就愈高，直接導致了氣動彈性研究計算量大、分析周期長的結果。因此，發展非定常流降階模型(ROM)已成為一個較為活躍的研究領域，基于CFD技術的非定常氣動力降階模型，可以快速求解系統對任意輸入的響應，不必重新運行CFD代碼［1］。

依據非定常CFD解的動態線性化假設，目前已經發展了很多的降階建模技術。其中包括:本征正交分解(POD)［2-3］、Volterra 級數方法［4-6］、神經網格模型、自回歸滑動平均(ARMA)模型方法［7］。ARMA模型作為系統辨識的一種通用模型結構，多用于辨識MIMO系統。ARMA模型形式簡潔，可準確地描述非線性系統，鑒于應用的準確性和便捷性，氣動力的ARMA/ROM在氣動彈性分析方面得到了較多的應用。Cowan等［8］利用 ARMA模型結構基于一組CFD的非定常數值解擬合了線性系統的系數。Gupta等［9］將ARMA/ROM應用于高超聲速飛行器復雜系統的顫振分析，體現了該方法在工程應用方面的潛力。國內方面，張偉偉等［10］進行了基于ARMA氣動力辨識的高效高精度顫振計算研究。

本文系統地給出氣動力ARMA/ROM的實現過程，并將之應用于AGARD445.6機翼的跨聲速顫振研究。非定常氣動力方面，采用位移輸入作為訓練輸入，分別使用Euler和Euler/附面層耦合的方法完成ROM訓練過程并獲得訓練輸出，基于Matlab編程實現ARMA模型的參數估計，并進行氣動力模型的驗證［11-12］。上述氣動力ROM的建立過程簡便可靠，一定程度上體現了工程應用的特點。

1 ARMA模型降階方法

基于動態線性化假設，以顫振為代表的氣動彈性穩定性問題中，非定常氣動力子系統可視為線性時不變系統。因此，可以將氣動彈性計算中的非定常CFD求解器視為一個待辨識的動力系統，以結構的模態位移作為輸入，廣義氣動力作為輸出。

1.1 過濾脈沖信號

使用脈沖信號激勵對系統的動態特性進行充分激勵，選擇了物理上易于實現的多正弦疊加信號。全階計算的氣動彈性方法用于計算脈沖信號激勵后的訓練數據。信號的表達式為:

對于AGARD445.6來說，脈沖信號的參數選擇為:α0= －0.1，ωc=3.05，τi=0。

使用Euler/BL全階氣動彈性方法，通過上述的多正弦疊加信號，逐級激勵前四階結構模態。AGARD445.6機翼(Ma=0.901)的前四階濾波脈沖輸入信號和系統的廣義氣動力的響應如圖1和圖2所示。由圖可見，ARMA/ROM的各階仿真輸入、輸出和CFD直接模擬結果吻合良好。這表明，對于涉及復雜振型的三維機翼，氣動力ARMA/ROM的精度仍然較高。需要說明的是，本文中對廣義氣動力和廣義位移的各變量及過程變量，均進行了無量綱化處理。

圖1 濾波脈沖信號輸入Fig.1 Staggered sequence of FIM input signals

圖2 廣義氣動力響應Fig.2 General aerodynamic force responses

1.2 非定常氣動力降階模型

以顫振為代表的氣動彈性問題中，非定常氣動力模型可表示為離散的線性狀態空間，其差分方程表達式為:

在非定常流體動力學中，廣義位移作為系統輸入，廣義氣動力作為系統輸出。對于結構的N階模態，結構的運動可表示為:

式中:φ為模態向量;η為模態位移。由式(2)和式(3)可得:

式中:q=1，…，N;yp，up分別為廣義氣動力和廣義位移［1］。

1.3 狀態空間模型

1.3.1 ARMA矩陣和變換

使用Hankel矩陣法將ARMA模型轉換得到一個不可約的狀態空間表達式［13］。式(4)可以用以下的向量寫成矩陣形式:

則:

其中:

將式(5)進行z變換可得:

1.3.2 Markov參數表達式

任何的線性動力系統采樣數據都可以寫成Markov參數表達式的形式。因此，式(6)可以等效寫為:

式中:Hj為 Markov 參數，Hj∈RN×N。

由式(6)和式(7)可得［12］:

可得Markov參數的計算表達式:

1.3.3 連續時間的狀態空間實現

對于ARMA模型和Markov參數模型的狀態空間實現并不是唯一的，通過Hankel矩陣得到狀態空間的不可約狀態［13］:

其中:

將Hankel矩陣T進行奇異值分解，得到:

據此得到一個不可約的離散狀態空間:

其中，輸入和輸出的矩陣表達式為:

得到氣動力對模態位移的響應離散狀態空間模型式(11)后，通過Laplace微分算子進行z變換，可以很方便地獲得等效的連續時間狀態空間模型。

其中，輸入輸出矩陣表達式為［12］:

1.4 線性顫振閉環分析

控制方程采用結構動力學方程如下:

為了便于時域求解式(13)給出的結構動力學方程，引入狀態變量，則式(13)可寫成狀態空間的形式:

其中:

式中:M，K分別為結構廣義質量和剛度矩陣;E，F分別為廣義坐標向量和廣義氣動力向量。

使用非定常氣動力的狀態空間式(12)，該ARMA/ROM取代CFD求解器與結構動力學狀態空間(14)進行耦合，進行高精度的氣動彈性分析，采用四階龍格-庫塔方法求解微分方程組。

2 算例與分析

本文使用AGARD445.6機翼進行了顫振計算。顫振分析主要取其前四階模態，即一階彎曲(mode1，9.57 Hz)、一階扭轉(mode2，38.17 Hz)、二階彎曲(mode3，48.35 Hz)和二階扭轉(mode4，91.55 Hz)。

Euler方程可以較準確地模擬跨聲速氣動力，對網格和計算資源的要求也相對適中;Euler/BL方法是以Euler方程為控制方程，通過邊界層方程計算邊界層特性［14-15］求解無粘流場。Euler/BL方程很好地補足了Euler方程無法考慮流體粘性的缺點，不論是定常還是非定常計算都具有很高的精度［14-15］。因此，文中的CFD直接模擬和ROM訓練過程均分別使用以Euler和Euler/BL方程為控制方程進行計算和辨識。

分別采用ARMA/ROM和CFD全階模擬方法對AGARD445.6機翼的顫振特性進行分析，圖3給出了Ma=0.901工況的Euler/BL辨識結果。對降階和全階計算結果與風洞試驗值之間的比較如圖4和圖5所示。

圖3 ARMA/ROM辨識結果Fig.3 Identification results of ARMA/ROM

圖4 全階與降階系統顫振速度比較Fig.4 Comparison for flutter speeds of full-order and reduced-order systems

圖5 全階與降階系統顫振頻率比較Fig.5 Comparison for flutter frequencies of full-order and reduced-order systems

對于亞聲速、跨聲速和超聲速三種工況，ARMA/ROM方法所得到的機翼顫振特性和CFD全階模擬差別均較小。兩種方法的精度相當，這也體現了氣動力ROM技術在復雜氣動彈性問題中的應用潛力。

通過比較全階氣動彈性計算方法和氣動力辨識法對AGARD445.6機翼的顫振分析結果，得出如下結論:

(1)兩種方法計算結果的趨勢一致，AGARD445.6機翼顫振“凹坑點”大約在 Ma=0.954附近;

(2)氣動力辨識法和全階氣動彈性的計算結果基本一致，但在計算效率上氣動力辨識法能提高3～4倍;

(3)從兩種方法與試驗值的比較來看，在“凹坑點”之前，本文采用Euler/BL方程辨識計算結果與試驗值符合性比N-S方程計算結果更好一些;在“凹坑點”之后N-S方程計算結果更接近于試驗結果。總體而言，二者的計算精度相差不大，但Euler/BL方程的辨識計算效率比N-S方程要高，這也正是本文采用Euler/BL方程作為流場控制方程的主因。

3 結束語

本文以狀態空間方程的形式提供了一種低維的流體模型，通過ARMA/ROM技術識別流體模態，采用系統實現理論基于Markov參數獲得模態系數相對應的狀態空間模型。結果表明，基于ARMA的降階模型能高效準確地模擬全階氣動彈性模型的特征，準確地預測流體系統的動態行為;相對于全階系統而言，自由度數量以及計算時間都大幅縮減，為下一步進行三維機翼的工程應用計算奠定了基礎。

［1］ 楊超，劉曉燕，吳志剛.基于POD-Observer技術的非定常氣動力建模方法［J］.中國科學:技術科學，2010，40(8):861-866.

［2］ Hall K C，Thomas JP，Dowell E H.Reduced-order modeling of unsteady small-disturbance flows using afrequencydomain proper orthogonal decomposition technique［R］.AIAA-99-0655，1999.

［3］ Hall K C，Thomas JP，Dowell E H.Proper orthogonal decomposition technique for transonic unsteady aerodynamic flows［J］.AIAA Journal，2000，38(2):1853-1862.

［4］ Silva W A，Raveh D E.Development of unsteady aerodynamic state-space models from CFD-based pulse responses［R］.AIAA-2001-1213，2001.

［5］ 徐敏，李勇，曾憲昂，等.基于Volterra級數的非定常氣動力降階模型［J］.強度與環境，2007，34(5):22-28.

［6］ 陳剛，徐敏，陳士櫓.基于Volterra級數的非線性非定常氣動力降階模型［J］.宇航學報，2004，25(5):492-495.

［7］ 崔鵬.基于CFD/CSD的機翼氣動彈性計算研究［D］.南京:南京航空航天大學，2011.

［8］ Cowan T J，Arena A S.Development of a discrete-time aerodynamic model for CFD-based aeroelastic analysis［R］.AIAA-99-765，1999.

［9］ Gupta K K，Voelker L S.CFD-based aeroelastic analysis of the X-43 hypersonic flight vehicle［R］.AIAA-2001-0712，2001.

［10］張偉偉，葉正寅.基于氣動力降階模型的跨音速氣動彈性穩定性分析［J］.計算力學學報，2007，24(6):768-772.

［11］ Lai K L，Tsai H M.Flutter simulation and prediction with CFD-based reduced-order model［R］.AIAA-2007-731，2007.

［12］ Lai K L，Lum K Y.Flutter-boundary prediction using system identification-based reduced-order aeroelasticity analysis［R］.AIAA-2012-1710，2012.

［13］劉黨輝.系統辨識方法及應用［M］.北京:國防工業出版社，2010:90-94.

［14］ Gao C，Luo S，Liu F，et al.Calculation of airfoil flutter by an Euler method with approximate boundary conditions［R］.AIAA-2003-3830，2003.

［15］馬艷峰，賀爾銘，李俊杰，等.大展弦比機翼的靜氣動彈性研究［J］.航空計算技術，2014，44(3):53-57.