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

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

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

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

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

1 計(jì)算方法

1.1 基于模態(tài)疊加的流固耦合方法

多自由度線性系統(tǒng)受迫振動(dòng)方程為：

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

式中：Φ 為結(jié)構(gòu)質(zhì)量歸一模態(tài)；d為廣義坐標(biāo)，則式(1)可化為：

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

在流固耦合計(jì)算中引入假設(shè)：只考慮垂直于機(jī)翼所處平面方向的變形，即機(jī)翼處于xoy平面，模態(tài)只考慮z方向分量。該假設(shè)在靜變形和動(dòng)態(tài)變形不太大時(shí)是合理的。

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

1)給系統(tǒng)施加阻尼，計(jì)算式(3)直至收斂得到靜氣彈結(jié)果；

2)在靜氣彈的基礎(chǔ)上施加廣義速度擾動(dòng)d˙0，同時(shí)取消阻尼，計(jì)算式(3)得到動(dòng)氣彈時(shí)域響應(yīng)。

RBFNN 由BROOMHEAD 等[27]提出，其結(jié)構(gòu)包含三層：輸入層、隱含層和輸出層，其中隱含層的激活函數(shù)為徑向基函數(shù)，輸入靠近基函數(shù)中心點(diǎn)時(shí)神經(jīng)元才會(huì)被激活，輸出層激活函數(shù)為線性函數(shù)。典型的RBFNN 結(jié)構(gòu)如圖1 所示。

圖1 徑向基神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig. 1 Structure of RBFNN

1.2 基于徑向基神經(jīng)網(wǎng)絡(luò)的非定常氣動(dòng)降階方法

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

式中：yc為輸出層第c個(gè)輸出；wcv為隱含層第v個(gè)神經(jīng)元到輸出層第c個(gè)神經(jīng)元的權(quán)重。RBFNN 訓(xùn)練過程中，在隱含層神經(jīng)元個(gè)數(shù)、基函數(shù)、中心向量和寬度確定后，權(quán)重的學(xué)習(xí)僅需求解線性方程，因而RBFNN 相比于傳統(tǒng)多層前饋神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)簡單、訓(xùn)練成本低、收斂速度快。另外，只要隱含層數(shù)量足夠大，RBFNN 能以任意精度逼近非線性函數(shù)。

RBFNN 一般有兩種：正則RBFNN 和廣義RBFNN。正則RBFNN 的隱含層數(shù)量與訓(xùn)練樣本量相同，隱含層中心向量與對(duì)應(yīng)的輸入向量相同，正則RBFNN 通常足夠精確，但考慮本文的問題，輸入樣本量較大導(dǎo)致神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)龐大，可能導(dǎo)致訓(xùn)練困難、矩陣出現(xiàn)病態(tài)，訓(xùn)練結(jié)果可能過擬合；本文選擇建立的廣義RBFNN，則是逐個(gè)添加隱含層神經(jīng)元，直到均方誤差達(dá)到設(shè)定要求，建立的模型則為滿足誤差要求的最簡RBFNN 結(jié)構(gòu)，詳細(xì)步驟為：

1)用神經(jīng)網(wǎng)絡(luò)擬合訓(xùn)練樣本，得到訓(xùn)練樣本中誤差最大的輸入向量；

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

3)優(yōu)化權(quán)重矩陣，返回步驟1)，直到均方誤差達(dá)到設(shè)定要求或隱含層神經(jīng)元數(shù)量達(dá)到設(shè)定上限。

1.3 基于氣動(dòng)降階模型的顫振預(yù)測方法

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

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

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

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

式中：L為參考長度；ρ 為密度；V為來流速度；ωα為扭轉(zhuǎn)模態(tài)角頻率；μ為質(zhì)量系數(shù)；AIC為頻域氣動(dòng)力系數(shù)矩陣；gs為阻尼；ω為結(jié)構(gòu)運(yùn)動(dòng)角頻率。AIC與減縮頻率和來流速度有關(guān)，因?yàn)闄C(jī)翼在大迎角時(shí)速度的改變會(huì)改變機(jī)翼的靜變形位置，從而對(duì)非定常氣動(dòng)特性產(chǎn)生影響，進(jìn)而影響AIC矩陣。通過降階模型預(yù)測某一減縮頻率、速度下的時(shí)域廣義力系數(shù)，再將結(jié)果轉(zhuǎn)化到頻域，即可得到該減縮頻率、速度下的AIC矩陣?；诮惦A模型，使用VG 法預(yù)測某一迎角顫振特性的步驟為：

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

2 基于流固耦合的數(shù)值算例與結(jié)果分析

2.1 流固耦合方法驗(yàn)證

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

圖2 VFE-2 計(jì)算網(wǎng)格Fig. 2 Mesh distribution of VFE-2

圖3 不同截面壓強(qiáng)系數(shù)分布CFD 與試驗(yàn)對(duì)比Fig. 3 Pressure coefficient at various sections by CFD and experiment

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

圖4 AGARD445.6 計(jì)算網(wǎng)格Fig. 4 Mesh distribution of AGARD445.6 wing

圖5 迎角10°90%展長處壓強(qiáng)系數(shù)分布對(duì)比Fig. 5 The comparison result of pressure coefficient at 90%span and angle of 10°

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

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

表3 試驗(yàn)與流固耦合計(jì)算結(jié)果對(duì)比Table 3 Flutter results of CFD and experiment

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

