基于非定常氣動力降階的AGARD445.6 硬機翼不同迎角顫振研究

容浩然，戴玉婷，許云濤，楊 超

(北京航空航天大學航空科學與工程學院，北京 100191)

現代戰斗機具備的高機動性，使之需要在大迎角條件下飛行；導彈機動時舵面處于大迎角狀態；飛機在遇到陣風時機翼可能會進入較大的迎角狀態；風力機葉片、渦輪葉片等機械部件通常也是在有穩態迎角狀態工作。非零迎角狀態機翼的顫振特性與零迎角時有一定區別，但傳統的氣動彈性試驗以及計算研究通常采用零迎角，針對三維機翼不同迎角的氣動彈性研究相對較少。ASHLEY[1]、EDWARD 等[2]、DOGGETT 等[3]和YATES 等[4]的研究表明：機翼迎角是跨音速段的顫振特性的重要影響參數，他們利用風洞試驗、理論方法等研究了跨音速段迎角對顫振邊界的影響，但以上研究都屬于小迎角范圍，且迎角對顫振邊界的影響主要是由跨音速段氣動非線性引起；YE 等[5]分別對矩形機翼和三角翼完成了不同迎角的顫振風洞試驗；張偉偉等[6]和全景閣等[7?8]對比不同迎角下大后掠三角翼的氣動彈性時域響應特性，研究了分離渦破裂前和破裂后三角翼的氣動彈性穩定性特性，得到了有益規律；翼型小迎角經典顫振和大迎角失速顫振方面，葉正寅等[9]使用數值計算方法研究NACA0012 翼型在不同迎角下的顫振特性，發現到達某一迎角后，顫振邊界突然下降；劉暢暢等[10]對NACA0012 翼型完成不同迎角的風洞試驗，顫振結果與葉正寅流固耦合計算結果基本一致；LI 等[11]使用數值計算、循環神經網絡降階方法，研究了NACA64A010 二元翼段在跨音速小迎角范圍內不同迎角的顫振邊界；另外，TALLEY 等[12]針對AGARD445.6 軟機翼，基于模態法使用數值計算、降階模型方法分析不同迎角、來流馬赫數下的氣動彈性時域響應，發現迎角增加后，顫振邊界降低。根據現有研究，在超過某一迎角后，機翼的顫振邊界可能下降，因此，有必要對不同迎角下三維機翼的顫振特性進行研究。

氣動彈性研究中，非定常氣動力計算方法有ONERA 模型[13?14]、Leishman-Beddoes(L-B)模型[15]等。隨著計算機技術的興起，計算流體力學(Computational fluid dynamics, CFD)開始應用于各項工程和研究[16?17]。但是基于CFD 的流固耦合(Fluid-structure interaction, FSI)方法的非定常建模通常計算量大、計算時間長，因而需要發展高效、高精度的非定常氣動降階模型(Reduced order model, ROM)，快速、準確計算非定常氣動力。非定常氣動降階方法主要包括系統辨識法和特征提取法，前者以自回歸滑動平均法(AutoRegressive moving average, ARMA)、Volterra 級數、機器學習等為代表，后者以本征正交分解(Proper orthogonal decomposition, POD)、動力學模態分解(Dynamic mode decomposition, DMD)為代表。在各種非定常氣動降階方法中，人工神經網絡(Artificial neural network, ANN)具有很強的非線性擬合能力[18]，且ANN 天然就具備對多輸入/多輸出系統的描述能力，這些優勢使得近年來ANN 在非定常氣動降階領域應用非常廣泛。MARQUES 等[19]建立了多層神經網絡，分別以單俯仰運動為輸入和俯仰運動、來流馬赫數耦合輸入，以翼型升力系數和力矩系數為輸出，辨識跨聲速段不同馬赫數下二維翼型非定常氣動力，由于只考慮了當前時刻輸入與當前輸出之間的關系，本質是一種定常模型；QIU等[20]建立了多層神經網絡，把時間序列數據作為輸入，預測了二維翼型的失速顫振；WINTER 等[21]基于反饋模糊神經網絡建立了AGARD445.6 機翼變馬赫數氣動降階模型，預測了不同馬赫數下給定廣義運動的廣義力，快速計算頻域氣動力系數矩陣，然后使用pk 法研究了機翼跨音速顫振特性；王博斌等[22]建立了帶反饋的徑向基神經網絡(Radial basis function neural network, RBFNN)氣動降階模型，辨識預測跨聲速不同頻率、不同振幅俯仰運動下二元翼段的氣動特性；KOU 等[23]在此基礎上其與自回歸(AutoRegressive eXogenous,ARX)方法并聯得到混合降階模型；LI 等[24]和KOU 等[25]建立了長短時記憶單元(Long short-term memory, LSTM)循 環 神 經 網 絡(Recurrent neural network, RNN)，考慮跨音速條件下，以俯仰、沉浮耦合和馬赫數作為輸入，辨識預測了不同馬赫數下翼段氣動特性，并將氣動降階模型用于顫振特性預測。目前的非定常氣動降階研究中，研究對象大多針對零迎角剛性翼型，針對非零迎角三維彈性機翼的非定常氣動降階研究較少。

