劉南，郭承鵬，白俊強(qiáng)

1. 中國(guó)航空工業(yè)空氣動(dòng)力研究院 氣動(dòng)研究與試驗(yàn)二部，沈陽(yáng) 110034 2. 高速高雷諾數(shù)氣動(dòng)力航空科技重點(diǎn)實(shí)驗(yàn)室，沈陽(yáng) 110034 3. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

顫振是一種典型的氣動(dòng)彈性不穩(wěn)定性現(xiàn)象，可能導(dǎo)致飛行器的結(jié)構(gòu)遭到災(zāi)難性破壞，在飛行包線內(nèi)絕對(duì)不允許發(fā)生顫振現(xiàn)象且留有一定的裕度。因此，準(zhǔn)確預(yù)測(cè)飛行器顫振特性對(duì)結(jié)構(gòu)設(shè)計(jì)和強(qiáng)度校核具有重要意義，傳統(tǒng)工程顫振計(jì)算采用基于諧振蕩假設(shè)的偶極子格網(wǎng)法，在跨聲速階段該方法的精度大幅下降。然而，飛行器的跨聲速顫振速壓邊界遠(yuǎn)低于其他速域，而且激波運(yùn)動(dòng)及附面層分離常引發(fā)各種復(fù)雜的非線性顫振現(xiàn)象，顯著地影響飛行安全和任務(wù)性能。

針對(duì)飛行器跨聲速顫振問(wèn)題，計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)方法能夠較準(zhǔn)確地預(yù)測(cè)跨聲速流場(chǎng)，因此耦合非定常CFD方法和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(Computational Structural Dynamics，CSD)成為跨聲速顫振分析的重要手段[1]。但是，三維流場(chǎng)的空間自由度通常達(dá)到106量級(jí)以上，利用CFD方法計(jì)算非定常氣動(dòng)力需要耗費(fèi)大量的計(jì)算資源和計(jì)算時(shí)間。例如波音公司的Hong等[2]利用CFD/CSD耦合方法分析HSCT構(gòu)型的顫振特性時(shí)，獲得每個(gè)狀態(tài)的時(shí)域響應(yīng)需要花費(fèi)長(zhǎng)達(dá)108個(gè)CPU小時(shí)。

為了提高分析效率，研究人員提出了若干種氣動(dòng)彈性高效分析手段，主要有Volterra級(jí)數(shù)、本征正交分解和諧波平衡(Harmonic Balance, HB)方法[3]，其中HB方法因其物理意義明確且適用于非線性問(wèn)題等優(yōu)點(diǎn)受到廣泛關(guān)注[4]。Hall等[5]進(jìn)一步改進(jìn)HB方法，不直接求解流場(chǎng)解的諧波系數(shù)，而是通過(guò)離散Fourier變換和逆變換建立諧波系數(shù)和一個(gè)周期內(nèi)等距分布的流場(chǎng)變量之間的關(guān)系，對(duì)后者進(jìn)行求解，不僅提高了分析效率，而且實(shí)現(xiàn)較為簡(jiǎn)單，這種方法通常被稱作高階諧波平衡(High-Order Harmonic Balance, HOHB)方法。該方法公式中顯式地出現(xiàn)運(yùn)動(dòng)頻率，因此對(duì)于頻率已知問(wèn)題，HOHB方法表現(xiàn)優(yōu)異。Ekici等[6]利用無(wú)黏Euler求解器Duke ROTOR建立HOHB分析方法，從而對(duì)旋翼前飛和懸停狀態(tài)進(jìn)行分析，表現(xiàn)出較高的可靠性。Woodgate和Badcock[7]構(gòu)造了一種隱式HOHB求解器，對(duì)NACA0012翼型和F-5機(jī)翼周期性強(qiáng)迫俯仰運(yùn)動(dòng)進(jìn)行分析，與非定常計(jì)算相比，在保證精度的前提下，所需CPU時(shí)間降低10倍以上。陳琦等[8]對(duì)現(xiàn)有計(jì)算程序進(jìn)行改造，實(shí)現(xiàn)HOHB方法，對(duì)跨聲速翼型和超聲速鈍錐周期性強(qiáng)迫俯仰振蕩進(jìn)行分析，結(jié)果表明，二階HOHB方法能夠得到時(shí)域響應(yīng)與非定常計(jì)算非常接近，計(jì)算時(shí)間降低80%。Ronch等[9]采用HOHB方法分析了跨聲速DLR-F12翼身組合體的動(dòng)導(dǎo)數(shù)，并與線性頻域法進(jìn)行對(duì)比，HOHB表現(xiàn)出較高的準(zhǔn)確性。貢伊明等[10]采用帶預(yù)處理的廣義極小殘差算法提高了HOHB方法的計(jì)算收斂性。

