基于POD的尾流激勵(lì)葉片氣動(dòng)力降階模型

劉 釗，王沐晨，李金玖，張 珺，李立州

(1.中北大學(xué) 機(jī)電工程學(xué)院， 太原 030051； 2.太原學(xué)院數(shù)學(xué)系， 太原 030001)

保證航空發(fā)動(dòng)機(jī)具有較長(zhǎng)使用壽命和較強(qiáng)可靠性，一直是現(xiàn)代航空發(fā)動(dòng)機(jī)設(shè)計(jì)過程中追求的目標(biāo)[1]。在發(fā)動(dòng)機(jī)內(nèi)部動(dòng)、靜葉之間的高速相對(duì)轉(zhuǎn)動(dòng)，使得下游葉片極容易受到周期性的上游尾流激勵(lì)，誘發(fā)下游葉片流場(chǎng)產(chǎn)生周期性振動(dòng)，進(jìn)而使得下游葉片表面氣動(dòng)力具有周期性。這種周期性的氣動(dòng)力有可能導(dǎo)致葉片強(qiáng)迫振動(dòng)，甚至出現(xiàn)嚴(yán)重的疲勞破壞，影響整個(gè)發(fā)動(dòng)機(jī)的可靠性和使用壽命。因此，研究尾流激勵(lì)下葉片的氣動(dòng)性能，能為航空發(fā)動(dòng)機(jī)葉片的氣動(dòng)優(yōu)化提供依據(jù)，對(duì)提高整個(gè)航空發(fā)動(dòng)機(jī)的質(zhì)量、使用壽命和可靠性具有重要意義。

近年來，計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)被廣泛用于研究葉片氣動(dòng)特性[3-9]，然而在工程實(shí)踐中明顯可以看到，CFD技術(shù)雖然計(jì)算精度高，但其最大的缺陷在于計(jì)算時(shí)間長(zhǎng)，當(dāng)研究需要反復(fù)迭代的問題時(shí)，其計(jì)算量更為巨大，工程應(yīng)用很不便。為有效地避免這個(gè)問題，Dowell[10]、Silva[11]等提出了基于CFD技術(shù)的非定常氣動(dòng)力降階模型(Reduced Order Model，ROM)。ROM作為一個(gè)可以表述非定常流體系統(tǒng)絕大部分動(dòng)力特性的簡(jiǎn)化數(shù)學(xué)模型[12]，其計(jì)算所消耗的內(nèi)存、時(shí)間和設(shè)計(jì)周期相對(duì)于傳統(tǒng)的CFD方法減少了十幾個(gè)甚至幾十個(gè)數(shù)量級(jí)[13]。現(xiàn)有的降階模型方法有Krylov子空間法、POD法、Volterra級(jí)數(shù)和諧波平衡(Harmonic balance，HB)等幾種[14]，其中POD法計(jì)算效率和精度較高，在工程上應(yīng)用范圍較廣。在流體力學(xué)領(lǐng)域，Lumely[15]引入POD方法用來分析紊流場(chǎng)的流動(dòng)特性，其理論方法是在尋找最小二乘意義下的最優(yōu)有序正交基。周強(qiáng)等[16]采用POD方法建立了降階的氣動(dòng)彈性系統(tǒng)狀態(tài)空間方程，實(shí)現(xiàn)了快速預(yù)測(cè)氣彈系統(tǒng)顫振邊界。Jian Z等[17]用POD法求解了曲面組合板非線性顫振的運(yùn)動(dòng)方程，與CFD/CSD結(jié)果吻合良好，節(jié)省大量計(jì)算時(shí)間。

為快速計(jì)算上游尾流激勵(lì)下葉片氣動(dòng)力，采用POD法提取主要降階基模態(tài)，構(gòu)建基于POD的尾流激勵(lì)下葉片氣動(dòng)力降階模型，并對(duì)葉片氣動(dòng)力數(shù)據(jù)重構(gòu)，實(shí)現(xiàn)快速分析上游尾流激勵(lì)下的葉片氣動(dòng)力。

1 POD降階基模型建立

(1)

為簡(jiǎn)化描述，采用矩陣u來表示流場(chǎng)氣動(dòng)力：

(2)

根據(jù)POD理論，可將尾流激勵(lì)下的葉片氣動(dòng)力分解為平均值和脈動(dòng)值的疊加[18]：

(3)

(4)

其中，αr為第r個(gè)POD基向量對(duì)應(yīng)的模態(tài)系數(shù)；φr為POD降階基向量；R表示POD基的個(gè)數(shù)，選擇合適的階次R就可以得到尾流激勵(lì)的葉片氣動(dòng)力降階模型。