本文首先介紹了基于模態疊加的流固耦合方法、基于RBFNN 的非定常氣動力ROM 建模方法和兩種基于降階模型的顫振預測方法；針對AGARD445.6 硬機翼[26]，完成不同初始迎角流固耦合顫振計算，總結了隨著初始迎角增加顫振邊界變化規律；選取10°初始迎角條件，建立基于RBFNN 的非定常氣動力ROM，預測不同速度、減縮頻率的非定常氣動力，預測10°迎角的顫振特性；基于10°迎角非定常氣動降階建模的研究結果，建立考慮初始迎角輸入的非定常氣動力ROM，快速預測機翼不同迎角的顫振特性；最后基于建立的ROM 分析顫振邊界隨迎角變化的機理。

1 計算方法

1.1 基于模態疊加的流固耦合方法

多自由度線性系統受迫振動方程為：

式中：M、C、K分別為系統的質量矩陣、阻尼矩陣、剛度矩陣；z為多自由度系統真實位移；Q為外部作用力。系統真實位移可以表示為廣義坐標的線性疊加：

式中：Φ 為結構質量歸一模態；d為廣義坐標，則式(1)可化為：

式中：I、Cp、Kp分別為廣義質量矩陣(單位陣)、廣義阻尼矩陣、廣義剛度矩陣(對角陣且對角元素為固有角頻率平方)；p0為來流動壓；f為時域廣義力系數。氣動彈性系統中通常忽略結構阻尼，則式(3)的微分方程組解耦。通過施加初始廣義速度擾動d˙0，使用CFD 方法求解氣動力Q，采用四階龍格-庫塔方法求解此微分方程組，即可得到該初始條件下的氣動彈性時間響應。

在流固耦合計算中引入假設：只考慮垂直于機翼所處平面方向的變形，即機翼處于xoy平面，模態只考慮z方向分量。該假設在靜變形和動態變形不太大時是合理的。

在計算具有初始迎角機翼的氣動彈性響應時，應先計算該迎角下的靜氣動彈性變形，在靜氣動彈性計算收斂后施加廣義速度擾動，得到動氣動彈性時間響應，具體步驟如下：

1)給系統施加阻尼，計算式(3)直至收斂得到靜氣彈結果；

2)在靜氣彈的基礎上施加廣義速度擾動d˙0，同時取消阻尼，計算式(3)得到動氣彈時域響應。

RBFNN 由BROOMHEAD 等[27]提出，其結構包含三層：輸入層、隱含層和輸出層，其中隱含層的激活函數為徑向基函數，輸入靠近基函數中心點時神經元才會被激活，輸出層激活函數為線性函數。典型的RBFNN 結構如圖1 所示。

圖1 徑向基神經網絡結構Fig. 1 Structure of RBFNN

1.2 基于徑向基神經網絡的非定常氣動降階方法