對(duì)于運(yùn)動(dòng)頻率未知問(wèn)題，例如圓柱繞流的漩渦脫落頻率、顫振頻率、極限環(huán)振蕩頻率等，需要添加額外的條件。其中一種是相位固定條件，通過(guò)相位對(duì)頻率的梯度進(jìn)行分析。利用這種方法，Thomas等[11-12]基于HOHB方法建立了一種非線性極限環(huán)求解器，對(duì)無(wú)黏NACA64A010翼型和黏性NLR7301翼型跨聲速極限環(huán)振蕩開展研究。Liu和Dowell[13]研究了HDHB方法在帶操縱面極限環(huán)非線性的二維翼型極限環(huán)振蕩問(wèn)題中的應(yīng)用。劉南等[14]研究了高階諧波平衡方法中的非物理解來(lái)源，通過(guò)擴(kuò)充非線性項(xiàng)的子時(shí)間層消除了非物理解且降低了計(jì)算所需的階數(shù)，并將之應(yīng)用于二維兩自由度立方非線性極限環(huán)振蕩問(wèn)題中。但是相位固定條件對(duì)于復(fù)雜問(wèn)題則非常繁瑣和耗時(shí)。Ekici和Hall[15]提出了一種基于殘差L2范數(shù)最小的頻率更新方法，在單自由度非線性極限環(huán)振蕩的數(shù)值模擬中得到較好的結(jié)果。Yao和Marques[16]指出Ekici方法在多自由度問(wèn)題中收斂性變差，并提出改進(jìn)措施，用于分析二維翼型和三維機(jī)翼的無(wú)黏極限環(huán)振蕩。但是基于殘差L2范數(shù)最小的頻率更新方法較為依賴給定的頻率初值，可靠性較差。

針對(duì)顫振頻率未知的問(wèn)題，利用HOHB方法結(jié)合結(jié)構(gòu)振型首先快速地獲得廣義氣動(dòng)力影響系數(shù)(generalized Aerodynamic force Influence Coefficient, AIC)矩陣，從而獲得系統(tǒng)阻尼和來(lái)流速度之間的關(guān)系，建立高效顫振分析流程，避免了在HOHB方法中進(jìn)行頻率求解的難題，有效地提高了顫振分析效率。

1 顫振時(shí)間推進(jìn)方法

基于CFD/CSD耦合的顫振時(shí)間推進(jìn)方法涉及到非定常流場(chǎng)求解、結(jié)構(gòu)有限元求解以及氣動(dòng)/結(jié)構(gòu)數(shù)據(jù)傳遞等。介紹本文所采用的流場(chǎng)求解、結(jié)構(gòu)運(yùn)動(dòng)方程求解、數(shù)據(jù)傳遞和流程建立等部分，最終利用AGARD445.6國(guó)際標(biāo)摸驗(yàn)證分析方法的精度和可靠性。

1.1 CFD求解方法

計(jì)算坐標(biāo)系下無(wú)源項(xiàng)的三維可壓縮RANS方程可以寫成：

(1)