圖1 葉片氣動(dòng)力系統(tǒng)示意圖

1.1 POD降階基提取

本征正交分解(POD)的原理就是將一個(gè)高維向量，通過一組最優(yōu)正交基投影到一個(gè)低階的向量空間中，保留主要特征，并且可以基本重構(gòu)出所對(duì)應(yīng)的原始高階向量。令ui，i=1，2，3，…，N是在N個(gè)瞬時(shí)提取的流場(chǎng)氣動(dòng)力數(shù)據(jù)(稱為“快照”)，可以通過數(shù)值模擬或者實(shí)驗(yàn)測(cè)量得到。POD的目標(biāo)在于找到一個(gè)正交基函數(shù)，使得：

(5)

(6)

令U={u1，u2，u3，…，uN}，對(duì)式(6)求偏導(dǎo)數(shù)有：

(7)

令Z=UUT，式(5)求最大值的問題就轉(zhuǎn)化為求POD核函數(shù)Z特征值和特征向量的問題：

Zφ=λφ

(8)

對(duì)式(8)使用奇異值分解(Singular Value Decomposition，SVD)方法，即可求出Z的特征向量以及特征值λi，i=1，2，3，…，N。定義能量百分比：

(9)

式中I(r)表示前R個(gè)特征值占總能量的份額，I(r)越接近1，表示特征向量包含的信息越完整，通常取I(r)=0.99，求得的前R個(gè)特征值對(duì)應(yīng)的特征向量就是最終的POD基φr(r=1，2，3，…，R)。

1.2 氣動(dòng)力數(shù)據(jù)重構(gòu)

對(duì)于任一時(shí)刻的尾流輸入總壓Pr，響應(yīng)的POD模態(tài)系數(shù)為αr，則有：

(10)

2 算例分析

二維葉片流場(chǎng)CFD計(jì)算模型如圖2所示，計(jì)算網(wǎng)格采用ICEM輸出的非結(jié)構(gòu)網(wǎng)格，網(wǎng)格單元數(shù)為164 606，節(jié)點(diǎn)數(shù)為92 040，葉片網(wǎng)格結(jié)構(gòu)如圖3所示，導(dǎo)入Fluent計(jì)算葉片氣動(dòng)力，流體定義為理想氣體，采用Spallart-Allmaras湍流模型，進(jìn)口設(shè)置為壓力進(jìn)口pressure-inlet。穩(wěn)態(tài)時(shí)，進(jìn)口總壓為 120 500 Pa，出口壓力101 325 Pa，溫度300 K。瞬態(tài)計(jì)算時(shí)，將UDF編譯的上游尾流加載在流場(chǎng)進(jìn)口最下方點(diǎn)，尾流速度大小10 m/s，方向?yàn)檫M(jìn)口y方向。

圖2 葉片CFD計(jì)算模型示意圖

圖3 葉片網(wǎng)格結(jié)構(gòu)示意圖

上游尾流輸入的進(jìn)口總壓波形如圖4所示，設(shè)置計(jì)算時(shí)間t=[0，0.08]，時(shí)間步長(zhǎng)為2e-5 s，進(jìn)口壓力隨時(shí)間呈現(xiàn)周期性變化，計(jì)算可得4 000個(gè)進(jìn)口總壓的葉片氣動(dòng)力數(shù)據(jù)并保存。

圖4 周期信號(hào)壓力波波形

將4 000個(gè)原氣動(dòng)力系統(tǒng)的總壓數(shù)據(jù)劃分為100個(gè)具有相同概率的非重疊區(qū)間，從每個(gè)區(qū)間中按均勻分布隨機(jī)地抽取一個(gè)點(diǎn)，得到100個(gè)時(shí)刻的葉片表面氣動(dòng)力，作為氣動(dòng)力快照矩陣，對(duì)快照矩陣做POD分解，得到100個(gè)POD基模態(tài)。圖5分別展示了樣本點(diǎn)計(jì)算的CFD結(jié)果與已知CFD結(jié)果的絕對(duì)誤差，誤差范圍分別為[0.01，0.027]與[-0.05，-0.027 5]，均在可接受范圍內(nèi)，說明樣本點(diǎn)CFD結(jié)果的準(zhǔn)確度高，可以忽略樣本點(diǎn)CFD計(jì)算數(shù)據(jù)對(duì)研究結(jié)果的影響。

圖5 CFD計(jì)算結(jié)果絕對(duì)誤差曲線

圖6 能量比值隨POD模態(tài)變化曲線