式中：隱含層的激活函數一般為高斯基函數，bv為隱含層第v個神經元輸出；cv為第v個神經元的中心向量；σv第v個神經元的寬度。隱含層到輸出層為線性變換：

式中：yc為輸出層第c個輸出；wcv為隱含層第v個神經元到輸出層第c個神經元的權重。RBFNN 訓練過程中，在隱含層神經元個數、基函數、中心向量和寬度確定后，權重的學習僅需求解線性方程，因而RBFNN 相比于傳統多層前饋神經網絡結構簡單、訓練成本低、收斂速度快。另外，只要隱含層數量足夠大，RBFNN 能以任意精度逼近非線性函數。

RBFNN 一般有兩種：正則RBFNN 和廣義RBFNN。正則RBFNN 的隱含層數量與訓練樣本量相同，隱含層中心向量與對應的輸入向量相同，正則RBFNN 通常足夠精確，但考慮本文的問題，輸入樣本量較大導致神經網絡結構龐大，可能導致訓練困難、矩陣出現病態，訓練結果可能過擬合；本文選擇建立的廣義RBFNN，則是逐個添加隱含層神經元，直到均方誤差達到設定要求，建立的模型則為滿足誤差要求的最簡RBFNN 結構，詳細步驟為：

1)用神經網絡擬合訓練樣本，得到訓練樣本中誤差最大的輸入向量；

2)增加隱含層神經元，其中心向量為上述誤差最大的輸入向量；

3)優化權重矩陣，返回步驟1)，直到均方誤差達到設定要求或隱含層神經元數量達到設定上限。

1.3 基于氣動降階模型的顫振預測方法

在有初始迎角的氣動彈性分析中，廣義位移、廣義力均可以表示為靜氣彈量和相對于靜氣彈量的疊加，即 ξ=ξs+ξt，其中：ξ 為廣義位移、廣義力、廣義力系數；下角標s 為靜氣彈量；下角標t 為相對于靜氣彈量。因此，顫振方程式(3)可以表示為：

為了描述方便，若無特別說明，后續所述廣義位移、廣義力、廣義力系數等均為相對于靜氣彈量，并省略下角標t。本文基于非定常氣動力ROM 的顫振預測涉及兩種方法，時域下的龍格庫塔法(RK 法)和頻域下的VG 法[28]。

RK 法中，使用降階模型預測顫振方程式(6)中的廣義力系數f，設定初始條件，使用四階龍格庫塔方法數值求解方程式(6)，即可得到氣動彈性時域響應。

VG 法中，臨界穩定的顫振方程可以化為式(7)。定義減縮頻率和無量綱來流速度，即減縮速度如式(8)：

式中：L為參考長度；ρ 為密度；V為來流速度；ωα為扭轉模態角頻率；μ為質量系數；AIC為頻域氣動力系數矩陣；gs為阻尼；ω為結構運動角頻率。AIC與減縮頻率和來流速度有關，因為機翼在大迎角時速度的改變會改變機翼的靜變形位置，從而對非定常氣動特性產生影響，進而影響AIC矩陣。通過降階模型預測某一減縮頻率、速度下的時域廣義力系數，再將結果轉化到頻域，即可得到該減縮頻率、速度下的AIC矩陣?；诮惦A模型，使用VG 法預測某一迎角顫振特性的步驟為：

1)給定一系列速度輸入V*，計算該速度下的AIC。

2 基于流固耦合的數值算例與結果分析

2.1 流固耦合方法驗證

使用ANSYS Fluent 軟件進行氣動力計算。首先驗證SA 湍流模型計算分離流的適用性，選取一大后掠角三角翼模型VFE-2[29]，三角翼前緣后掠角為65°，平均氣動弦長為0.436 m，展長為0.610 m，模型詳細幾何參數見參考文獻[29]。其前緣可以更換，本文選取帶尖前緣的VFE-2 三角翼模型。計算條件為雷諾數6×106、來流馬赫數0.4、迎角13.3°，網格如圖2 所示，網格數量1.13×106，最大y+為24。計算不同展長截面的壓力系數分布，與試驗結果對比后發現趨勢基本一致，如圖3 所示。