顫振計(jì)算涉及到氣動(dòng)模型的幾何變形，在非定常計(jì)算物理時(shí)間步后均需調(diào)用網(wǎng)格變形子程序更新空間網(wǎng)格。目前，研究人員提出多種網(wǎng)格變形方法，但是這些方法在變形效率、變形能力和變形后網(wǎng)格質(zhì)量方面表現(xiàn)出明顯的優(yōu)缺點(diǎn)[21]。

因此，本文利用背景網(wǎng)格的思想，結(jié)合徑向基函數(shù)[22](Radial Basis Function，RBF)和無(wú)限插值[23](Trans Finite Interpolation，TFI)兩種網(wǎng)格變形方法的優(yōu)點(diǎn)，建立了一套高效高魯棒性網(wǎng)格變形方法[24]。

1.2 結(jié)構(gòu)有限元求解

結(jié)構(gòu)運(yùn)動(dòng)方程可以統(tǒng)一寫成：

(2)

式中：M、D和K分別為結(jié)構(gòu)質(zhì)量陣、阻尼陣和剛度陣；q為結(jié)構(gòu)變形；f為氣動(dòng)力。采用Patran對(duì)結(jié)構(gòu)進(jìn)行建模，通過(guò)Nastran SOL103求解器可以獲得結(jié)構(gòu)振型φm,m=1,2,…及對(duì)應(yīng)的特征頻率ωm,m=1,2,…。僅考慮前M階結(jié)構(gòu)振型，組裝成矩陣：

(3)

結(jié)構(gòu)變形q可以通過(guò)結(jié)構(gòu)振型與廣義位移的線性疊加得到，

(4)

利用結(jié)構(gòu)振型可以建立廣義結(jié)構(gòu)運(yùn)動(dòng)方程：

(5)

式中：廣義質(zhì)量陣、阻尼陣和剛度陣的表達(dá)式為

(6)

廣義力的表達(dá)式為

(7)

1.3 氣動(dòng)/結(jié)構(gòu)數(shù)據(jù)傳遞

能量守恒是氣動(dòng)/結(jié)構(gòu)交界面數(shù)據(jù)傳遞需要滿足的一種重要條件[22]，利用虛功原理得到

(8)

式中：δW為虛功；δu為虛位移；f為力向量，下標(biāo)a代表氣動(dòng)節(jié)點(diǎn)，s為結(jié)構(gòu)節(jié)點(diǎn)。定義一個(gè)耦合矩陣，氣動(dòng)節(jié)點(diǎn)和結(jié)構(gòu)節(jié)點(diǎn)的位移之間的關(guān)系為

ua=Hus

(9)

通過(guò)式(8)和式(9)，力插值公式可以寫成：

fs=HTfa

(10)

通過(guò)RBF方法[22]可以構(gòu)造矩陣H，還可以滿足力和力矩守恒性、插值外形光滑性等要求，此處不再贅述。

對(duì)于結(jié)構(gòu)振型，將結(jié)構(gòu)振型利用矩陣H插值到氣動(dòng)表面得

Φa=HΦs

(11)

式中：Φa、Φs分別是結(jié)構(gòu)振型在氣動(dòng)物面網(wǎng)格點(diǎn)和結(jié)構(gòu)網(wǎng)格點(diǎn)上的表達(dá)。式(11)結(jié)合式(4)和式(7) 即可滿足能量守恒要求。因此，僅需將結(jié)構(gòu)振型插值到氣動(dòng)表面網(wǎng)格上獲得Φa作為計(jì)算輸入文件，在各個(gè)物理時(shí)間步通過(guò)式(12)和式(13) 獲得氣動(dòng)力和位移：

(12)

(13)

1.4 時(shí)域計(jì)算流程

將廣義結(jié)構(gòu)運(yùn)動(dòng)方程轉(zhuǎn)化為狀態(tài)空間方程形式：

(14)

式中：

(15)

圖1 氣動(dòng)彈性時(shí)域分析流程Fig.1 Process of aeroelastic time marching analysis

