周 萌，高國柱

(中國電子科技集團(tuán)公司第三十八研究所 浮空平臺部，合肥 230088)

0 引言

為了解決CFD/CSD計算效率低下的問題，在20世紀(jì)90年代，杜克大學(xué)的Dowell[1]和美國國家航天航空局的Silva[2]等人將降階模型(ROM)方法應(yīng)用到氣動彈性領(lǐng)域中。構(gòu)造降階模型的目的主要有兩點(diǎn)[3]：一是降階模型技術(shù)大大減小原數(shù)值模型系統(tǒng)的階數(shù)，且在保證較為精確的系統(tǒng)動力學(xué)特征描述的同時大大減少計算耗時；二是降階模型為學(xué)者們研究系統(tǒng)動力學(xué)特征提供了強(qiáng)有力的工具。目前，降階模型主要分為兩類：一類是對流場特征結(jié)構(gòu)進(jìn)行分析的降階模型，其中以基于本征正交分解(Proper Orthogonal Decomposition，POD)技術(shù)[4]的降階模型和基于諧波平衡(Harmonic Balance，HB)方法[5]的降階模型為代表；另一類是基于系統(tǒng)外部輸入和輸出響應(yīng)的系統(tǒng)辨識方法的降階模型，這種方法以Silver的Volterra模型[6]、基于代理模型[7]和基于時間序列的非線性帶外部輸入的自回歸滑動平均模型[8](NARMAX)的降階模型為代表。HB方法是一種在頻域內(nèi)求解非線性周期流場的方法。POD方法是一種獲得高維空間在低維空間內(nèi)近似表達(dá)的方法，通過CFD計算快照構(gòu)造相關(guān)矩陣，獲得特征模態(tài)形成降階子空間，將控制方程投影到系統(tǒng)特征來建立ROM。HB方法和POD方法是在全階系統(tǒng)的基礎(chǔ)上，建立新的控制方程并對原系統(tǒng)數(shù)值求解源代碼進(jìn)行修改，進(jìn)而提高分析效率。系統(tǒng)辨識方法相較于HB方法和POD方法處理方式簡單，適用性強(qiáng)，能夠適用于線性和非線性問題。

目前，國內(nèi)外研究人員基于不同的ROM構(gòu)造了多種分析方法，并與標(biāo)準(zhǔn)算例進(jìn)行了對比。雖然針對計算條件變化(如改變初始迎角[9]、CFD模型[10]、粘性[11]和結(jié)構(gòu)參數(shù)[12]等)對二維翼型顫振邊界的影響開展了一定的研究，但是在跨聲速內(nèi)，翼型厚度變化對顫振邊界的影響則尚未給出研究規(guī)律。

本文基于CFD技術(shù)，采用系統(tǒng)辨識方法中的ARMA模型構(gòu)造非定常氣動力降階模型[7]，耦合結(jié)構(gòu)運(yùn)動方程，建立一套頻域/時域氣動彈性ROM分析系統(tǒng)。然后采用該方法詳細(xì)研究在翼型最大厚度所在位置保持不變時翼型厚度對顫振邊界的影響。

1 氣動彈性系統(tǒng)降階模型

1.1 氣動力降階模型

基于ARMA模型的非定常氣動力降階模型是一種基于時間序列的降階模型[3]，可以表述為如下方程：

其中，fa(k)為第k步得到的廣義氣動力；na、nb分別為系統(tǒng)輸出和輸入的延遲階數(shù)；q表示輸入的廣義位移；Ai、Bi是模型中待辨識的參數(shù)。采用基于過濾高斯白噪聲(FWGN)的模態(tài)激勵信號進(jìn)行模型辨識，該信號具有頻帶寬且幅值范圍廣的特點(diǎn)。

在構(gòu)建氣動力降階模型時，采用MIMO (Multiple Input Multiple Output，MIMO) 系統(tǒng)方法，即一次輸入n個位移信號同時激勵n個模態(tài)的方式，通過最小二乘法擬合CFD求解出來的曲線，辨識出Ai、Bi矩陣。

定義狀態(tài)向量xa(k)如下：

將離散差分方程轉(zhuǎn)化為離散狀態(tài)空間方程形式，則氣動力降階模型的狀態(tài)空間方程和輸出方程可以表示為：

其中，

1.2 氣動彈性降階模型

氣動彈性系統(tǒng)的結(jié)構(gòu)運(yùn)動方程對應(yīng)的狀態(tài)空間方程可以表達(dá)為如下形式：

其中，

式中，q表示位移，M表示質(zhì)量矩陣，K表示剛度矩陣，D表示阻尼矩陣。

將離散狀態(tài)的氣動力降階模型轉(zhuǎn)化為連續(xù)時間狀態(tài)空間方程，并忽略靜平衡氣動力，即

式中，各字母含義與式(2)中一致，表示在連續(xù)時間狀態(tài)下的氣動系統(tǒng)矩陣系數(shù)。

耦合結(jié)構(gòu)運(yùn)動方程和氣動力降階模型方程，得到氣動彈性系統(tǒng)的狀態(tài)空間形式的控制方程：

其中，Q表示動壓，狀態(tài)空間矩陣如下：

可以對狀態(tài)空間矩陣求解不同動壓下的特征值，并采用模態(tài)跟蹤技術(shù)[13]畫出根軌跡，根據(jù)根軌跡分析系統(tǒng)的穩(wěn)定性。

將上式轉(zhuǎn)化為狀態(tài)空間方程的形式通過時間推進(jìn)求解[14]，如下：

其中，