圖2 VFE-2 計算網格Fig. 2 Mesh distribution of VFE-2

圖3 不同截面壓強系數分布CFD 與試驗對比Fig. 3 Pressure coefficient at various sections by CFD and experiment

AGARD445.6 機翼展長為0.762 m，參考長度L取根部半弦長0.28 m，機翼詳細幾何參數見參考文獻[21]。計算網格采用六面體網格劃分，如圖4所示，網格數量為1.38×106，第一層網格厚度為10?5m 量級，動網格變形方法為擴散光順，湍流模型采用SA 模型，時間步長0.0002 s。選取來流馬赫數0.451，迎角10°，使用上述網格(最大y+為25)和細網格(最大y+為6)計算機翼90%展長處的壓強分布如圖5，同時計算翼面法向力系數(垂直翼面方向)對比如表1，可以發現兩套網格計算結果基本相同，這是由于ANSYS Fluent 軟件中對SA 模型進行了y+去敏感化處理[30]。為了節省計算量并避免出現動網格負體積問題，使用最大y+為25 的網格進行后續計算。

圖4 AGARD445.6 計算網格Fig. 4 Mesh distribution of AGARD445.6 wing

圖5 迎角10°90%展長處壓強系數分布對比Fig. 5 The comparison result of pressure coefficient at 90%span and angle of 10°

表1 迎角10°法向力系數Cn 對比Table 1 The comparison result of normal force coefficient Cn at angle of 10°

然后，采用文獻[21]的AGARD445.6 軟機翼3 號模型的顫振試驗數據驗證流固耦合方法。由于機翼顫振形式為前兩階模態耦合，計算取用前兩階模態，頻率如表2。編寫用戶自定義函數(User defined functions, UDF)完成流固耦合數值仿真。選取來流馬赫數0.499、空氣密度為0.428 kg/m3的風洞試驗條件進行0°顫振計算，通過改變大氣壓強和溫度，保證密度不變的同時改變來流速度，不同條件下求得兩階模態的廣義位移響應如圖6。最終本文計算得到顫振速度與文獻試驗結果[26]對比如表3，本文采用的流固耦合方法計算結果與試驗結果基本一致。

表3 試驗與流固耦合計算結果對比Table 3 Flutter results of CFD and experiment

圖6 0°迎角不同來流條件廣義位移響應Fig. 6 Time responses of wing tip displacement at various velocities (angle of 0°)

表2 AGARD445.6 軟機翼前兩階模態頻率Table 2 Mode frequency of AGARD445.6 weakened wing

2.2 不同初始迎角下的機翼顫振邊界規律

本節計算機翼0°～10°迎角的顫振特性。由于軟機翼在大迎角下靜變形過大不符合本文假設，使用硬機翼模型進行后續研究，前兩階模態頻率如表4[26]。

表4 AGARD445.6 硬機翼前兩階模態頻率Table 4 Mode frequency of AGARD445.6 solid wing

選取來流馬赫數0.451、空氣密度1.113 kg/m3。為了減少計算量的同時便于得到迎角與顫振特性的規律，本文設計了不同速度和初始迎角的計算狀態點，計算得到的顫振特標準差誤差棒如圖7所示。通過計算每個狀態點的廣義位移時域響應即可得到該迎角下顫振邊界和顫振頻率范圍，以10°迎角為例，通過數值仿真計算發現V*=0.459 時響應收斂、V*=0.472 時響應發散，如圖8 所示，則10°迎角的顫振邊界和顫振頻率比在兩者之間。從圖7 中不同狀態的顫振計算結果中可以看出，在迎角0°～7°，顫振邊界減小，而7°～10°初始迎角時顫振邊界反而增加。

圖7 計算狀態點和顫振特性標準差誤差棒Fig. 7 Calculation conditions and standard error bar of flutter boundary and frequency ratio

圖8 10°迎角不同來流條件廣義位移時域響應Fig. 8 Aeroelastic responses at different reduced velocities (angle of 10°)

3 基于降階模型的顫振預測與結果分析