顫振時(shí)間推進(jìn)方法的流程如圖1所示，圖中展示了流場(chǎng)和結(jié)構(gòu)的時(shí)間推進(jìn)和數(shù)據(jù)傳遞路線，是一種典型的串行處理方法。

1.5 標(biāo)模驗(yàn)證

AGARD445.6機(jī)翼于1961年在NASA蘭利中心TDT(Transonic Dynamics Tunnel)風(fēng)洞進(jìn)行試驗(yàn)[25]，其中三號(hào)軟模型成為顫振分析程序驗(yàn)證的標(biāo)準(zhǔn)模型[26-27]。利用本文建立的顫振分析方法預(yù)測(cè)AGARD445.6機(jī)翼的無(wú)黏和黏性顫振邊界，試驗(yàn)結(jié)果與本文計(jì)算結(jié)果的對(duì)比如圖2所示，其中EXP、Inviscid和Viscous分別代表試驗(yàn)結(jié)果、無(wú)黏計(jì)算結(jié)果和黏性計(jì)算結(jié)果，橫坐標(biāo)為來(lái)流馬赫數(shù)，縱坐標(biāo)Qinf和Freq分別代表來(lái)流動(dòng)壓和振動(dòng)頻率。由圖可見(jiàn)：本文所建立顫振時(shí)間推進(jìn)方法得到的黏性結(jié)果與試驗(yàn)結(jié)果吻合較好，無(wú)黏結(jié)果在低超聲速工況顫振速壓偏高，與參考文獻(xiàn)[26-27]十分吻合，表明了本文方法的可靠性。

圖2 AGARD445.6機(jī)翼顫振邊界結(jié)果對(duì)比Fig.2 Comparison of flutter boundary result of AGARD445.6 wing

2 高效顫振分析方法

2.1 HOHB流場(chǎng)求解方法

式(1)可以寫成

(16)

其中

(17)

如果已知式(16)的解具有周期性且基頻為ω，則Q(t)可以通過(guò)Fourier級(jí)數(shù)展開：

(18)

式中：NH為諧波階數(shù)。將式(18)對(duì)時(shí)間求導(dǎo)：

(19)

式(16)方程右端殘差項(xiàng)R的各階Fourier系數(shù)可以通過(guò)式(20)～式(22)推導(dǎo)或數(shù)值計(jì)算獲得：

(20)

(21)

(22)

再根據(jù)方程式(16)左右兩側(cè)相等的原則，可以得到

(23)

寫成矩陣形式：

(24)

其中

(25)

其中

Q(t2NH)]T

為等距分布的各子時(shí)間層ti,i=0,1,2,…,2NH上的時(shí)域解，ti分布為

(26)

矩陣E為離散Fourier變換矩陣：

(27)

將式(25)代入式(24)，得到：

(28)

式(28)兩邊同乘以E-1：

(29)

直接將式(29)中第2項(xiàng)簡(jiǎn)化為

(30)

其中

R(t2NH)]T

將式(30)代入式(29)，可得

(31)

式中：U=E-1JE。引入偽時(shí)間導(dǎo)數(shù)項(xiàng)驅(qū)動(dòng)式(31) 殘差趨于0，構(gòu)造求解方程：

(32)

2.2 頻域結(jié)構(gòu)運(yùn)動(dòng)方程

對(duì)于顫振等周期性問(wèn)題，模態(tài)位移可以表示為

(33)

(34)

將式(33)和式(34)代入結(jié)構(gòu)運(yùn)動(dòng)方程式(5)可得：

(35)

(36)

其中

假設(shè)結(jié)構(gòu)模態(tài)位移的運(yùn)動(dòng)形式為

(37)

式中：p=g+jω，ω為振蕩頻率(為了簡(jiǎn)單起見(jiàn)，本文沒(méi)有采用減縮頻率k，而是直接采用頻率ω)，g為系統(tǒng)阻尼，是一種更準(zhǔn)確的物理衰減率的近似。將式(37)代入廣義結(jié)構(gòu)運(yùn)動(dòng)方程式(5)可得

