基于POD的尾流激勵葉片氣動力降階模型

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

(1.中北大學 機電工程學院， 太原 030051； 2.太原學院數學系， 太原 030001)

保證航空發動機具有較長使用壽命和較強可靠性，一直是現代航空發動機設計過程中追求的目標[1]。在發動機內部動、靜葉之間的高速相對轉動，使得下游葉片極容易受到周期性的上游尾流激勵，誘發下游葉片流場產生周期性振動，進而使得下游葉片表面氣動力具有周期性。這種周期性的氣動力有可能導致葉片強迫振動，甚至出現嚴重的疲勞破壞，影響整個發動機的可靠性和使用壽命。因此，研究尾流激勵下葉片的氣動性能，能為航空發動機葉片的氣動優化提供依據，對提高整個航空發動機的質量、使用壽命和可靠性具有重要意義。

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

為快速計算上游尾流激勵下葉片氣動力，采用POD法提取主要降階基模態，構建基于POD的尾流激勵下葉片氣動力降階模型，并對葉片氣動力數據重構，實現快速分析上游尾流激勵下的葉片氣動力。

1 POD降階基模型建立

(1)

為簡化描述，采用矩陣u來表示流場氣動力：

(2)

根據POD理論，可將尾流激勵下的葉片氣動力分解為平均值和脈動值的疊加[18]：

(3)

(4)

其中，αr為第r個POD基向量對應的模態系數；φr為POD降階基向量；R表示POD基的個數，選擇合適的階次R就可以得到尾流激勵的葉片氣動力降階模型。

圖1 葉片氣動力系統示意圖

1.1 POD降階基提取

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

(5)

(6)

令U={u1，u2，u3，…，uN}，對式(6)求偏導數有：

(7)

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

Zφ=λφ

(8)

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

(9)

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

1.2 氣動力數據重構

對于任一時刻的尾流輸入總壓Pr，響應的POD模態系數為αr，則有：

(10)

2 算例分析

二維葉片流場CFD計算模型如圖2所示，計算網格采用ICEM輸出的非結構網格，網格單元數為164 606，節點數為92 040，葉片網格結構如圖3所示，導入Fluent計算葉片氣動力，流體定義為理想氣體，采用Spallart-Allmaras湍流模型，進口設置為壓力進口pressure-inlet。穩態時，進口總壓為 120 500 Pa，出口壓力101 325 Pa，溫度300 K。瞬態計算時，將UDF編譯的上游尾流加載在流場進口最下方點，尾流速度大小10 m/s，方向為進口y方向。

圖2 葉片CFD計算模型示意圖

圖3 葉片網格結構示意圖

上游尾流輸入的進口總壓波形如圖4所示，設置計算時間t=[0，0.08]，時間步長為2e-5 s，進口壓力隨時間呈現周期性變化，計算可得4 000個進口總壓的葉片氣動力數據并保存。

圖4 周期信號壓力波波形

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

圖5 CFD計算結果絕對誤差曲線

圖6 能量比值隨POD模態變化曲線

降階模型計算結果與CFD計算結果如圖7所示。圖中，紅色實線為CFD計算結果，藍色實線為POD降階模型的計算結果。從圖中可以看出，降階模型計算結果與CFD計算結果一致，因此可以采用POD方法構建尾流激勵氣動力降階模型。同時CFD計算結果與降階模型計算結果之間存在一定的誤差。

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

(11)

(12)

式中：i為第i個計算結果的氣動力數據，Oi為POD降階模型得到的氣動力數據，Fi為CFD求解的氣動力數據。分別將數據代入式(11)和式(12)可以得到系統的均方根誤差和平均絕對誤差，如表1。

圖7 上游尾流激勵下的葉片氣動力響應隨時間的變化曲線

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

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

3 驗證POD降階模型

為了進一步說明降階模型可以計算不同尾流激勵下的葉片氣動力，將圖4中的將尾流壓力波換成一組正弦信號波，波形如圖8所示，葉片計算的穩態邊界條件保持不變，將正弦尾流壓力信號輸入到下游葉片的入口處，完成CFD計算，得到正弦信號激勵下葉片氣動力數據。

將壓力信號數據直接輸入到已經建立的葉片氣動力降階模型，得到葉片氣動力響應值并保存，所得計算結果與CFD計算結果如圖9所示。利用POD計算結果和CFD計算結果，得到降階模型的均方根誤差和平均絕對誤差，如表2。

圖8 正弦信號壓力波波形圖

圖9 游尾流激勵下的葉片氣動力響應隨時間的變化曲線

表1 正弦信號輸入的均方根誤差和平均絕對誤差

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

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

4 結論

1) 當更換輸入的壓力信號以后，POD降階模型與CFD計算結果誤差很小。

2) 本文建立的POD氣動力降階模型能快速準確的分析上游尾流激勵下的葉片氣動力。

3) 尾流輸入的穩態值一定時，尾流壓力振幅影響葉片氣動力降階模型精度。