3.1 10°迎角氣動降階與顫振預測

機翼在大迎角時，速度的不同會影響機翼的靜變形量，從而導致非定常氣動力特性的改變。為了探究ROM 中速度輸入的影響，在10°迎角條件下建立兩種氣動力ROM：ROM1 僅有廣義位移輸入，不考慮速度的輸入，輸入輸出形式如式(9)；ROM2以速度和廣義位移為輸入，輸入輸出形式如式(10)。

式中：fij為第j階廣義運動引起的第i階時域廣義氣動力系數；Δt為時間步長；n為時間步延遲數；m為樣本量。

靜氣彈基礎上設計正弦疊加形式廣義位移的訓練信號如式(11)：

式中：ηi為廣義運動幅值的控制系數；正弦疊加的幅值As為隨機數；頻率φs從5 Hz 開始，間隔10 Hz 直到65 Hz。分別給定兩階訓練信號，計算1000 個時間步，對應的廣義力系數作為ROM1 訓練集。ROM1 訓練信號圖像如圖9。

圖9 ROM1 訓練信號Fig. 9 Training signal of ROM1

ROM2 的訓練集中，選取4 個速度條件，在不同速度下，分別給定各階廣義位移以式(11)的訓練信號，計算對應的廣義力系數作為ROM2 訓練集。由于給定了4 個速度階梯，ROM2 訓練集的CFD 計算量為ROM1 的4 倍。ROM2 訓練信號和速度條件如圖10。

圖10 ROM2 訓練信號Fig. 10 Training signal of ROM2

設計不同速度、減縮頻率條件正弦形式廣義位移測試信號，并得到頻域廣義力系數矩陣AIC作為測試集。由于大迎角下存在的氣動力非線性，測試信號的正弦運動選取小幅值以減少氣動力非線性對頻域廣義力系數矩陣AIC和顫振特性預測的影響，一階廣義運動幅值取0.0003，二階廣義運動幅值取0.0002。兩種ROM 訓練集和測試集包含的減縮頻率、速度信息如圖11，測試集中，選取了4 個訓練集范圍外的減縮頻率和速度條件，測試降階模型的泛化能力。

圖11 訓練集、測試集包含的減縮頻率、速度信息Fig. 11 Ranges of reduced frequencies and reduced velocities of training data set and test data set

兩種ROM 辨識結果如圖12 和圖13，辨識誤差如表5，其中定義誤差e為式(12)的形式，其中：yCFD為CFD 計算結果向量；yROM為降階模型計算結果向量。

圖12 ROM1 辨識結果Fig. 12 Identification results of ROM1

圖13 ROM2 辨識結果Fig. 13 Identification results of ROM2

表5 ROM1 和ROM2 辨識誤差Table 5 Identification error of ROM1 and ROM2

使用得到的兩種ROM，預測速度V*=0.428，不同減縮頻率下的頻域氣動力系數矩陣AIC，并與預測集CFD 計算結果對比如圖14，速度不變時，兩種降階模型預測集結果均與CFD 計算結果基本一致，且訓練集之外部分的預測結果也較好，證明了ROM 的泛化能力；基于兩種ROM 預測減縮頻率k=0.4，不同速度下的頻域氣動力系數矩陣AIC，并與預測集CFD 計算結果對比如圖15，顯然ROM2 的預測與CFD 結果吻合更好。

圖14 V*=0.428，不同減縮頻率k 的AIC 分量Fig. 14 AIC components at V*=0.428 and various reduced frequencies

圖15 k=0.4，不同速度V*的AIC 分量Fig. 15 AIC components at k=0.4 and various reduced velocities

基于建立的兩種ROM，使用RK 法和VG 法的顫振預測如圖16～圖19，顫振預測結果對比如表6 所示，結果仍然是ROM2 的顫振邊界預測精度更高，且能較準確的預測不同速度下的廣義位移時域響應?？梢钥闯?，考慮速度輸入的ROM2的不同速度的AIC預測和顫振預測精度更高，但增加速度輸入后，ROM 訓練集需要的CFD 計算代價也將大量增加；不考慮速度輸入的ROM1 雖然不同速度的AIC預測和顫振預測誤差相對較大，但是訓練集需要的CFD 計算代價較低。另外，RK 法和VG 法預測結果基本一致。