(38)

根據(jù)式(38)，可得

(39)

(40)

其中式(39)是系統(tǒng)靜變形方程，可以獲得系統(tǒng)靜平衡位置，式(40)就是顫振頻域控制方程。將式(36)代入(40)，可得

(41)

參考Chen等提出的g方法，引入變量：

(42)

將上述變換代入式(41)，可得

(43)

綜合式(42)和式(43)得到

Ly=gy

(44)

式中：

L=

至此，就通過(guò)廣義力影響系數(shù)矩陣將顫振邊界問(wèn)題轉(zhuǎn)化為復(fù)特征值求解問(wèn)題。具體的分析流程如下：

1) 利用HOHB求解器快速獲得一系列頻率{ωi,i=1,2,…,Nω}下的廣義氣動(dòng)力影響系數(shù)矩陣F(jω)。

2) 在某頻率ωi下，得到一系列來(lái)流速壓Qi,i=1,2,…,NQ下的矩陣L，并采用復(fù)矩陣特征值分析得到不同來(lái)流速壓下的特征值。

3) 根據(jù)特征值隨來(lái)流速度的變化趨勢(shì)，得到特征值虛部為0時(shí)特征值實(shí)部g和來(lái)流速壓Q。

4) 得到特征值實(shí)部g和來(lái)流速壓Q隨頻率的變化趨勢(shì)，當(dāng)特征值實(shí)部g=0時(shí)，對(duì)應(yīng)的ω和Q分別就是顫振臨界頻率和動(dòng)壓。

上述就是本文建立的跨聲速顫振高效預(yù)測(cè)方法。綜上可見(jiàn)，本文避免了直接求解顫振頻率ω的難題，而是利用HOHB方法獲得多個(gè)頻率下的廣義氣動(dòng)力系數(shù)矩陣，通過(guò)頻域法計(jì)算顫振邊界。

3 二維翼型顫振邊界預(yù)測(cè)

對(duì)于僅具有沉浮和俯仰兩自由度的二維翼型，廣義力和模態(tài)位移之間的關(guān)系為

(45)

圖3 Ma=0.80時(shí)通過(guò)強(qiáng)迫簡(jiǎn)諧俯仰運(yùn)動(dòng)得到的廣義氣動(dòng)力影響系數(shù)Fig.3 Generalized aerodynamic force influence coefficients obtained by harmonic pitch motion at Ma=0.80

圖4 Ma=0.85時(shí)通過(guò)強(qiáng)迫簡(jiǎn)諧俯仰運(yùn)動(dòng)得到的廣義氣動(dòng)力影響系數(shù)Fig.4 Generalized aerodynamic force influence coefficients obtained by harmonic pitch motion at Ma=0.85

Ma=0.80和0.85時(shí)沉浮和俯仰模態(tài)阻尼系數(shù)隨來(lái)流速度的變化趨勢(shì)如圖5所示，橫坐標(biāo)是來(lái)流速度，縱坐標(biāo)是模態(tài)阻尼系數(shù)。在Ma=0.80時(shí)，沉浮模態(tài)是顫振臨界模態(tài)，黏性對(duì)阻尼曲線影響較??；在Ma=0.85時(shí)，無(wú)黏結(jié)果在V=33.5 m/s附近沉浮模態(tài)的阻尼系數(shù)變成正數(shù)，系統(tǒng)發(fā)生顫振，然后在V=140 m/s附近阻尼系數(shù)又重新變成負(fù)數(shù)，系統(tǒng)又回到穩(wěn)定狀態(tài)，最后在V=148.8 m/s 附近俯仰模態(tài)的阻尼變成正數(shù)，系統(tǒng)再次發(fā)生顫振，共有3個(gè)臨界點(diǎn)；而黏性結(jié)果僅在V=128 m/s附近沉浮模態(tài)阻尼變成正數(shù)，表明系統(tǒng)發(fā)生顫振。