采用葉正寅和張偉偉等人[15]構(gòu)建的四階隱式Adams線性多步法，通過預(yù)估-校正技術(shù)求解該隱式方程，提高氣動彈性ROM的分析效率和魯棒性[14]。

綜上，建立基于系統(tǒng)辨識技術(shù)的氣動彈性時域分析系統(tǒng)的耦合流程具體步驟如圖1所示。

1.3 模型驗證

為驗證本文氣動彈性ROM的精度，選取氣動彈性分析的標(biāo)準(zhǔn)算例——Isogai案例A模型，該翼型由于在跨聲速狀態(tài)下呈現(xiàn)復(fù)雜的穩(wěn)定性特性而被廣泛應(yīng)用于氣動彈性程序驗證。Isogai模型的翼型為NACA64A010翼型，是Isogai[16]提出計算顫振的二維翼型。計算中CFD求解器選取k-ω SST湍流模型，按照文獻(xiàn)[17]，在所有馬赫數(shù)下，基于弦長的雷諾數(shù)均給定為6×106。

在無粘條件下計算結(jié)果呈現(xiàn)為“S”形顫振邊界，與參考文獻(xiàn)結(jié)果一致。粘性顫振邊界數(shù)值結(jié)果表明，跨聲速凹坑最低點(diǎn)對應(yīng)的馬赫數(shù)大約為0.83，隨后無因次顫振速度急劇增大；在Ma=0.87時無因次顫振速度開始減小，在Ma=0.9時達(dá)到第二個凹坑最低點(diǎn)，隨后隨著馬赫數(shù)增加無因次顫振速度繼續(xù)增大；當(dāng)Ma≥0.93時，無因次顫振速度隨馬赫數(shù)的增加變化很小，基本保持在3.8左右。

圖1基于系統(tǒng)辨識技術(shù)的氣動彈性時域分析系統(tǒng)的耦合流程

上述結(jié)果是采用Intel(R) Xeon(R) CPU E5-2680 @2.50 GHz共19個進(jìn)程并行計算，表1給出了直接采用CFD/CSD方法和基于ROM的頻域分析方法的計算時間對比?？梢缘玫?，直接采用CFD/CSD方法計算耗時約7200 min，采用ROM預(yù)測顫振邊界耗時246.4 min，計算效率提高1～2個數(shù)量級。

(a)無粘條件

(b)粘性條件

圖2 采用CFD和ROM得到的Isogai翼型無因次顫振速度結(jié)果對比變化趨勢

2 翼型厚度對顫振邊界影響

翼型厚度影響著翼型的氣動性能，是飛行器精細(xì)化設(shè)計時的重要參數(shù)，因此有必要研究翼型外形參數(shù)對顫振特性的影響。本小節(jié)研究了在最大厚度位置保持不變時，翼型厚度對顫振邊界的影響。選取的翼型分別為NACA64A008、NACA64A009、NACA64A010、NACA64A011和NACA64A012，采用粘性條件，計算所用的氣動參數(shù)和結(jié)構(gòu)參數(shù)與1.3小節(jié)中的保持一致。

翼型厚度對顫振邊界的影響如圖3所示，圖中展示了在最大厚度位置保持不變時，不同厚度翼型的無因次顫振速度邊界和頻率比邊界。不同表面壓力分布和摩阻系數(shù)分布對比圖如圖4～8所示。

數(shù)值結(jié)果表明：

(1)當(dāng)馬赫數(shù)Ma<0.80時，隨著翼型厚度增加，顫振速度略有減小，但減小幅度在10 %以內(nèi)。從圖4可以看出，此時由于厚度的增加，導(dǎo)致激波強(qiáng)度逐漸增加，進(jìn)而使得顫振速度逐漸減小。

(2)當(dāng)馬赫數(shù)Ma>0.91時，隨著翼型厚度增加，顫振速度逐漸增加；且同一翼型隨著馬赫數(shù)的增加，顫振速度基本保持不變。從圖8可以看出隨著厚度的增加，分離區(qū)逐漸增加，激波強(qiáng)度基本保持不變，因此顫振速度逐漸增加。

(3)跨聲速凹坑隨著翼型厚度的增加而逐漸左移。從圖5～7可以看出，在分離區(qū)作用下，跨聲速凹坑左移，也就是第一次顫振速度增加。當(dāng)進(jìn)一步增加馬赫數(shù)時，分離區(qū)域達(dá)到翼型后緣且隨著馬赫數(shù)增加而分離區(qū)減小，此時在激波的作用下，無因次顫振速度邊界出現(xiàn)第二個臺階。

3 結(jié)語

本文采用ARMA模型構(gòu)建了非定常氣動力ROM，耦合結(jié)構(gòu)運(yùn)動方程，建立了頻域/時域的氣動彈性ROM。采用該方法研究Isogai翼型的跨聲速顫振問題，并與CFD/CSD計算結(jié)果對比，結(jié)果表明該方法求解結(jié)果可靠。研究了最大厚度位置保持不變時，翼型厚度對顫振邊界的影響，數(shù)值結(jié)果表明：當(dāng)Ma<0.80時，隨著翼型厚度增加，顫振速度略有減??；當(dāng)馬赫數(shù)Ma>0.91時，隨著翼型的厚度增加，顫振速度逐漸增加?？缏曀侔伎与S著翼型厚度的增加而逐漸左移。

綜上所述，采用ROM可以精確快速地預(yù)測顫振邊界。當(dāng)翼型最大厚度所在位置保持不變時，為了達(dá)到提高顫振速度的目標(biāo)，可以嘗試采取改變機(jī)翼翼型厚度的方式提高機(jī)翼對飛行環(huán)境的適應(yīng)能力。