圖16 基于ROM1，RK 法的顫振預測Fig. 16 Flutter prediction by ROM1, RK method

圖17 基于ROM2，RK 法的顫振預測Fig. 17 Flutter prediction by ROM2 and RK method

表6 基于降階模型的10°迎角顫振預測與流固耦合對比Table 6 Flutter prediction by ROM and FSI, α = 10°

圖18 基于ROM1，VG 法的顫振預測Fig. 18 Flutter prediction by ROM1 and VG method

圖19 基于ROM2，VG 法的顫振預測Fig. 19 Flutter prediction by ROM2 and VG method

3.2 考慮初始迎角輸入的氣動降階與顫振預測

為了進行不同迎角的顫振預測和顫振特性分析，需要建立考慮初始迎角α 輸入的非定常氣動力ROM。為了降低CFD 計算量同時得到定性正確的規律，降階模型ROM3 輸入為不同時間步的廣義位移和初始迎角，而不考慮速度輸入，輸入形式如式(13)。

為了對迎角7°顫振邊界拐點附近的非定常氣動力特性進行更精細的建模，ROM3 分為兩部分：ROM3A 模 型 的 迎 角 輸 入 為[0°, 2°, 4°]，ROM3B 模型的迎角輸入為[4°, 6°, 8°, 10°]。后續使 用ROM3A 完成迎 角[0°, 1°, 2°, 3°, 4°]、使用ROM3B 完成迎角[4°, 5°, 6°, 7°, 8°, 9°, 10°]的頻域氣動力系數矩陣AIC和顫振預測。廣義位移訓練信號依然采用式(11)，訓練信號和初始迎角條件如圖20 所示；測試集選擇計算減縮頻率k=0.4，迎角α=[0°, 3°, 5°, 7°, 10°]下的頻域氣動力系數矩陣AIC。訓練集、測試集計算的速度V*=0.443。降階模型的辨識誤差如表7 所示。

圖20 ROM3 訓練信號Fig. 20 Training signal of ROM3

表7 ROM3 辨識誤差Table 7 Identification error of ROM3

基于ROM3 預測速度V*=0.443，減縮頻率k=0.4，不同初始迎角下的頻域氣動力系數矩陣AIC，并與測試集CFD 計算結果對比如圖21，ROM3 能準確預測AIC隨迎角變化的規律；使用RK 法和VG 法預測得到的不同迎角顫振特性如圖22 所示，兩種預測方法的結果基本一致。雖然隨著迎角的增加，來流速度對機翼靜變形量的影響增大，但ROM3 預測的顫振邊界隨迎角變化趨勢與流固耦合結果基本一致，說明ROM3 雖然沒有考慮速度輸入，但仍能得到正確的顫振邊界變化規律。之后ROM3 將用于迎角對顫振特性影響的機理分析中。

圖21 k=0.4，V*=0.443，不同迎角α 的AIC 分量Fig. 21 AIC components at k=0.4, V*=0.443 and various angles of attack

圖22 不同迎角ROM3 顫振預測與流固耦合結果對比Fig. 22 Flutter prediction by ROM3 and CFD at various angles of attack

4 迎角對顫振邊界影響的機理分析

將圖21 中ROM3 預測的頻域氣動力系數矩陣AIC轉化為時域廣義力系數fij，其與廣義位移dj的幅值比和相位差隨迎角的變化規律如圖23 所示。發現5°迎角內，相位差基本不變，廣義力系數的幅值比隨迎角的增加略有增大；7°迎角后，隨著迎角的增加，4 個廣義力系數的相位差開始變化，而廣義力系數f21和f22的幅值比下降。即迎角小于5°時，靜氣彈基礎上的相同運動下，廣義力的幅值隨迎角增加略有增大；迎角大于7°后，靜氣彈基礎上的相同運動下，二階廣義力的幅值隨迎角增加開始下降。

