基于降階模型的不同厚度翼型顫振邊界預(yù)測(cè)

周 萌，高國(guó)柱

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

0 引言

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

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

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

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

1.1 氣動(dòng)力降階模型

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

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

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

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

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

其中，

1.2 氣動(dòng)彈性降階模型

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

其中，

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

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

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

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

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

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

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

其中，

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

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

1.3 模型驗(yàn)證

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

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

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

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

(a)無(wú)粘條件

(b)粘性條件

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

2 翼型厚度對(duì)顫振邊界影響

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

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

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

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

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

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

3 結(jié)語(yǔ)

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

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