基于高階諧波平衡的跨聲速顫振高效預測方法

劉南，郭承鵬，白俊強

1. 中國航空工業空氣動力研究院 氣動研究與試驗二部，沈陽 110034 2. 高速高雷諾數氣動力航空科技重點實驗室，沈陽 110034 3. 西北工業大學 航空學院，西安 710072

顫振是一種典型的氣動彈性不穩定性現象，可能導致飛行器的結構遭到災難性破壞，在飛行包線內絕對不允許發生顫振現象且留有一定的裕度。因此，準確預測飛行器顫振特性對結構設計和強度校核具有重要意義，傳統工程顫振計算采用基于諧振蕩假設的偶極子格網法，在跨聲速階段該方法的精度大幅下降。然而，飛行器的跨聲速顫振速壓邊界遠低于其他速域，而且激波運動及附面層分離常引發各種復雜的非線性顫振現象，顯著地影響飛行安全和任務性能。

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

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

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

針對顫振頻率未知的問題，利用HOHB方法結合結構振型首先快速地獲得廣義氣動力影響系數(generalized Aerodynamic force Influence Coefficient, AIC)矩陣，從而獲得系統阻尼和來流速度之間的關系，建立高效顫振分析流程，避免了在HOHB方法中進行頻率求解的難題，有效地提高了顫振分析效率。

1 顫振時間推進方法

基于CFD/CSD耦合的顫振時間推進方法涉及到非定常流場求解、結構有限元求解以及氣動/結構數據傳遞等。介紹本文所采用的流場求解、結構運動方程求解、數據傳遞和流程建立等部分，最終利用AGARD445.6國際標摸驗證分析方法的精度和可靠性。

1.1 CFD求解方法

計算坐標系下無源項的三維可壓縮RANS方程可以寫成：

(1)

顫振計算涉及到氣動模型的幾何變形，在非定常計算物理時間步后均需調用網格變形子程序更新空間網格。目前，研究人員提出多種網格變形方法，但是這些方法在變形效率、變形能力和變形后網格質量方面表現出明顯的優缺點[21]。

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

1.2 結構有限元求解

結構運動方程可以統一寫成：

(2)

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

(3)

結構變形q可以通過結構振型與廣義位移的線性疊加得到，

(4)

利用結構振型可以建立廣義結構運動方程：

(5)

式中：廣義質量陣、阻尼陣和剛度陣的表達式為

(6)

廣義力的表達式為

(7)

1.3 氣動/結構數據傳遞

能量守恒是氣動/結構交界面數據傳遞需要滿足的一種重要條件[22]，利用虛功原理得到

(8)

式中：δW為虛功；δu為虛位移；f為力向量，下標a代表氣動節點，s為結構節點。定義一個耦合矩陣，氣動節點和結構節點的位移之間的關系為

ua=Hus

(9)

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

fs=HTfa

(10)

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

對于結構振型，將結構振型利用矩陣H插值到氣動表面得

Φa=HΦs

(11)

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

(12)

(13)

1.4 時域計算流程

將廣義結構運動方程轉化為狀態空間方程形式：

(14)

式中：

(15)

圖1 氣動彈性時域分析流程Fig.1 Process of aeroelastic time marching analysis

顫振時間推進方法的流程如圖1所示，圖中展示了流場和結構的時間推進和數據傳遞路線，是一種典型的串行處理方法。

1.5 標模驗證

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

圖2 AGARD445.6機翼顫振邊界結果對比Fig.2 Comparison of flutter boundary result of AGARD445.6 wing

2 高效顫振分析方法

2.1 HOHB流場求解方法

式(1)可以寫成

(16)

其中

(17)

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

(18)

式中：NH為諧波階數。將式(18)對時間求導：

(19)

式(16)方程右端殘差項R的各階Fourier系數可以通過式(20)～式(22)推導或數值計算獲得：

(20)

(21)

(22)

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

(23)

寫成矩陣形式：

(24)

其中

(25)

其中

Q(t2NH)]T

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

(26)

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

(27)

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

(28)

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

(29)

直接將式(29)中第2項簡化為

(30)

其中

R(t2NH)]T

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

(31)

式中：U=E-1JE。引入偽時間導數項驅動式(31) 殘差趨于0，構造求解方程：

(32)

2.2 頻域結構運動方程

對于顫振等周期性問題，模態位移可以表示為

(33)

(34)

將式(33)和式(34)代入結構運動方程式(5)可得：

(35)

(36)

其中

假設結構模態位移的運動形式為

(37)

式中：p=g+jω，ω為振蕩頻率(為了簡單起見，本文沒有采用減縮頻率k，而是直接采用頻率ω)，g為系統阻尼，是一種更準確的物理衰減率的近似。將式(37)代入廣義結構運動方程式(5)可得