表2 AGARD445.6 軟機(jī)翼前兩階模態(tài)頻率Table 2 Mode frequency of AGARD445.6 weakened wing

2.2 不同初始迎角下的機(jī)翼顫振邊界規(guī)律

本節(jié)計(jì)算機(jī)翼0°～10°迎角的顫振特性。由于軟機(jī)翼在大迎角下靜變形過大不符合本文假設(shè)，使用硬機(jī)翼模型進(jìn)行后續(xù)研究，前兩階模態(tài)頻率如表4[26]。

表4 AGARD445.6 硬機(jī)翼前兩階模態(tài)頻率Table 4 Mode frequency of AGARD445.6 solid wing

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

圖7 計(jì)算狀態(tài)點(diǎn)和顫振特性標(biāo)準(zhǔn)差誤差棒Fig. 7 Calculation conditions and standard error bar of flutter boundary and frequency ratio

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

3 基于降階模型的顫振預(yù)測與結(jié)果分析

3.1 10°迎角氣動(dòng)降階與顫振預(yù)測

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

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

靜氣彈基礎(chǔ)上設(shè)計(jì)正弦疊加形式廣義位移的訓(xùn)練信號(hào)如式(11)：

式中：ηi為廣義運(yùn)動(dòng)幅值的控制系數(shù)；正弦疊加的幅值A(chǔ)s為隨機(jī)數(shù)；頻率φs從5 Hz 開始，間隔10 Hz 直到65 Hz。分別給定兩階訓(xùn)練信號(hào)，計(jì)算1000 個(gè)時(shí)間步，對(duì)應(yīng)的廣義力系數(shù)作為ROM1 訓(xùn)練集。ROM1 訓(xùn)練信號(hào)圖像如圖9。

圖9 ROM1 訓(xùn)練信號(hào)Fig. 9 Training signal of ROM1

ROM2 的訓(xùn)練集中，選取4 個(gè)速度條件，在不同速度下，分別給定各階廣義位移以式(11)的訓(xùn)練信號(hào)，計(jì)算對(duì)應(yīng)的廣義力系數(shù)作為ROM2 訓(xùn)練集。由于給定了4 個(gè)速度階梯，ROM2 訓(xùn)練集的CFD 計(jì)算量為ROM1 的4 倍。ROM2 訓(xùn)練信號(hào)和速度條件如圖10。

圖10 ROM2 訓(xùn)練信號(hào)Fig. 10 Training signal of ROM2

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

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

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

圖12 ROM1 辨識(shí)結(jié)果Fig. 12 Identification results of ROM1

圖13 ROM2 辨識(shí)結(jié)果Fig. 13 Identification results of ROM2

表5 ROM1 和ROM2 辨識(shí)誤差Table 5 Identification error of ROM1 and ROM2

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

圖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 法的顫振預(yù)測如圖16～圖19，顫振預(yù)測結(jié)果對(duì)比如表6 所示，結(jié)果仍然是ROM2 的顫振邊界預(yù)測精度更高，且能較準(zhǔn)確的預(yù)測不同速度下的廣義位移時(shí)域響應(yīng)?？梢钥闯?，考慮速度輸入的ROM2的不同速度的AIC預(yù)測和顫振預(yù)測精度更高，但增加速度輸入后，ROM 訓(xùn)練集需要的CFD 計(jì)算代價(jià)也將大量增加；不考慮速度輸入的ROM1 雖然不同速度的AIC預(yù)測和顫振預(yù)測誤差相對(duì)較大，但是訓(xùn)練集需要的CFD 計(jì)算代價(jià)較低。另外，RK 法和VG 法預(yù)測結(jié)果基本一致。

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

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

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

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

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

3.2 考慮初始迎角輸入的氣動(dòng)降階與顫振預(yù)測

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

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

圖20 ROM3 訓(xùn)練信號(hào)Fig. 20 Training signal of ROM3

表7 ROM3 辨識(shí)誤差Table 7 Identification error of ROM3

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

圖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 顫振預(yù)測與流固耦合結(jié)果對(duì)比Fig. 22 Flutter prediction by ROM3 and CFD at various angles of attack

4 迎角對(duì)顫振邊界影響的機(jī)理分析

將圖21 中ROM3 預(yù)測的頻域氣動(dòng)力系數(shù)矩陣AIC轉(zhuǎn)化為時(shí)域廣義力系數(shù)fij，其與廣義位移dj的幅值比和相位差隨迎角的變化規(guī)律如圖23 所示。發(fā)現(xiàn)5°迎角內(nèi)，相位差基本不變，廣義力系數(shù)的幅值比隨迎角的增加略有增大；7°迎角后，隨著迎角的增加，4 個(gè)廣義力系數(shù)的相位差開始變化，而廣義力系數(shù)f21和f22的幅值比下降。即迎角小于5°時(shí)，靜氣彈基礎(chǔ)上的相同運(yùn)動(dòng)下，廣義力的幅值隨迎角增加略有增大；迎角大于7°后，靜氣彈基礎(chǔ)上的相同運(yùn)動(dòng)下，二階廣義力的幅值隨迎角增加開始下降。

圖23 k=0.4，V*=0.443，不同迎角α 下的廣義力系數(shù)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

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

圖24 機(jī)翼上表面壓強(qiáng)云圖和流線圖Fig. 24 Pressure contour and streamline of the upper surface of the wing

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

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

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

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

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

5 結(jié)論

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

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

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

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