圖23 k=0.4，V*=0.443，不同迎角α 下的廣義力系數fij 與廣義位移dj 的幅值比和相位差Fig. 23 Amplitude ratio and phase between general force coefficient and general displacement at k=0.4, V*=0.443 and various angles of attack by ROM3

為了探究上述規律的原因，選取來流馬赫數0.42，迎 角0°、5°、7°、10°，速 度V*=0.427 的狀態進行研究。后三個非零迎角條件下，靜氣彈收斂后的機翼上表面壓強云圖和流線圖如圖24 所示?？梢园l現5°迎角前機翼表面基本為附著流；隨著迎角的增加，展向渦導致機翼上表面的分離區范圍擴大，改變了上翼面的壓強分布。然后在靜氣彈基礎上分別給予兩階模態相同幅值、減縮頻率k=0.42 的正弦運動。得到的廣義力系數fij可以表示為式(14)：

圖24 機翼上表面壓強云圖和流線圖Fig. 24 Pressure contour and streamline of the upper surface of the wing

式中：ψi為第i階模態的z方向分量；S為面積；δCp為相對于靜氣彈的上下翼面壓強系數之差。在分別給兩階模態正弦運動后，選取廣義力系數正弦的第一個峰值時刻，對比展向90%截面的δCp弦向分布如圖25 所示?？梢钥闯?，小迎角時δCp主要集中于前緣，隨著迎角的增加運動引起的δCp逐漸后移；大迎角后展向渦出現，導致δCp逐漸平均到整個翼弦。

圖25 不同迎角展向90%截面的δCp 弦向分布Fig. 25 δCp distribution at 90% span section for different angles

δCpψi體現了第i階廣義力系數在翼面上的分布，其于90%展長截面的弦向分布如圖26 所示。可以發現，5°迎角前廣義力系數分布與δCp分布類似，隨著迎角的增加，其分布逐漸后移；5°迎角后，由于彎曲模態ψ1為正值，一階廣義力系數分布皆為正，而扭轉模態ψ2沿弦向從正到負，導致二階廣義力系數分布出現負值。經過式(14)積分后得到該截面不同迎角的廣義力系數如圖27 所示。隨著迎角的增加，小迎角時截面廣義力系數略有增加，大迎角時截面二階廣義力系數f21和f22下降。

圖26 90%展長截面廣義力系數分布Fig. 26 Generalized force coefficients at 90% span section

圖27 通過式(14)積分得到90%展長截面的廣義力系數Fig. 27 Generalized force coefficients at 90% span section is obtained by integrating equation (14)

5 結論

本文發展了流固耦合方法和ROM 方法計算初始迎角對顫振邊界的影響。針對AGARD445.6 硬機翼模型，計算不同迎角和速度下的氣動彈性響應并總結出了顫振特性隨迎角的變化規律；在10°迎角下研究了考慮速度輸入和不考慮速度的非定常氣動力ROM，并預測了10°迎角的顫振特性；建立了考慮初始迎角輸入的非定常氣動力ROM 并預測了0°～10°不同迎角的顫振特性，最后分析了顫振邊界隨迎角變化的機理。本文主要結論如下：

(1)對于本文的AGARD445.6 硬機翼，當初始迎角小于5°時，隨著迎角的增大，廣義力系數幅值略微增加，導致顫振邊界逐漸下降；初始迎角7°～10°，隨著迎角的增大，展向渦的產生改變了機翼上翼面壓強分布，導致二階廣義力系數的幅值下降，最終導致顫振邊界的增加。

(2)基于RBFNN 的非定常氣動力ROM 可以準確預測不同速度、減縮頻率的非定常氣動力，預測10°迎角的顫振邊界時，考慮速度輸入后ROM 的預測結果更精確，另外顫振預測的時域RK 法和頻域VG 法預測結果基本一致。

(3)不考慮速度輸入而考慮迎角輸入的ROM能正確預測的顫振邊界隨迎角變化的趨勢，與流固耦合計算不同迎角的顫振特性規律基本一致。