基于ARMA/ROM的顫振邊界模擬和快速預(yù)測(cè)

馬艷峰，賀爾銘，李俊杰，曾憲昂，唐長(zhǎng)紅

(1.西北工業(yè)大學(xué) 航空學(xué)院，陜西 西安710072;2.西安飛機(jī)設(shè)計(jì)研究所 強(qiáng)度設(shè)計(jì)研究所，陜西 西安710089;3.西安飛機(jī)設(shè)計(jì)研究所 總師辦，陜西 西安710089)

0 引言

使用CFD技術(shù)對(duì)氣動(dòng)彈性進(jìn)行分析能夠顯著提升計(jì)算精度和應(yīng)用范圍。CFD技術(shù)在時(shí)空維上對(duì)氣體流動(dòng)的刻畫(huà)愈細(xì)致，則非定常氣動(dòng)力系統(tǒng)的維數(shù)就愈高，直接導(dǎo)致了氣動(dòng)彈性研究計(jì)算量大、分析周期長(zhǎng)的結(jié)果。因此，發(fā)展非定常流降階模型(ROM)已成為一個(gè)較為活躍的研究領(lǐng)域，基于CFD技術(shù)的非定常氣動(dòng)力降階模型，可以快速求解系統(tǒng)對(duì)任意輸入的響應(yīng)，不必重新運(yùn)行CFD代碼［1］。

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

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

1 ARMA模型降階方法

基于動(dòng)態(tài)線性化假設(shè)，以顫振為代表的氣動(dòng)彈性穩(wěn)定性問(wèn)題中，非定常氣動(dòng)力子系統(tǒng)可視為線性時(shí)不變系統(tǒng)。因此，可以將氣動(dòng)彈性計(jì)算中的非定常CFD求解器視為一個(gè)待辨識(shí)的動(dòng)力系統(tǒng)，以結(jié)構(gòu)的模態(tài)位移作為輸入，廣義氣動(dòng)力作為輸出。

1.1 過(guò)濾脈沖信號(hào)

使用脈沖信號(hào)激勵(lì)對(duì)系統(tǒng)的動(dòng)態(tài)特性進(jìn)行充分激勵(lì)，選擇了物理上易于實(shí)現(xiàn)的多正弦疊加信號(hào)。全階計(jì)算的氣動(dòng)彈性方法用于計(jì)算脈沖信號(hào)激勵(lì)后的訓(xùn)練數(shù)據(jù)。信號(hào)的表達(dá)式為:

對(duì)于AGARD445.6來(lái)說(shuō)，脈沖信號(hào)的參數(shù)選擇為:α0= －0.1，ωc=3.05，τi=0。

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

圖1 濾波脈沖信號(hào)輸入Fig.1 Staggered sequence of FIM input signals

圖2 廣義氣動(dòng)力響應(yīng)Fig.2 General aerodynamic force responses

1.2 非定常氣動(dòng)力降階模型

以顫振為代表的氣動(dòng)彈性問(wèn)題中，非定常氣動(dòng)力模型可表示為離散的線性狀態(tài)空間，其差分方程表達(dá)式為:

在非定常流體動(dòng)力學(xué)中，廣義位移作為系統(tǒng)輸入，廣義氣動(dòng)力作為系統(tǒng)輸出。對(duì)于結(jié)構(gòu)的N階模態(tài)，結(jié)構(gòu)的運(yùn)動(dòng)可表示為:

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

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

1.3 狀態(tài)空間模型

1.3.1 ARMA矩陣和變換

使用Hankel矩陣法將ARMA模型轉(zhuǎn)換得到一個(gè)不可約的狀態(tài)空間表達(dá)式［13］。式(4)可以用以下的向量寫(xiě)成矩陣形式:

則:

其中:

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

1.3.2 Markov參數(shù)表達(dá)式

任何的線性動(dòng)力系統(tǒng)采樣數(shù)據(jù)都可以寫(xiě)成Markov參數(shù)表達(dá)式的形式。因此，式(6)可以等效寫(xiě)為:

式中:Hj為 Markov 參數(shù)，Hj∈RN×N。

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

可得Markov參數(shù)的計(jì)算表達(dá)式:

1.3.3 連續(xù)時(shí)間的狀態(tài)空間實(shí)現(xiàn)

對(duì)于ARMA模型和Markov參數(shù)模型的狀態(tài)空間實(shí)現(xiàn)并不是唯一的，通過(guò)Hankel矩陣得到狀態(tài)空間的不可約狀態(tài)［13］:

其中:

將Hankel矩陣T進(jìn)行奇異值分解，得到:

據(jù)此得到一個(gè)不可約的離散狀態(tài)空間:

其中，輸入和輸出的矩陣表達(dá)式為:

得到氣動(dòng)力對(duì)模態(tài)位移的響應(yīng)離散狀態(tài)空間模型式(11)后，通過(guò)Laplace微分算子進(jìn)行z變換，可以很方便地獲得等效的連續(xù)時(shí)間狀態(tài)空間模型。

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

1.4 線性顫振閉環(huán)分析

控制方程采用結(jié)構(gòu)動(dòng)力學(xué)方程如下:

為了便于時(shí)域求解式(13)給出的結(jié)構(gòu)動(dòng)力學(xué)方程，引入狀態(tài)變量，則式(13)可寫(xiě)成狀態(tài)空間的形式:

其中:

式中:M，K分別為結(jié)構(gòu)廣義質(zhì)量和剛度矩陣;E，F(xiàn)分別為廣義坐標(biāo)向量和廣義氣動(dòng)力向量。

使用非定常氣動(dòng)力的狀態(tài)空間式(12)，該ARMA/ROM取代CFD求解器與結(jié)構(gòu)動(dòng)力學(xué)狀態(tài)空間(14)進(jìn)行耦合，進(jìn)行高精度的氣動(dòng)彈性分析，采用四階龍格-庫(kù)塔方法求解微分方程組。

2 算例與分析

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

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

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

圖3 ARMA/ROM辨識(shí)結(jié)果Fig.3 Identification results of ARMA/ROM

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

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

對(duì)于亞聲速、跨聲速和超聲速三種工況，ARMA/ROM方法所得到的機(jī)翼顫振特性和CFD全階模擬差別均較小。兩種方法的精度相當(dāng)，這也體現(xiàn)了氣動(dòng)力ROM技術(shù)在復(fù)雜氣動(dòng)彈性問(wèn)題中的應(yīng)用潛力。

通過(guò)比較全階氣動(dòng)彈性計(jì)算方法和氣動(dòng)力辨識(shí)法對(duì)AGARD445.6機(jī)翼的顫振分析結(jié)果，得出如下結(jié)論:

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

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

(3)從兩種方法與試驗(yàn)值的比較來(lái)看，在“凹坑點(diǎn)”之前，本文采用Euler/BL方程辨識(shí)計(jì)算結(jié)果與試驗(yàn)值符合性比N-S方程計(jì)算結(jié)果更好一些;在“凹坑點(diǎn)”之后N-S方程計(jì)算結(jié)果更接近于試驗(yàn)結(jié)果?？傮w而言，二者的計(jì)算精度相差不大，但Euler/BL方程的辨識(shí)計(jì)算效率比N-S方程要高，這也正是本文采用Euler/BL方程作為流場(chǎng)控制方程的主因。

3 結(jié)束語(yǔ)

本文以狀態(tài)空間方程的形式提供了一種低維的流體模型，通過(guò)ARMA/ROM技術(shù)識(shí)別流體模態(tài)，采用系統(tǒng)實(shí)現(xiàn)理論基于Markov參數(shù)獲得模態(tài)系數(shù)相對(duì)應(yīng)的狀態(tài)空間模型。結(jié)果表明，基于ARMA的降階模型能高效準(zhǔn)確地模擬全階氣動(dòng)彈性模型的特征，準(zhǔn)確地預(yù)測(cè)流體系統(tǒng)的動(dòng)態(tài)行為;相對(duì)于全階系統(tǒng)而言，自由度數(shù)量以及計(jì)算時(shí)間都大幅縮減，為下一步進(jìn)行三維機(jī)翼的工程應(yīng)用計(jì)算奠定了基礎(chǔ)。

［1］ 楊超，劉曉燕，吳志剛.基于POD-Observer技術(shù)的非定常氣動(dòng)力建模方法［J］.中國(guó)科學(xué):技術(shù)科學(xué)，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í)數(shù)的非定常氣動(dòng)力降階模型［J］.強(qiáng)度與環(huán)境，2007，34(5):22-28.

［6］ 陳剛，徐敏，陳士櫓.基于Volterra級(jí)數(shù)的非線性非定常氣動(dòng)力降階模型［J］.宇航學(xué)報(bào)，2004，25(5):492-495.

［7］ 崔鵬.基于CFD/CSD的機(jī)翼氣動(dòng)彈性計(jì)算研究［D］.南京:南京航空航天大學(xué)，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］張偉偉，葉正寅.基于氣動(dòng)力降階模型的跨音速氣動(dòng)彈性穩(wěn)定性分析［J］.計(jì)算力學(xué)學(xué)報(bào)，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］劉黨輝.系統(tǒng)辨識(shí)方法及應(yīng)用［M］.北京:國(guó)防工業(yè)出版社，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ī)翼的靜氣動(dòng)彈性研究［J］.航空計(jì)算技術(shù)，2014，44(3):53-57.