本文建立的HOHB高效顫振分析方法與時(shí)間推進(jìn)方法得到的無(wú)黏和黏性顫振邊界對(duì)比分別如圖6和圖7所示，其中Time marching和HOHB分別是時(shí)間推進(jìn)和HOHB兩種方法的計(jì)算結(jié)果，橫坐標(biāo)是來(lái)流馬赫數(shù)，圖6(a)和圖6(b)中縱坐標(biāo)分別是顫振速度和頻率。由圖可見(jiàn)，本文方法具有較高的預(yù)測(cè)精度，能夠準(zhǔn)確捕捉無(wú)黏的“S”形顫振邊界[30]、高馬赫數(shù)(Ma≥0.90)時(shí)的“鎖定”現(xiàn)象[31]以及黏性結(jié)果中Ma≥0.85時(shí)由激波誘導(dǎo)附面層分離導(dǎo)致的顫振速度增大現(xiàn)象[32]。

圖5 Ma=0.80和0.85時(shí)模態(tài)阻尼系數(shù)隨來(lái)流速度的變化趨勢(shì)Fig.5 Variation of modal damping coefficients with freestream velocity at Ma=0.80 and Ma=0.85

兩個(gè)方法計(jì)算所用時(shí)間對(duì)比如表1所示。由表可見(jiàn)，HOHB方法所用時(shí)間約為時(shí)間推進(jìn)方法的1/6，但是值得注意的是表中估算是較為保守的，僅采用5次非定常計(jì)算很難得到較為準(zhǔn)確的顫振邊界，而計(jì)算10個(gè)頻率下的廣義力影響系數(shù)矩陣偏多。

圖6 時(shí)間推進(jìn)方法和HOHB法計(jì)算得到的無(wú)黏顫振速度和頻率對(duì)比Fig.6 Comparison of inviscid flutter velocities and frequencies obtained with time marching and HOHB methods

圖7 時(shí)間推進(jìn)方法和HOHB法計(jì)算得到的黏性顫振速度和頻率對(duì)比Fig.7 Comparison of viscous flutter velocities and frequencies obtained with time marching and HOHB methods

方法計(jì)算時(shí)間加速比時(shí)間推進(jìn)5×33.95=169.75HOHB10×2×1.43+0.58=29.185.82

4 三維機(jī)翼顫振邊界預(yù)測(cè)

然后以三維Goland+機(jī)翼為研究對(duì)象，具體氣動(dòng)外形和結(jié)構(gòu)參數(shù)參見(jiàn)Beran等[33]。選取前四階結(jié)構(gòu)模態(tài)(包括一階彎曲模態(tài)、一階扭轉(zhuǎn)模態(tài)、二階彎曲模態(tài)和二階扭轉(zhuǎn)模態(tài))，頻域內(nèi)廣義力和模態(tài)位移之間的關(guān)系為

(46)

令各個(gè)結(jié)構(gòu)模態(tài)做一系列頻率的簡(jiǎn)諧運(yùn)動(dòng)qm(t)=0.001sin(ωit),m=1,2,3,4，從而獲得廣義力影響系數(shù)矩陣，然后利用復(fù)特征值分析得到模態(tài)阻尼隨來(lái)流速度的變化趨勢(shì)，如圖8～圖10所示，其中Mode1、Mode2、Mode3和Mode4分別是前1至4階結(jié)構(gòu)模態(tài)阻尼系數(shù)，橫坐標(biāo)是來(lái)流速度，縱坐標(biāo)是模態(tài)阻尼，分別是無(wú)外掛物構(gòu)型的無(wú)黏阻尼系數(shù)矩陣以及外掛物構(gòu)型的無(wú)黏和黏性阻尼系數(shù)矩陣。結(jié)果表明：Ma=0.80時(shí)顫振臨界模態(tài)是一階彎曲模態(tài)，而Ma=0.92時(shí)顫振臨界模態(tài)由一階彎曲模態(tài)轉(zhuǎn)移至一階扭轉(zhuǎn)模態(tài)，這是由于激波引起壓心后移導(dǎo)致的。