降階模型計(jì)算結(jié)果與CFD計(jì)算結(jié)果如圖7所示。圖中，紅色實(shí)線為CFD計(jì)算結(jié)果，藍(lán)色實(shí)線為POD降階模型的計(jì)算結(jié)果。從圖中可以看出，降階模型計(jì)算結(jié)果與CFD計(jì)算結(jié)果一致，因此可以采用POD方法構(gòu)建尾流激勵(lì)氣動(dòng)力降階模型。同時(shí)CFD計(jì)算結(jié)果與降階模型計(jì)算結(jié)果之間存在一定的誤差。

為了評(píng)估降階模型精度，引入了均方根誤差RMSE和平均絕對(duì)誤差MAE，其基本公式分別為：

(11)

(12)

式中：i為第i個(gè)計(jì)算結(jié)果的氣動(dòng)力數(shù)據(jù)，Oi為POD降階模型得到的氣動(dòng)力數(shù)據(jù)，F(xiàn)i為CFD求解的氣動(dòng)力數(shù)據(jù)。分別將數(shù)據(jù)代入式(11)和式(12)可以得到系統(tǒng)的均方根誤差和平均絕對(duì)誤差，如表1。

圖7 上游尾流激勵(lì)下的葉片氣動(dòng)力響應(yīng)隨時(shí)間的變化曲線

表1 周期信號(hào)輸入的均方根誤差(RMSE)和平均絕對(duì)誤差(MAE)

由表1可知：誤差均在可接受范圍內(nèi)，故葉片氣動(dòng)力降階模型具有很高精度。

3 驗(yàn)證POD降階模型

為了進(jìn)一步說明降階模型可以計(jì)算不同尾流激勵(lì)下的葉片氣動(dòng)力，將圖4中的將尾流壓力波換成一組正弦信號(hào)波，波形如圖8所示，葉片計(jì)算的穩(wěn)態(tài)邊界條件保持不變，將正弦尾流壓力信號(hào)輸入到下游葉片的入口處，完成CFD計(jì)算，得到正弦信號(hào)激勵(lì)下葉片氣動(dòng)力數(shù)據(jù)。

將壓力信號(hào)數(shù)據(jù)直接輸入到已經(jīng)建立的葉片氣動(dòng)力降階模型，得到葉片氣動(dòng)力響應(yīng)值并保存，所得計(jì)算結(jié)果與CFD計(jì)算結(jié)果如圖9所示。利用POD計(jì)算結(jié)果和CFD計(jì)算結(jié)果，得到降階模型的均方根誤差和平均絕對(duì)誤差，如表2。

圖8 正弦信號(hào)壓力波波形圖

圖9 游尾流激勵(lì)下的葉片氣動(dòng)力響應(yīng)隨時(shí)間的變化曲線

表1 正弦信號(hào)輸入的均方根誤差和平均絕對(duì)誤差

對(duì)于本文第二次輸入的正弦信號(hào)，圖9(a)、圖9(b)分別為POD降階模型和CFD輸出的阻力和升力的對(duì)比圖。從圖中可以得出，POD降階模型與CFD求解結(jié)果一致。誤差數(shù)值基本在0左右波動(dòng)，進(jìn)一步說明降階模型的準(zhǔn)確性。

對(duì)比表1和表2數(shù)據(jù)可知，兩組輸入狀態(tài)下，降階模型精度相差10倍，差異較大的原因是輸入尾流的振幅不同。由于葉片氣動(dòng)力系統(tǒng)為弱非線性系統(tǒng)，算例中選擇的穩(wěn)態(tài)壓力值為120 500 Pa，周期信號(hào)的振幅為[-5 000 Pa，1 000 Pa]，正弦信號(hào)的振幅為[-175 Pa，175 Pa]，振幅越小弱非線性特征越明顯，精度越高。綜合表1和表2的數(shù)據(jù)，阻力、升力和力矩最大平均絕對(duì)誤差分別為5.658 0、3.728 5和0.133 0，誤差均滿足要求，說明當(dāng)輸入尾流信號(hào)穩(wěn)態(tài)值一定，振幅不同時(shí)，基于POD方法建立的尾流激勵(lì)狀態(tài)下葉片氣動(dòng)力降階模型可以用于葉片氣動(dòng)力快速分析，表明基于POD方法的降階模型具有實(shí)用性和可行性。

4 結(jié)論

1) 當(dāng)更換輸入的壓力信號(hào)以后，POD降階模型與CFD計(jì)算結(jié)果誤差很小。

2) 本文建立的POD氣動(dòng)力降階模型能快速準(zhǔn)確的分析上游尾流激勵(lì)下的葉片氣動(dòng)力。

3) 尾流輸入的穩(wěn)態(tài)值一定時(shí)，尾流壓力振幅影響葉片氣動(dòng)力降階模型精度。