(38)

根據式(38)，可得

(39)

(40)

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

(41)

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

(42)

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

(43)

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

Ly=gy

(44)

式中：

L=

至此，就通過廣義力影響系數矩陣將顫振邊界問題轉化為復特征值求解問題。具體的分析流程如下：

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

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

3) 根據特征值隨來流速度的變化趨勢，得到特征值虛部為0時特征值實部g和來流速壓Q。

4) 得到特征值實部g和來流速壓Q隨頻率的變化趨勢，當特征值實部g=0時，對應的ω和Q分別就是顫振臨界頻率和動壓。

上述就是本文建立的跨聲速顫振高效預測方法。綜上可見，本文避免了直接求解顫振頻率ω的難題，而是利用HOHB方法獲得多個頻率下的廣義氣動力系數矩陣，通過頻域法計算顫振邊界。

3 二維翼型顫振邊界預測

對于僅具有沉浮和俯仰兩自由度的二維翼型，廣義力和模態位移之間的關系為

(45)

圖3 Ma=0.80時通過強迫簡諧俯仰運動得到的廣義氣動力影響系數Fig.3 Generalized aerodynamic force influence coefficients obtained by harmonic pitch motion at Ma=0.80

圖4 Ma=0.85時通過強迫簡諧俯仰運動得到的廣義氣動力影響系數Fig.4 Generalized aerodynamic force influence coefficients obtained by harmonic pitch motion at Ma=0.85

Ma=0.80和0.85時沉浮和俯仰模態阻尼系數隨來流速度的變化趨勢如圖5所示，橫坐標是來流速度，縱坐標是模態阻尼系數。在Ma=0.80時，沉浮模態是顫振臨界模態，黏性對阻尼曲線影響較小；在Ma=0.85時，無黏結果在V=33.5 m/s附近沉浮模態的阻尼系數變成正數，系統發生顫振，然后在V=140 m/s附近阻尼系數又重新變成負數，系統又回到穩定狀態，最后在V=148.8 m/s 附近俯仰模態的阻尼變成正數，系統再次發生顫振，共有3個臨界點；而黏性結果僅在V=128 m/s附近沉浮模態阻尼變成正數，表明系統發生顫振。

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

圖5 Ma=0.80和0.85時模態阻尼系數隨來流速度的變化趨勢Fig.5 Variation of modal damping coefficients with freestream velocity at Ma=0.80 and Ma=0.85

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

圖6 時間推進方法和HOHB法計算得到的無黏顫振速度和頻率對比Fig.6 Comparison of inviscid flutter velocities and frequencies obtained with time marching and HOHB methods

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

方法計算時間加速比時間推進5×33.95=169.75HOHB10×2×1.43+0.58=29.185.82

4 三維機翼顫振邊界預測

然后以三維Goland+機翼為研究對象，具體氣動外形和結構參數參見Beran等[33]。選取前四階結構模態(包括一階彎曲模態、一階扭轉模態、二階彎曲模態和二階扭轉模態)，頻域內廣義力和模態位移之間的關系為

(46)

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

圖8 無外掛物Goland+機翼無黏阻尼系數隨來流速度的變化趨勢Fig.8 Variation of inviscid damping coefficients with freestream velocity of Goland+ wing without store

圖11至圖13是顫振速度和頻率對比，其中Time marching是時間推進方法計算結果，HOHB nmodes=4和HOHB nmodes=2分別是考慮前四階模態和僅考慮前兩階模態的HOHB方法計算結果，橫坐標是來流馬赫數。在時間推進方法中，結構模態數量對計算效率基本沒有影響(在顫振時間推進流程中廣義結構運動方程的計算消耗相比非定常流場計算微乎其微)，但是在HOHB方法中需要使各個結構模態做簡諧運動，

圖11 無外掛物Goland+機翼無黏顫振邊界對比Fig.11 Comparison of inviscid flutter boundaries of Goland+ wing without store

圖12 帶外掛物Goland+機翼無黏顫振邊界對比Fig.12 Comparison of inviscid flutter boundaries of Goland+ wing with store

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

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

5 結 論

1) 基于高精度CFD方法和結構模態運動方程耦合求解，建立了顫振時間推進方法，通過對AGARD445.6機翼標摸進行測試，計算與試驗結果吻合較好，表現出較高的計算精度和可靠性。

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

3) 基于HOHB建立高效顫振頻域計算方法，并將之應用于二維和三維跨聲速顫振邊界預測，與時間推進方法計算結果吻合較好，表現出較高的效率和可靠性。