圖8 無(wú)外掛物Goland+機(jī)翼無(wú)黏阻尼系數(shù)隨來(lái)流速度的變化趨勢(shì)Fig.8 Variation of inviscid damping coefficients with freestream velocity of Goland+ wing without store

圖11至圖13是顫振速度和頻率對(duì)比，其中Time marching是時(shí)間推進(jìn)方法計(jì)算結(jié)果，HOHB nmodes=4和HOHB nmodes=2分別是考慮前四階模態(tài)和僅考慮前兩階模態(tài)的HOHB方法計(jì)算結(jié)果，橫坐標(biāo)是來(lái)流馬赫數(shù)。在時(shí)間推進(jìn)方法中，結(jié)構(gòu)模態(tài)數(shù)量對(duì)計(jì)算效率基本沒(méi)有影響(在顫振時(shí)間推進(jìn)流程中廣義結(jié)構(gòu)運(yùn)動(dòng)方程的計(jì)算消耗相比非定常流場(chǎng)計(jì)算微乎其微)，但是在HOHB方法中需要使各個(gè)結(jié)構(gòu)模態(tài)做簡(jiǎn)諧運(yùn)動(dòng)，

圖11 無(wú)外掛物Goland+機(jī)翼無(wú)黏顫振邊界對(duì)比Fig.11 Comparison of inviscid flutter boundaries of Goland+ wing without store

圖12 帶外掛物Goland+機(jī)翼無(wú)黏顫振邊界對(duì)比Fig.12 Comparison of inviscid flutter boundaries of Goland+ wing with store

圖13 帶外掛物Goland+機(jī)翼黏性顫振邊界對(duì)比Fig.13 Comparison of viscous flutter boundaries of Goland+ wing with store

從而獲得廣義力影響系數(shù)矩陣，計(jì)算效率與模態(tài)數(shù)成反比。結(jié)果表明：① 考慮所有結(jié)構(gòu)模態(tài)的HOHB法結(jié)果與時(shí)間推進(jìn)方法吻合很好，僅考慮前兩階模態(tài)也能夠較好地預(yù)測(cè)跨聲速顫振邊界的變化趨勢(shì)；② 對(duì)于不同的構(gòu)型，高階模態(tài)的影響有所不同，例如無(wú)外掛物構(gòu)型，忽略二階彎曲和扭轉(zhuǎn)模態(tài)后HOHB法預(yù)測(cè)的顫振速度平均誤差在1.0%以內(nèi)，而帶外掛物構(gòu)型平均誤差在3.0%左右；③ HOHB法所用時(shí)間僅為時(shí)間推進(jìn)方法的1/6，能夠顯著提高顫振邊界預(yù)測(cè)效率。

5 結(jié) 論

1) 基于高精度CFD方法和結(jié)構(gòu)模態(tài)運(yùn)動(dòng)方程耦合求解，建立了顫振時(shí)間推進(jìn)方法，通過(guò)對(duì)AGARD445.6機(jī)翼標(biāo)摸進(jìn)行測(cè)試，計(jì)算與試驗(yàn)結(jié)果吻合較好，表現(xiàn)出較高的計(jì)算精度和可靠性。

2) 對(duì)于顫振等頻率未知的問(wèn)題，很難利用HOHB方法直接求解。采取一種間接的方式——利用HOHB方法快速地獲得廣義氣動(dòng)力系數(shù)矩陣，進(jìn)而求解顫振邊界，建立高效顫振頻域分析流程。

3) 基于HOHB建立高效顫振頻域計(jì)算方法，并將之應(yīng)用于二維和三維跨聲速顫振邊界預(yù)測(cè)，與時(shí)間推進(jìn)方法計(jì)算結(jié)果吻合較好，表現(xiàn)出較高的效率和可靠性。